gamma

Gamma function.

SYNOPSIS:

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

DESCRIPTION:

Returns the gamma function of the argument. The result is correctly signed, and the sign (+1 or -1) is also returned in a global (extern) variable named sgngam. This variable is also filled in by the logarithmic gamma function lgam().

Arguments |x| <= 34 are reduced by recurrence and the function approximated by a rational function of degree 6/7 in the interval (2,3). Large arguments are handled by Stirling's formula. Large negative arguments are made positive using a reflection formula.

ACCURACY:
                      Relative error:
arithmetic   domain     # trials      peak         rms
    DEC      -34, 34      10000       1.3e-16     2.5e-17
    IEEE    -170,-33      20000       2.3e-15     3.3e-16
    IEEE     -33, 33     20000       9.4e-16     2.2e-16
    IEEE      33, 171.6   20000       2.3e-15     3.2e-16

Error for arguments outside the test range will be larger owing to error amplification by the exponential function.