Gamma Function Using Spouge's Method

[lcwright1964]

Thanks, Dieter, this is excellent--and very compact too.

Now that I am revisiting the gamma/factorial for calculators issue, it is coming back to me why Spouge and Lanczos, though theoretically attractive because they are convergent and accurate for small real arguments, have seemed less attractive for our calculator work. It is that the just don't do better than the old fashioned asymptotic formulae, even if in theory they should. In Spouge especially, on computes coefficients on the fly and it is no more accurate.

I don't think I ever managed to improve on J-M Baillard's GAM+ implementation in the old HP41 program archives. His 70-step program has some key features--it uses synthetic registers when storage is needed, does most everything on the stack, cleverly rearranges the Stirling series to simplify the constants embedded in the code, and uses the calculator's FACT function if the argument is an integer. It does the Pochhammer shift for arguments less than 5, which seems to suffice in a 10-digit environment to give at least eight good digits and usually nine for most anything.

I also rediscovered Peter Luschny's superb page:

http://www.luschny.de/math/factorial/app...Cases.html

I think I was once drawn to this pseudocode (edited by me to include Luschny's suggested simplification in the formula) which implements the Stieltjes continued fraction:

Stieltjes3Factorial(x)
y = x+1; p = 1;
while y < 7 do p = p*y; y = y+1; od;
r = exp(y*(log(y)-1)+1/(2*(6*y+1/(5*(y+1/(4*y))))));
if x < 7 then r = r/p fi;
return r*sqrt(2*Pi / y) end;

However, despite it's compactness and elegance and guarantee of at least 8 digits, I don't think I was ever able to better my favourite Baillard program.

Fun times!

Les