Gamma Function Using Spouge's Method

[lcwright1964]

Well I did find that Maple worksheet. Looks like I basically ported Viktor's C++ code for Godfrey's method and took advantage of Maple's easy to use linear algebra capabilities.

Unfortunately I did NOT expand things out and collect terms as Viktor did (and you used) so my result looks a bit different.

My preliminary conclusion is that going with Pugh's recommendation for the arbitrary parameter one can keep the relative error the same (at least 10 digits good for positive real arguments as long as one computes in 13-digit precision), with one fewer terms in the series.

Here is my result with the coefficients presented to 13 significant digits:

(0.9446965718570e-2 + 1.367185253712/z - 1.839535625335/(z+1) + .6838535775175/(z+2) - 0.6542905498252e-1/(z+3) + 0.6572071053721e-3/(z+4)) * ((z+5.081000210223)/exp(1))^(z-1/2)

Input has been adjust to compute Gamma(z), not Gamma(z+1) or z!.

The scaling by e in the second factor before exponentiation should prevent overflow when other formulas tend to for larger arguments, but I have not rigourously examined that.

In this case the arbitrary parameter recommended by Pugh is 5.58100021022. In the original paper Lanczos does not present an example with n = 5 (so 6 coefficients) for comparison.

In the version you use now you embed 7 12-digit coefficients and one 2-digit parameter--5.5 (from 5 + .5). Here you embed 6 13-digit coefficients and one 13-digit parameter. And there is one less series term. So this should be faster and smaller, at least a little.

I should work on doing what Viktor does so as to produce a polynomial divided by a product, etc. In the meantime, this is what I have.

I work in Mathematica more these days so I should see what happens there.

Les