lgam

Natural logarithm of gamma function.

SYNOPSIS:

double x, y, lgam();
y = lgam(x);

DESCRIPTION:

Returns the base e (2.718...) logarithm of the absolute value of the gamma function of the argument. The sign (+1 or -1) of the gamma function is returned in a global (extern) variable named sgngam.

For arguments greater than 13, the logarithm of the gamma function is approximated by the logarithmic version of Stirling's formula using a polynomial approximation of degree 4. Arguments between -33 and +33 are reduced by recurrence to the interval [2,3] of a rational approximation. The cosecant reflection formula is employed for arguments less than -33.

Arguments greater than MAXLGM return MAXNUM and an error message.

MAXLGM = 2.035093e36 for DEC arithmetic or 2.556348e305 for IEEE arithmetic.

ACCURACY:

arithmetic      domain        # trials     peak         rms
    DEC     0, 3                  7000     5.2e-17     1.3e-17
    DEC     2.718, 2.035e36       5000     3.9e-17     9.9e-18
    IEEE    0, 3                 28000     5.4e-16     1.1e-16
    IEEE    2.718, 2.556e305     40000     3.5e-16     8.3e-17

The error criterion was relative when the function magnitude was greater than one but absolute when it was less than one.
This test used the relative error criterion, though at certain points the relative error could be much higher than indicated.

    IEEE    -200, -4             10000     4.8e-16     1.3e-16