This function returns a random variate from the gamma
distribution. The distribution function is,
\[p(x) dx = {1 \over \Gamma(a) b^a} x^{a-1} e^{-x/b} dx\]
for \(x > 0\).
The gamma distribution with an integer parameter a
is known as
the Erlang distribution.
The variates are computed using the Marsaglia-Tsang fast gamma method.