Gamma Function by Stieltjes Continued Fraction

[Dieter]

I now did some more tests with the Stieltjes algorithm, and I think it is a very nice method that can deliver excellent results without much effort. The resulting accuracy depends on different factors: The number of terms (for our purposes three will usually do, four give even better results), the restriction of integer coefficients (results improve if real values are allowed) and the threshold below which the shift-and-divide method is applied (the higher the better, but a high threshold means more multiplications that may cause some loss of accuracy for small arguments).

Here are some examples coded in VBA. The following examples are not optimized (e.g. there is no provision to prevent overflow), they only are intended to show the basic idea.

First of all a here is a three-term implementation. It uses the simplification that originally led to three integer coefficients 4, 5 and 6.

Code:

Function st3gamma(x)

sqrt2pi = 2.50662827463100

a = 4.024278
b = 5.000135   ' two coefficients
t = 7          ' threshold for shift-and-divide

s = 1
Do While x < t
   s = s * x
   x = x + 1
Loop

q = 1 / 2 / (1 / (b * (x + 1 / (a * x))) + 6 * x)

st3gamma = sqrt2pi * x^(x-0.5) * Exp(q-x) / s

End Function

Two real coefficients, threshold = 7 => error ~ 1,1 E–11.

Set a = 4, b = 5.00003671 and t = 8 => error ~ 1 E–11.
Set a = 4, b = 5.00005017 and t = 9 => error ~ 2 E–12.

With four terms the results get even more accurate. Here the original Stieltjes CF method is used, thus the first Stieltjes coefficients 1/12 and 1/30 show up in the code.

Code:

Function st4gamma(x)

sqrt2pi = 2.50662827463100

a = 2
b = 3.963421   ' two coefficients
t = 7          ' threshold for shift-and-divide

s = 1
Do While x < t
   s = s * x
   x = x + 1
Loop

q = 1 / (12 * (x + 1 / (30 * (x + 1 / (b * (x + 1 / (a * x)))))))

st4gamma = sqrt2pi * x^(x-0.5) * Exp(q-x) / s

End Function

One real coefficient, threshold = 7 => error ~ 7 E–13.
Set a = 2, b = 3.963274 and t = 8 => error ~ 4 E–13.
Set a = 1.957, b = 3.96266 and t = 8 => error ~ 7 E–14.

All this is more accurate than the original Stieltjes continued fraction method with its rational coefficients. The versions above simply trade the extremely high accuracy for large arguments (which is lost) for higher accuracy with smaller arguments down to 7.

Now all we need is a sufficiently accurate calculator that can take advantage of the examples above. ;-) I just tried the second last version with error < 4 E–13 on my trusted 35s and the results I got seem to be off only in the last digit – which is not caused by the algorithm but due to roundoff errors during the calculation.

In an actual implementation there are two crucial points: avoid overflow and evaluate exp(q–x) as accurately as possible. The calculator programs earlier in this thread handle both problems quite well.

Dieter