  A fast and numerically stable algorithm is described for computing the
discrete Hankel transform of order $0$ as well as evaluating Schl\"{o}milch and
Fourier--Bessel expansions in $\mathcal{O}(N(\log N)^2/\log\!\log N)$
operations. The algorithm is based on an asymptotic expansion for Bessel
functions of large arguments, the fast Fourier transform, and the Neumann
addition formula. All the algorithmic parameters are selected from error bounds
to achieve a near-optimal computational cost for any accuracy goal. Numerical
results demonstrate the efficiency of the resulting algorithm.
