  We present the singular Euler--Maclaurin expansion, a new method for the
efficient computation of large singular sums that appear in long-range
interacting systems in condensed matter and quantum physics. In contrast to the
traditional Euler--Maclaurin summation formula, the new method is applicable
also to the product of a differentiable function and a singularity. For
suitable non-singular functions, we show that the approximation error decays
exponentially in the expansion order and polynomially in the characteristic
length scale of the non-singular function, where precise error estimates are
provided. The sum is approximated by an integral plus a differential operator
acting on the non-singular function factor only. The singularity furthermore is
included in a generalisation of the Bernoulli polynomials that form the
coefficients of the differential operator. We demonstrate the numerical
performance of the singular Euler--Maclaurin expansion by applying it to the
computation of the full non-linear long-range forces inside a macroscopic
one-dimensional crystal with $2\times 10^{10}$ particles. A reference
implementation in Mathematica is provided online.
