Gamma Function Using Spouge's Method

[lcwright1964]

Thanks for all this. I have to admit I have obviously been out of touch. The last time I paid attention to SandMath I think it consisted of just lower and upper halves, and that's about it--basically, what is still included with Diego's most recent ClonixConfig software. I took a peak at the most recent offerings on TOS, including the simpler SandMath-IV module that doesn't use Library #4, and modfile.exe tells me that there is a lot more going on with these things in the three or four years since.

And thanks for the configuration suggestion. I would probably just burn the SandMath 4x4 constituents with Library #4, and leave SandMatrix and the others be for now.

Of course, you have tempted me to order another of the Clonix family from Diego

Getting back to the original topic, I have used that very same Maple worksheet and Pugh's recommendations to give this n = 4 (i.e. 5 coefficient) version for a 10-digit FOCAL environment:

(0.3264892413e-1 + 1.188688981/z - 1.068410309/(z+1) + .1887119902/(z+2) - 0.2745021709e-2/(z+3)) * ((z+ 3.840881909)/exp(1))^(z-1/2)

In an ideal situation (i.e., a couple of guard digits kicking around) this is supposed to give between 9.5 and 11 edd in the z = 0 to 70 range, but of course we don't have guard digits so 8 to 10 is about the best I think it would do. Given this I think JMB's implementation of the Stirling series in GAM+, with its simpler coefficients, will still do the same or better. That said, I would like to see how this fares in FOCAL code for the HP41. It would even work for the 67/97 where one could store the program on one side of the card (it would be very short) and the six 10-digit constants as register data on the other side . Without Pochhammer shifting for small arguments this should be faster than a Stirling approach, and of course the Spouge approach uses more terms and calculates more.

If I actually turn this into program steps for the 41 or the 67/97 I will post.

Les

P.S. I notice that the coefficients (q0 to q6) from Viktor's page are only presented to 12 digits, not 13. Is that what you are using, or did he or you recompute 13-digit versions for use in your MCODE? That does seem possible given that the original Lanczos coefficients, p0 to p6, are given just above them to up to 16 digits.