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\times(n-1)!$, and conversely, $(n-1)!=n!/n$, it is possible to define 0! as 1!/1, or 1.

If you plot the factorial of integers on a graph, you'll notice that the points you plotted will be along what appears to be some kind of an exponential curve. In any case, it looks very much possible to connect the points with a smooth line. This begs the question: might it be possible to extend the definition of the factorial and create a smooth function that has a value for all real numbers, not just integers, but for integers, its value is equal to that of the factorial itself?

The answer is a resounding yes; such a function indeed exists. In fact, the so-called Gamma function extends the definition of the factorial to the entire complex plane. In other words, while the factorial is defined only for non-negative integers, the Gamma function is defined for any complex number z as long as z is not a negative integer.

Definition

The Gamma function is defined by the following integral:

\[\Gamma(z)=\int\limits_0^\infty t^{z-1}e^{-t}~dt.\]

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

\[\Gamma(z+1)=z\Gamma(z).\]

Consequently, for positive integers:

\[n!=\Gamma(n+1).\]

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:

\[\Gamma(-z)=\frac{-\pi}{z\Gamma(z)\sin\pi z}.\]

But how do you actually compute the Gamma function? The definition integral is not very useful; to produce an accurate result, an insanely high number of terms would have to be added during some numerical integration procedure. There must be some other way; there indeed is another way, several other ways in fact, each yielding useful algorithms for handheld calculators.

The Numerical Recipes Solution

Back in 1979, when I was possessed with the desire to squeeze a Gamma function implementation into the 72 programming steps of my tiny calculator, I could not find any ready reference on numerical approximations of this function. Eventually, I came across a 1955 Hungarian book on differential and integral calculus, which provided a usable approximation formula. Being theoretical in nature, the book offered not only no hints on efficient computation of the function, it didn't even provide numerical approximations of the coefficients in the formula; that was all left to the reader to work out (or to not work out, as might be the case, since theoretical mathematicians rarely sink to the level of actually working out a numeric result). Work it out I did, making my little calculator run for several days from its wall charger, while it calculated some of those coefficients. Eventually I had my result: a 71-step program that provided 8 or more digits of precision for real $z$'s in the range of $0.5\lt z\lt 1.5$. From these values, it is possible to compute the Gamma function of any real number using the recurrence relation.

As it turns out, while my approximation is not a bad one, there's a much better method for calculating the Gamma function. Numerical Recipes in C (2nd ed. Cambridge University Press, 1992) demonstrates a much more efficient algorithm, the so-called Lanczos approximation for computing the Gamma function for any positive argument (indeed, for any complex argument with a nonnegative real part) to a high level of accuracy. Here is their magic formula, computing the Gamma function with an error that's smaller than $|\epsilon|\lt 2\times 10^{-10}$ for any complex $z$ for which ${\rm Re}~z\gt 0$:

\[\Gamma(z)=\left[\frac{\sqrt{2\pi}}{z}\left(p_0+\sum\limits_{n=1..6}\frac{p_n}{z+n}\right)\right](z+5.5)^{z+0.5}e^{-(z+5.5)},\]

where

\begin{align}p_0&=1.000000000190015,\\
p_1&=76.18009172947146,\\
p_2&=-86.50532032941677,\\
p_3&=24.01409824083091,\\
p_4&=-1.231739572450155,\\
p_5&=1.208650973866179\times 10^{-3},\\
p_6&=-5.395239384953\times 10^{-6}.\end{align}

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

\[p_0+\sum\limits_{n=1..6}\frac{p_n}{z+n}=\frac{p_0(z+1)(z+2)(z+3)(z+4)(z+5)(z+6)+\sum\limits_{n=1..6}\frac{p_n(z+1)(z+2)(z+3)(z+4)(z+5)(z+6)}{z+n}}{(z+1)(z+2)(z+3)(z+4)(z+5)(z+6)}.\]

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 \(\sqrt{2\pi}\). The result is the following formula:

\[\Gamma(z)=\frac{\sum\limits_{n=0..N}q_nz^n}{\prod\limits_{n=0..N}(z+n)}(z+5.5)^{z+0.5}e^{-(z+5.5)},\]

\begin{align}q_0&=75122.6331530,\\
q_1&=80916.6278952,\\
q_2&=36308.2951477,\\
q_3&=8687.24529705,\\
q_4&=1168.92649479,\\
q_5&=83.8676043424,\\
q_6&=2.50662827511.\end{align}

The Lanczos Approximation

Numerical Recipes neglects to mention an important fact: that the Lanczos approximation (A Precision Approximation of the Gamma Function. J. SIAM Numer. Anal. Ser. B, Vol. 1 1964. pp. 86-96) can be used with other choices of values. It is possible to derive simpler versions with fewer coefficients; it is also possible to develop extremely accurate versions. As a matter of fact, the Lanczos approximation can be used to compute the Gamma function to arbitrary precision. I am grateful to Paul Godfrey for providing me with the details that are presented here.

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

\[\ln\Gamma(z+1)=\ln{\bf\mathrm{ZP}}+(z+0.5)\ln(z+g+0.5)-(z+g+0.5),\]

where $\vec{Z}$ is an $n$-dimensional row vector constructed from the function argument, $z$:

\[{\bf\mathrm{Z}}=\left[1~~\frac{1}{z+1}~~\frac{1}{z+2}~...~\frac{1}{z+n-1}\right],\]

and $\vec{P}$ is an $n$-dimensional column vector constructed as the product of several $n\times n$ matrices and the $n$-dimensional column vector $\vec{F}$:

\[{\bf\mathrm{P}}={\bf\mathrm{D}}\cdot{\bf\mathrm{B}}\cdot{\bf\mathrm{C}}\cdot{\bf\mathrm{F}},\]

where

\[{\bf\mathrm{B}}_{ij}=\left\{\begin{matrix}1&{\rm if}~i=0,\\-1^{j-i}\begin{pmatrix}i+j-1\\j-i\end{pmatrix}&{\rm if}~i>0,~j\ge i,\\0&{\rm otherwise}.\end{matrix}\right.\]

$\vec{C}$ is a matrix containing coefficients from Chebyshev polynomials:

\[{\bf\mathrm{C}}_{ij}=\left\{\begin{matrix}1/2&{\rm if}~i=j=0,\\0&{\rm if}~j>i,\\-1^{i-j}\sum\limits_{k=0}^{i}\begin{pmatrix}2i\\2k\end{pmatrix}\begin{pmatrix}k\\k+j-i\end{pmatrix}&{\rm otherwise},\end{matrix}\right.\]

and

\[{\bf\mathrm{D}}_{ij}=\left\{\begin{matrix}0&{\rm if}~i\ne j\\1&{\rm if}~i=j=0,\\-1&{\rm if}~i=j=1,\\\frac{{\bf{\mathrm{D}}_{i-1,j-1}}2(2i-1)}{i-1}&{\rm otherwise},\end{matrix}\right.\]

\[{\bf\mathrm{F}}_{i}=\frac{(2i)!e^{i+g+0.5}}{i!2^{2i-1}(i+g+0.5)^{i+0.5}},\]

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:

\[{\bf\mathrm{E}}={\bf\mathrm{C}}\cdot{\bf\mathrm{F}},\]

\[|\epsilon|=\left|\frac{\pi}{2\sqrt{2e}}\left(e^g\sqrt{\pi}-\sum\limits_{i=0}^{n-1}-1^{i}{\bf\mathrm{E}}_i\right)\right|.\]

I endeavored to follow up on Godfrey's work and wrote an arbitrary precision implementation of the Lanczos approximation in the C++ language, utilizing the GMP multi-precision library and my custom C++ wrapper classes for the functions contained therein. This arbitrary precision implementation enabled me to compute, for instance, the Gamma function of 101 (i.e., 100!) to full integral precision, i.e., to a precision better than 158 significant digits.

On a more mundane level, my implementation also allowed me to experiment with various sets of coefficients, trying to find the most accurate and economical version of the Lanczos approximation for pocket calculators. I indeed found a set of only 4 coefficients that has a computed precision of $|\epsilon|\lt 2\times 10^{-7}$, but in practice, appears to yield more than 8 significant digits for all positive real arguments. The precision is further increased when 5 or 6 coefficients are used:

The simplicity of the formula for $n=4$ makes it suitable for use with most programmable calculators, including some of the simplest models; and that includes my old OEM Commodore machine. Which means the end of my quest after more than 20 years of programming: I finally have a "perfect" method of computing the Gamma function, to the maximum precision possible, on my first programmable calculator.

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:

\[x!=x^xe^{-x}\sqrt{2\pi x}.\]

Unfortunately, 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}\sqrt{2\pi}\)).

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:

\[x!=x^xe^{-x}\sqrt{2\pi x}\left(1+\frac{1}{12x}\right).\]

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:

\[\ln\Gamma(z)=\left(z-\frac{1}{2}\right)\ln z-z+\frac{\ln 2\pi}{2}+\frac{B_2}{1\cdot 2}\frac{1}{z}+\frac{B_4}{3\cdot 4}\frac{1}{z^3}+...+\frac{B_{2k}}{(2k-1)2k}\frac{1}{z^{2k-1}}+...\]

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):

\begin{align}1+2B_1&=0,\\ 1+3B_1+3B_2&=0,\\ 1+4B_1+6B_2+4B_3&=0,\\ 1+5B_1+10B_2+10B_3+5B_4&=0,\end{align}

and so on. The first few $B_k$'s are: $B_1=-1/2$, $B_2=1/6$, $B_3=0$, $B_4=-1/30$, $B_5=0$, $B_6=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_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\times 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 ${\rm 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:

\[\ln\Gamma(z)\simeq z\ln z-z+\ln\sqrt{\frac{2\pi}{z}}+\frac{1}{1188z^9}-\frac{1}{1680z^7}+\frac{1}{1260z^5}-\frac{1}{360z^3}+\frac{1}{12z}.\]

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:

\[\ln\Gamma(x)-\left[\frac{\ln 2\pi}{2}+\left(x-\frac{1}{2}\right)\ln x-x\right]=\frac{1}{12x}-\frac{1}{360x^3}+\frac{1}{1260x^5}-\frac{1}{1680x^7}+...,\]

which he then rewrote as follows:

\[\ln\left[\left(\frac{\Gamma(x)}{\sqrt{\frac{2\pi}{x}}}\right)^{1/x}\frac{e}{x}\right]=\frac{1}{12x^2}-\frac{1}{360x^4}+\frac{1}{1260x^6}-\frac{1}{1680x^8}+...=\ln M(x).\]

Then, through a power series expansion, he obtained:

\[e^{2\ln M(x)}=1+\frac{1}{6x^2}+\frac{1}{120x^4}+\frac{13}{9072x^6}+...\]

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

\[x\sinh \frac{1}{x}=1+\frac{1}{6x^2}+\frac{1}{120x^4}+\frac{1}{5040x^6}+...\]

This led to an approximation formula for the Gamma function that is fairly accurate, and requires no stored constants except for \(\sqrt{2\pi}\) (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):

\[\left(\frac{x}{e}\sqrt{x\sinh\frac{1}{x}\left[+\frac{1}{810x^6}\right]}\right)^x\sqrt{\frac{2\pi}{x}}\simeq\Gamma(x).\]

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.

And Another!

In early July, 2006, I received an e-mail from Hungary, from Gergő Nemes, who sent me several formulas for the Gamma function that he developed. Two of these caught my attention especially, because they were surprisingly accurate yet simple. Then, a few days later, Gergő wrote to me again, this time with a formula that approximates the Gamma function with an error term proportional to $n^{-8}$, which is quite remarkable! Without further ado, here are his more interesting formulae (updated November 16, 2006):

1. $n!=e^{-n}\sqrt{2\pi n}(n+1/12n+1/1440n^3+239/362880n^5)^n$,
2. $n!=n^n\sqrt{2\pi n}\exp(1/(12n+2/5n)-n)(1+{\cal O}(n^{-5}))$,
3. $n!=n^n\sqrt{2\pi n}\exp(1/[12n+2/(5n+53/42n)]-n)(1+{\cal O}(n^{-8}))$.

Gergő says that he was able to prove that his formulae are correct. I have not seen his proofs, but I did check the formulae for selected values and confirmed that they do indeed appear to behave as promised. (Since our initial exchange of e-mails, Gergő has informed me that he has identified the second and third of these three formulae as the initial approximations of the Gamma function using continued fractions.)

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:

\[\gamma(a,x)=\int\limits_0^x t^{a-1}e^{-t}dt.\]

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

\[\gamma(a,x)=x^ae^{-x}\sum\limits_{n=0}^\infty\frac{x^n}{a(a+1)...(a+n)}.\]

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\lt 50$, $x=30$, 8-digit accuracy is achieved.