  An exact expression for the leading-order (LO) gluon distribution function
$G(x,Q^2)=xg(x,Q^2)$ from the DGLAP evolution equation for the proton structure
function $F_2^{\gamma p}(x,Q^2)$ for deep inelastic $\gamma^* p$ scattering has
recently been obtained [M. M. Block, L. Durand and D. W. McKay, Phys. Rev.
D{\bf 79}, 014031, (2009)] for massless quarks, using Laplace transformation
techniques. Here, we develop a fast and accurate numerical inverse Laplace
transformation algorithm, required to invert the Laplace transforms needed to
evaluate $G(x,Q^2)$, and compare it to the exact solution. We obtain accuracies
of less than 1 part in 1000 over the entire $x$ and $Q^2$ spectrum. Since no
analytic Laplace inversion is possible for next-to-leading order (NLO) and
higher orders, this numerical algorithm will enable one to obtain accurate NLO
(and NNLO) gluon distributions, using only experimental measurements of
$F_2^{\gamma p}(x,Q^2)$.
