  Hadron production at low transverse momenta in semi-inclusive deep inelastic
scattering can be described by transverse momentum dependent (TMD)
factorization. This formalism has also been widely used to study the Drell-Yan
process and back-to-back hadron pair production in $e^+e^-$ collisions. These
processes are the main ones for extractions of TMD parton distribution
functions and TMD fragmentation functions, which encode important information
about nucleon structure and hadronization. One of the most widely used TMD
factorization formalism in phenomenology formulates TMD observables in
coordinate $b_\perp$-space, the conjugate space of the transverse momentum. The
Fourier transform from $b_\perp$-space back into transverse momentum space is
sufficiently complicated due to oscillatory integrands that it requires a
careful and computationally intensive numerical treatment in order to avoid
potentially large numerical errors. Within the TMD formalism, the azimuthal
angular dependence is analytically integrated and the two-dimensional $b_\perp$
integration reduces to a one-dimensional integration over the magnitude
$b_\perp$. In this paper we develop a fast numerical Hankel transform algorithm
for such a $b_\perp$-integration that improves the numerical accuracy of TMD
calculations in all standard processes. Libraries for this algorithm are
implemented in Python 2.7 and 3, C++, as well as FORTRAN77. All packages are
made available open source.
