**Calculators and the Gamma Function**

No wonder
mathematicians find numbers to be the passion of a lifetime. Deceptively simple
things can lead to such amazing complexity, to intriguing links between
seemingly unconnected concepts. Take the simple concept of the
factorial. The factorial of a positive integer, *n*, is defined as the
product of all the integers between 1 and *n* inclusive. Thus, for
instance, the factorial of 5 is 5!=1×2×3×4×5, or 120.
Because *n*!=*n*×(*n*-1)!, and conversely, (*n*-1)!=*n*!/*n*,
it is possible to define 0! as 1!/1, or 1.

The Gamma function is defined by the following
integral:

For all complex numbers *z*, the following
recurrence relation is true:

Consequently, for positive integers:

Some very helpful
relationships exist between the Gamma function values for various arguments.
For instance, the following relationship makes it possible to compute the Gamma
function for a negative argument easily:

**The Numerical
Recipes Solution**

where:

*p*_{0} = 1.000000000190015; *p*_{1} = 76.18009172947146; *p*_{2} = -86.50532032941677; *p*_{3} = 24.01409824083091

*p*_{4} = -1.231739572450155 ; *p*_{5} = 1.208650973866179 · 10^{-3}; *p*_{6} = -5.395239384953 · 10^{-6}^{}

This formula can be rearranged in a form more suitable
for the limited program capacity of programmable calculators using the
following identity:

This way after like terms are collected, the numerator
becomes a simple polynomial, and the denominator can be calculated using a
simple loop. The calculation can be further simplified by multiplying each
constant with .
The result is the following formula:

*q*_{0} = 5122.6331530; *q*_{1} = 80916.6278952; *q*_{2}= 6308.2951477; *q*_{3}= 8687.24529705

*q*_{4}= 1168.92649479; *q*_{5}= 83.8676043424; *q*_{6}= 2.50662827511

The *Lanczos approximation* can, in essence, be
reduced to a simple vector inner product and then some additional elementary
computations:

Where **Z** is an *n*-dimensional row vector
constructed from the function argument, z:

And **P** is an *n*-dimensional column vector
constructed as the product of several *n*×*n* matrices and the *n*-dimensional
column vector **F**: **P** = **D**·**B**·**C**·**F**

Where

**C** is a matrix
containing coefficients from Chebyshev polynomials:

and

and *g* is an arbitrary real parameter, *n* is
an integer. For suitable choices of *g* and *n*, very good precision
can be obtained. The typical error can be computed as follows: **E** = **C**·**F**

When great accuracy isn't required, a much simpler
approximation exists for the generalized factorial: Stirling's formula. In its
original form, this is what this formula looks like:

This formula isn't particularly accurate (a slightly more
accurate version, often called Burnside's formula, gets rid of the *x* under
the square root and replaces all remaining *x*'s with *x*+1/2 in the
right-hand expression: *x*!=(*x*+1/2)^{x}^{+1/2}*e ^{
-x-1/2}*).

A simple modification, however, can make Stirling's
formula accurate to three or four significant digits for most arguments, making
it suitable for many applications:

These variants may make you wonder: are there more
accurate versions of Stirling's formula? Indeed there are. Here is a
generalized version for the logarithm of the Gamma function:

In this formula, the coefficients *B** _{n}* are derived from the so-called Bernoulli
polynomials. They can easily be calculated from the following equations (in
which you'll no doubt recognize those good old ordinary binomial coefficients):

1 + 2*B*_{1} = 0

1 + 3*B*_{1} + 3*B*_{2} = 0

1 + 4*B*_{1} + 6*B*_{2} + 4*B*_{3} = 0

1 + 5*B*_{1} + 10*B*_{2} + 10*B*_{3} + 5*B*_{4} = 0

...

And so on. The first few *B** _{k}*'s are:

Stirling's formula is an odd bird. What you see above may
suggest that the absolute values of *B*'s get ever smaller as the index
increases, and that therefore, Stirling's formula converges. This is not so;
after *B*_{6} in fact
the absolute values of these *B*'s begins to increase (e.g., *B*_{8}=-1/30) . After a while, the *B*'s will
get very large, and subsequent sums in the approximation will swing around ever
more wildly. For instance, *B*_{98} is approximately 1.1318×10^{76}! Still, Stirling's formula can be used to
approximate the Gamma function quite nicely if you don't sum too many terms of
the series. In fact, the higher the real part of *z* is, the more accurate
the formula becomes... so one trick in the case when Re* z* is small is to
evaluate Stirling's formula for *z*+*n*, and then divide the result
by *z*(*z*+1)(*z*+2)...(*z*+*n*-1).

One particular set of coefficients is easy to remember,
produces a relatively compact version of Stirling's formula, and yields very
good accuracy for arguments over 5:

While the *Lanczos approximation* may be somewhat
more accurate and compact, it requires several floating point constants that
may be difficult to remember, and take up precious (program or registry) memory
on your calculator. So in many cases, using Stirling's formula may be
preferable. However, it should not be forgotten that even with corrections, it
remains inaccurate for small arguments.

**A Curious
Result**

In November 2002, I received an e-mail from Robert H.
Windschitl who noticed a curious coincidence. He wrote up the extended
Stirling's formula as this:

which he then rewrote as follows:

then, through a power series expansion, he obtained:

At this point, he noticed that this expansion is very
similar to a known function expansion:

This led to an approximation formula for the Gamma
function that is fairly accurate, and requires no stored constants except for (easily
computed on most calculators) and, optionally, the integer value 810 (the part
below that's within square brackets can be omitted at a small cost in terms of
precision loss):

For values greater than 8, this approximation formula
yields 8+ digits of precision even without the correction term that is enclosed
in square brackets. On calculators with limited program or register memory,
this method may make it possible to implement the Gamma function even when
other methods fail.

I have not previously seen this approximation in
literature, although Windschitl modestly suggests that it's only a question of
searching in the right places.

A close relative to the Gamma function is the incomplete
Gamma function. Its name is due to the fact that it is defined with the same
integral expression as the Gamma function, but the infinite integration limit
is replaced by a finite number:

The incomplete Gamma function can be approximated fairly rapidly
using an iterative method:

The magnitude of successive terms quickly diminishes as *n*
increases; 8-12 digit accuracy is achieved for most values of *x*and *a*
after no more than a few hundred terms (often less than a hundred) are added. The incomplete Gamma function can also be used to approximate
the value of the Gamma function itself, by selecting a suitably high
integration limit *x*. For *a*<50, *x*=30, 8-digit accuracy
is achieved.