  This paper describes a trapezoidal quadrature method for the discretization
of singular and hypersingular boundary integral operators (BIOs) that arise in
solving boundary value problems for elliptic partial differential equations.
The quadrature is based on a uniform grid in parameter space coupled with the
standard punctured Trapezoidal rule. A key observation is that the error
incurred by the singularity in the kernel can be expressed exactly using
generalized Euler-Maclaurin formulae that involve the Riemann zeta function in
2D and the Epstein zeta functions in 3D. These expansions are exploited to
correct the errors via local stencils at the singular point using a novel
systematic moment-fitting approach. This new method provides a unified
treatment of all common BIOs (Laplace, Helmholtz, Stokes, etc.). We present
numerical examples that show convergence of up to 32nd-order in 2D and
9th-order in 3D with respect to the mesh size.
