Summary: This paper presents a theoretical discussion as well as novel solution algorithms for problems of scattering on smooth two-dimensional domains under Zaremba boundary conditions, for which Dirichlet and Neumann conditions are specified on various portions of the domain boundary. The theoretical basis of the proposed numerical methods, which is provided for the first time in the present contribution, concerns detailed information about the singularity structure of solutions of the Helmholtz operator under boundary conditions of Zaremba type. The new numerical method is based on the use of Green functions and integral equations, and it relies on the Fourier continuation method for regularization of all smooth-domain Zaremba singularities as well as newly derived quadrature rules which give rise to high-order convergence, even around Zaremba singular points. As demonstrated in this paper, the resulting algorithms enjoy high-order convergence, and they can be used to efficiently solve challenging Helmholtz boundary value problems and Laplace eigenvalue problems with high-order accuracy.