  We extend our previous derivation of 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 to include the effects of
heavy-quark masses. We derive the equation for $G(x,Q^2)$ in two different
ways, first using our original differential-equation method, and then using a
new method based on Laplace transforms. The results do not require the use of
the gluon evolution equation, or, to good approximation, knowledge of the
individual quark distributions. Given an analytic expression that successfully
reproduces the known experimental data for $F_2^{\gamma p}(x,Q^2)$ in a domain
${\cal D}(x,Q^2)$--where $x_{\rm min}(Q^2) \le x \le x_{\rm max}(Q^2)$,
$Q^2_{\rm min}\le Q^2\le Q^2_{\rm max}$ of the Bjorken variable $x$ and the
virtuality $Q^2$--$G(x,Q^2)$ is uniquely determined in the same domain. As an
application of the method, we construct a new global parametrization of the
complete set of ZEUS data on $F_2^{\gamma p}(x,Q^2)$, and use this to determine
the 5 quark gluon distribution, $G(x,Q^2)$, for massless $u, d, s$ and massive
$c, b$ quarks and discuss the mass effects evident in the result. We compare
these results to the gluon distributions for CTEQ6L, and in the domain ${\cal
D}(x,Q^2)$ where they should agree, they do not; the discrepancy is due to the
fact that the CTEQ6L results do not give an accurate description of the ZEUS
$F_2^{\gamma p}(x,Q^2)$ experimental data. We emphasize that our method for
obtaining the LO gluon distribution connects $G(x,Q^2)$ {\em directly} to the
proton structure function without either the need for individual parton
distributions or the gluon evolution equation.
