Browse Prior Art Database

# Fast and Accurate Computation of the Gamma Function

IP.com Disclosure Number: IPCOM000116572D
Original Publication Date: 1995-Oct-01
Included in the Prior Art Database: 2005-Mar-30
Document File: 2 page(s) / 40K

IBM

## Related People

Slishman, GR: AUTHOR

## Abstract

A fast and accurate method to compute the gamma function over real arguments is disclosed. Use the following definitions: 1. relative error is defined as ( x' - x ) / x, where x' is an approximation of x. 2. roundoff error is defined as x' - x, where x' is the closest machine number to a real number x. (Larger x has larger roundoff error, as the value of each machine fraction bit increases with x.) 3. lgamma(x) equals log(abs(gamma(x))) for real x. 4. exp(x) equals e to the power x.

This text was extracted from an ASCII text file.
This is the abbreviated version, containing approximately 75% of the total text.

Fast and Accurate Computation of the Gamma Function

A fast and accurate method to compute the gamma function over
real arguments is disclosed.  Use the following definitions:
1.  relative error is defined as ( x' - x ) / x, where x' is an
approximation of x.
2.  roundoff error is defined as x' - x, where x' is the closest
machine number to a real number x.  (Larger x has larger
roundoff error, as the value of each machine fraction bit
increases with x.)
3.  lgamma(x) equals log(abs(gamma(x))) for real x.
4.  exp(x) equals e to the power x.

The standard algorithm to compute gamma of real x computes
exp(lgamma(x)) and then affixes the appropriate sign.  Since
dexp(y)/exp(y) = dy, the composite  relative error bound is
commensurate with the  roundoff error from y=lgamma(x).  As larger x
admits larger lgamma(x)  roundoff error,  the larger lgamma(x)
roundoff error admits larger  relative error
in exp(lgamma(x)).

To compute gamma with better accuracy and speed, the disclosed
method incorporates a logarithm computation which is accurate to
working precision plus a certain number of extra bits.  After
computing lgamma(x) in-line to extra precision as log_hi + log_lo,
the method then computes accurate gamma(x) as exp(log_hi + log_lo),
which is well-approximated by exp(log_hi) * (1 + log_lo).  The
in-line extra-precision logarithm executes as fast as a
working-precision
l...