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