# Gamma function with Spouge's formula (Mathematica)

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,

$z! = (z+a)^{z+1/2} e^{-(z+a)} \sqrt{2 \pi} \left[ c_0 + \sum_{k=1}^{a-1} \frac{c_k}{z+k} + \epsilon(z) \right],$

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

$c_0 = 1\,$
$c_k = \frac{1}{\sqrt{2 \pi}} \frac{(-1)^{k-1}}{(k-1)!} (-k+a)^{k-1/2} e^{-k+a} \quad k=1,2,\dots, a-1.$

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

$\epsilon \le a^{-1/2} (2 \pi)^{-(a+1/2)}$

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

$a = -\frac{\log \epsilon}{\log 2 \pi} \approx -1.84 \, \log \epsilon.$

If we take $\epsilon = 10^{-n}\,\!$, this becomes

$a = \frac{\log 10}{\log 2 \pi} n \approx 1.26 n.$

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.

hijacker
hijacker
hijacker
hijacker