  The use of integral equation methods for the efficient numerical solution of
PDE boundary value problems requires two main tools: quadrature rules for the
evaluation of layer potential integral operators with singular kernels, and
fast algorithms for solving the resulting dense linear systems. Classically,
these tools were developed separately. In this work, we present a unified
numerical scheme based on coupling Quadrature by Expansion, a recent quadrature
method, to a customized Fast Multipole Method (FMM) for the Helmholtz equation
in two dimensions. The method allows the evaluation of layer potentials in
linear-time complexity, anywhere in space, with a uniform, user-chosen level of
accuracy as a black-box computational method.
  Providing this capability requires geometric and algorithmic considerations
beyond the needs of standard FMMs as well as careful consideration of the
accuracy of multipole translations. We illustrate the speed and accuracy of our
method with various numerical examples.
  Keywords: Layer Potentials; Singular Integrals; Quadrature; High-order
accuracy; Integral equations; Helmholtz equation; Fast multipole method.
