Gamma function with Spouge's formula (Mathematica)

From LiteratePrograms

Jump to: navigation, search

We will implement an arbitrary-precision approximation for the gamma function due to John L. Spouge. The main formula has the same general form as Lanczos' approximation,

where a is an arbitrary positive parameter that controls the accuracy (here taken to be an integer). However, Spouge develops a different set of coefficients ck, given by

This series converges more slowly than Lanczos', but has two important advantages: the coefficients are much easier to calculate, and the error is easier to estimate. It fact, it is easy to choose the parameter a to obtain any desired accuracy.


Choice of parameter

Spouge gives the bound

for the relative error, under the assumption that Re(z) > 0 and a > 2. Conservatively omitting the factor a − 1 / 2 and rearranging gives us the formula

If we take , this becomes

We may accordingly define a function to calculate the parameter a, given a desired number of digits of accuracy:

<<parameter calculation>>=
parameter[n_] := Ceiling[n Log[10]/Log[2 Pi]]

Series evaluation

The formula for the coefficients can be expressed directly in Mathematica:

<<generation of coefficients>>=
c[0] := 1
c[k_] := (2 Pi)^(-1/2) (-1)^(k-1)/(k-1)! (-k+a)^(k-1/2) Exp[-k + a]

The same applies for the series:

<<series evaluation>>=
ser[z_] := c[0] + Sum[c[k]/(z + k), {k, 1, a-1}]

Main function

The following function will evaluate the gamma function of an argument with positive real part. Note that the argument must be shifted by -1 since z! = Γ(z+1).

<<positive case>>=
gammapos[z_] := (z-1+a)^(z-1+1/2) Exp[-z+1-a] Sqrt[2 Pi] ser[z-1]

To make the formula valid in the negative complex half-plane, we use the reflection formula:

<<main function>>=
gamma[z_] := If[Re[z] > 0, gammapos[z], Pi/(Sin[Pi z] gamma[1 - z])]

Wrapping it up, we have:

parameter calculation
generation of coefficients
series evaluation
positive case
main function


We may calculate the relative error against Mathematica's built-in Gamma for a few different arguments:

cases = {1, 2, 1/2, 5037/2793, 5, 123, 4+3 I, -6/7, -13+(17/19)I};
maxerror[A_] := CompoundExpression[a=parameter[A];
  Max[Map[Function[x, Abs[(N[gamma[x], 2 A]-Gamma[x])/Gamma[x]]], cases]]];

These should return, in respective order, 2.5×10-15, 8.2×10-29, 1.6×10-52, and 1.8×10-106. There is a fairly big margin; some speed could be gained by tightening the estimate in the parameter function based on empirical data. Nevertheless, the routine is fairly fast as-is.

In fact, it appears to be faster than Mathematica's built-in gamma function for some inputs. Taking

Timing[N[Gamma[1/4], 1000]]
Timing[N[gamma[1/4], 1000]]

gave 4.4 seconds for Mathematica's Gamma and 3.1 seconds for the implementation presented here. Even better speed, especially for multiple evaluations, could be attained by storing the coefficients numerically.


  • Spouge, John L. "Computation of the gamma, digamma, and trigamma functions", SIAM Journal on Numerical Analysis 31 (1994), no. 3, 931-944.
Download code
Personal tools