  Boundary integral equations are an efficient and accurate tool for the
numerical solution of elliptic boundary value problems. The solution is
expressed as a layer potential; however, the error in its evaluation grows
large near the boundary if a fixed quadrature rule is used. Firstly, we analyze
this error for Laplace's equation with analytic density and the global periodic
trapezoid rule, and find an intimate connection to the complexification of the
boundary parametrization. Our main result is then a simple and efficient scheme
for accurate evaluation up to the boundary for single- and double-layer
potentials for the Laplace and Helmholtz equations, using surrogate local
expansions about centers placed near the boundary. The scheme---which also
underlies the recent QBX Nystr\"om quadrature---is asymptotically exponentially
convergent (we prove this in the analytic Laplace case), requires no
adaptivity, generalizes simply to three dimensions, and has O(N) complexity
when executed via a locally-corrected fast multipole sum. We give an example of
high-frequency scattering from an obstacle with perimeter 700 wavelengths long,
evaluating the solution at $2\times 10^5$ points near the boundary with
11-digit accuracy in 30 seconds in MATLAB on a single CPU core.
