Gamma Function Using Spouge's Method

[Dieter]

(08-12-2015 04:03 AM)lcwright1964 Wrote:  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.

I assume that's the calculated upper bound of the error.
But things are different on real calculators. ;-)

(08-12-2015 04:03 AM)lcwright1964 Wrote:  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.

That's exactly the crucial point. The theoretical error that would be possible if (!) all calculations were done with infinite (or at least sufficient) precision is much lower than what we can expect from a real calculator with, here, 10 digit working precision. The three additional guard digits you mentioned are only used internally, so they do not matter here. On the '41 or '67, all we got are ten digits.

Now let's take a look at a few different values of the variable "a" and see what we get:

Mathematically, a value as low as a=9 is sufficient for a 10-digit result over the whole (positive) 41/67 working range, i.e. for input between 0 and approx. 71 (overflow limit). Starting at x=0 the error increases with x, reaches a maximum and then decreases again as x gets larger. The largest relative error is 5,33 E–11 near x=32,8. Which is fine for 10 digits (maybe +1 ULP).

However, on a real '41 or '67 the last two or even three digits are off. Increasing the value of "a" does not neccessarily help here. For instance, a=13 would be good for a largest relative error of 2,72 E–15 (near x=77,9). But in real life the error gets even larger (!) than for a=9.

A value of a=15 theoretically yields an error less than 4 E–17, i.e. full 15–16 digit IEEE double precision. I tried this in Excel VBA as well as on a WP34s in SP mode and got about 10 valid digits near the error maximum.

So indeed the working precision has to be greater than the desired accuracy of the result. I assume that 13 digits would suffice for a 10 digit result. But this cannot be accomplished in a simple user program. On the '41, a 13-digit MCode implementation would be possible.

As long as no additional guard digits are available, a careful choice of "a" is crucial. I would suggest a=9 for the 10-digit machines. Larger values do not seem to provide better accuracy. About 8 digits is all we get.

Addendum: This also means that a value as high as a=12,5 does not make much sense here. Try Namir's program with a=7 (replace 12.5 with 7 and 1.012 with 1.006). My tests returned results at least as accurate as with a=12,5 (mostly even better) while the program runs about twice as fast. BTW there is an error in the original listing: line 40 should read RCL 06.

Dieter