Gamma Function Using Spouge's Method

[Dieter]

(08-25-2015 11:29 PM)lcwright1964 Wrote:  This is excellent. I have learned that Free42 is an excellent place to observe the behaviour of HP41 programs. The programs are completely compatible without modification...

I agree that Free42 is an excellent piece of software, but I would not recommend it for testing HP41 programs. The different precision and working range may and will lead to different results. For instance numbers may not round to zero while they do on a physical '41, so x=0? and x=y? tests may behave differently and iterations may not terminate. Since the '42 supports complex arithmetics, some operations (e.g. the log of a negative number) will return a complex result instead of an error message (which may be the intended behaviour of the original HP41 program). The same is true for the larger working range for which the program may not have been designed. In general, Free42 does not use HP code, so possibly even simple 12-digit results may differ from those of a hardware '42.

But of course there is a solution: V41 essentially behaves like a real HP41 and so I think this is the best option for testing programs for this calculator.

(08-25-2015 11:29 PM)lcwright1964 Wrote:  ...12-digits are displayed, and everything is done internally to 25 digits (not 34 like you say--that is DP mode in the WP34s).

This may be true for the original Free42, but there also are newer versions that use 34 digit BCD arithmetics ("Free42 Decimal"). And so indeed 3 [sqrt] [x²] 3 [–] returns –1E–33. Take a look at the Free42 website, where it says:

"Free42 Decimal uses the Intel Decimal Floating-Point Math Library; it uses IEEE 754-2008 quadruple precision decimal floating-point, which consumes 16 bytes per number, and gives 34 decimal digits of precision, with exponents ranging from -6143 to +6144."

(08-25-2015 11:29 PM)lcwright1964 Wrote:  Indeed, working with JMB's GAM+ routine I have concluded there is likely little benefit to adding a term to the Stirling series. With ample available precision on Free42 GAM+ unedited computes the Gamma of 5.000000001, just a tad about the shift-divide cutoff, to about 10.4 edd. A further term can't improve on this on the HP41--indeed, even more computations will add to rounding error. And there is little to be gained by lowering the shift-divide cutoff, as adding a term likely adds more operations than what might be saved by lowering the cutoff.

This essentially matches my results. I think the only way to get 10 digits out of an HP41 is 13-digit MCode. Here my approximations may be useful since they deliver the same accuracy level as the current Sandmath GAMMA implementation, but with two terms less. I tried the n=4 approximation on my trusted 12-digit HP35s and I got 10 valid digits, here and there the last one rounded the wrong way (i.e. ±1 ULP). Maybe Angel is listening. ;-)

BTW I now have an optimized approximation for n=5. The largest relative error for x = 0...70 (i.e. 1...71 for Gamma) now is ±3,9 E–13. I did this on Free42 with 34-digit support, so I'm not sure what the required working precision might be. I guess something like 16 digits would get you more or less close to this.

(08-25-2015 11:29 PM)lcwright1964 Wrote:  For your Lanczos formula your reciprocals start at 1/z+1 instead of 1/z in my versions. Does this permit you to make sure the error curve passes through z = 0?

That's simply because my approximations calculate x! resp. Gamma(x+1). Yes, the error is zero at z=0, but that's simply because that's one of the nodes in the Remez simulation. ;-)

Dieter