  We present a novel integral-equation algorithm for evaluation of Zaremba
eigenvalues and eigenfunctions}, that is, eigenvalues and eigenfunctions of the
Laplace operator with mixed Dirichlet-Neumann boundary conditions; of course,
(slight modifications of) our algorithms are also applicable to the pure
Dirichlet and Neumann eigenproblems. Expressing the eigenfunctions by means of
an ansatz based on the single layer boundary operator, the Zaremba eigenproblem
is transformed into a nonlinear equation for the eigenvalue $\mu$. For smooth
domains the singular structure at Dirichlet-Neumann junctions is incorporated
as part of our corresponding numerical algorithm---which otherwise relies on
use of the cosine change of variables, trigonometric polynomials and, to avoid
the Gibbs phenomenon that would arise from the solution singularities, the
Fourier Continuation method (FC). The resulting numerical algorithm converges
with high order accuracy without recourse to use of meshes finer than those
resulting from the cosine transformation. For non-smooth (Lipschitz) domains,
in turn, an alternative algorithm is presented which achieves high-order
accuracy on the basis of graded meshes. In either case, smooth or Lipschitz
boundary, eigenvalues are evaluated by searching for zero minimal singular
values of a suitably stabilized discrete version of the single layer operator
mentioned above. (The stabilization technique is used to enable robust
non-local zero searches.) The resulting methods, which are fast and highly
accurate for high- and low-frequencies alike, can solve extremely challenging
two-dimensional Dirichlet, Neumann and Zaremba eigenproblems with high
accuracies in short computing times---enabling, in particular, evaluation of
thousands of eigenvalues and corresponding eigenfunctions for a given smooth or
non-smooth geometry with nearly full double-precision accuracy
