# Gamma function with Spouge's formula (Mathematica)

### From LiteratePrograms

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 *c*_{k}, 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.

## Contents |

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

<<gammaspouge.nb>>=parameter calculation generation of coefficients series evaluation positive case main function

## Testing

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]]]; maxerror[10] maxerror[20] maxerror[40] maxerror[80]

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.

## References

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

Download code |