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.

Definition

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₀ = 1.000000000190015; p₁ = 76.18009172947146; p₂ = -86.50532032941677; p₃ = 24.01409824083091

p₄ =  -1.231739572450155 ; p₅ = 1.208650973866179 · 10^-3; p₆ = -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₀ = 5122.6331530; q₁ = 80916.6278952; q₂= 6308.2951477; q₃= 8687.24529705

q₄= 1168.92649479; q₅= 83.8676043424; q₆= 2.50662827511

The Lanczos Approximation

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

Stirling's Formula

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/2e ^-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 + 2B₁ = 0
1 + 3B₁ + 3B₂ = 0
1 + 4B₁ + 6B₂ + 4B₃ = 0
1 + 5B₁ + 10B₂ + 10B₃ + 5B₄ = 0
...

And so on. The first few B_k's are: B₁=-1/2, B₂=1/6, B₃=0, B₄=-1/30, B₅=0, B₆=1/42,...

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₆ in fact the absolute values of these B's begins to increase (e.g., B₈=-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₉₈ is approximately 1.1318×10⁷⁶! 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.

The Incomplete Gamma Function

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 xand 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.