🚀 go-pugleaf

It is safer to compute the natural logarithm of the Gamma density first,
because the Gamma function T(a) can easily overflow for moderate large value
of a. Writing g(x;a,b) = (x/b)**(a-1) exp(-x/b) / b / T(a), the logarithm of
this is ln(g(x;a,b)) = (a-1) * ln(x/b) - x/b - ln(b) - lnT(a). This formula
can be further optimized depending on how frequent x, a, and b will take new
values. Anyway, there is a site in Japan
http://momonga.t.u-tokyo.ac.jp/~ooura/
which provides a Fortran 90 codes to compute the logarithm of T(a). This is
very fast (because it uses "classical" Chebyshev approximation) and the
routine offers up to 14 to 16 significant digits.

Leonard Howell wrote:

> Hi all --
>
> I need to compute the gamma probability density function  g(x;a,b)=C(a,b)
> x**(a-1) exp(-x/b) on the order of 10**10 calls (on many different
> occasions) and thus would like to use the FASTEST algorithm available.  It
> must also be double precision with 14 or 15 decimal places of accuracy.
> Advice/coded FORTRAN solutions greatly appreciated.  Other language
> solutions also welcome, if the speed and accuracy is achieved.
>
> Thanks, Leonard