Gamma Function Using Spouge's Method

[lcwright1964]

Namir,

I looked at the original Spouge paper from 1994 (attached). There is heavy theoretical math in the proofs, and I don't follow it all, but the max relative error is give by a fairly simple formula in Theorem 1.3.1 at the top of page 934. Basically, the larger a is, the smaller the relative error, and the number of terms N = Ceiling(a) - 1. You use a =12.5, but a = 13 would do better (in theory) than your choice of 12.5, and you still only do the 12 loops.

To compare, Spouge's relative error for a = 12.5 is about 1.2e-11, whereas it is about 4.7e-12 for a =13.

That said it may not make a lot of difference in HP41 or the HP67. I understand that once you start wracking up the terms the promised theoretical precision can get eclipsed by rounding error if one doesn't have a lot of guard digits on hand. I know this is a consideration in arbitrary precision environments (e.g., I read someplace that if you want 45 digits accurate you need to work with 70), but maybe the impact is less with the 13 internal digits of our vintage machines.

I don't know if you made your choice of a based on the original paper or another source, but I thought this was interesting.

Les