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,

$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 $c k$ , 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.

1 Choice of parameter
2 Series evaluation
3 Main function
4 Testing
5 References

[edit]

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]]

[edit]

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}]

[edit]

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

[edit]

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.

[edit]