Gamma Function Using Spouge's Method

[Dieter]

(08-26-2015 12:50 PM)lcwright1964 Wrote:  But of course, and that's kind of my point. What I should have been clear about is that it lets me test various Gamma implementations, written in FOCAL, in a high precision environment where 12 digits--two more--are show. So I can play around with extra terms and shift-divide cutoffs and see where it is I get 10 digit accuracy where this is ample working precision.

OK – that's of course a useful application. I often do this kind of stuff with the 34s. Here you can even choose between 16 and 34 digits. ;-)

(08-26-2015 12:50 PM)lcwright1964 Wrote:  For example, in the GAM+ program of JMB that I love so much, it is evident that he worked it up in a higher precision environment, because the number of terms used and the selection of 5 as the shift/divide threshold gives at LEAST 10 good digits when there is ample working precision. This will be degraded somewhat in the HP41 10 digit environment, but adding another term, as contemplated, will NOT improve it there. It will give more digits when they are available, such as on the Free42, but won't confer any benefit on the HP41. That's what I mean.

That's essentially the same point I mentioned in a post earlier in this thread regarding the Spouge method, when I said that a=7 instead of 12.5 will yield the same or even better accuracy on a real '41.

(08-26-2015 12:50 PM)lcwright1964 Wrote:  And thank you for correcting me on the number of internal digits in Free42 Decimal. It used to be about 25-digit precision when Thomas used the BCD-10000 library. I didn't know he switched to the Intel quadruple precision library. I do suspect that not all of the Free42's functions have been extended to the higher precision.
(...)
Testing this against the Free42's built in GAM shows the latter to do no better than about 23 or 24 digits for a few arguments, so the older code is still there despite the more precise quadruple precision numbers.

Ah, yes, this may be possible. This is were the WP34s performs better. Usually you get 30+ digits in DP mode. Although only the 16 digits in SP mode are granted.

Finally, here are the coefficients for the n=5 case, optimized for x=0..70 (resp. 1...71 for Gamma). They differ slightly from what you posted, but the largest relative error here is as small as 3,75 E–13.

Code:

 c =  5,081
a0 =  9,44696770446255924982608 E-3
a1 =  1,36718522540552597224102
a2 = -1,83953548282705973597554
a3 =  6,83853459467545887340920 E-1
a4 = -6,54290299376343116198247 E-2
a5 =  6,57205667880559284049346 E-4

The values are rounded to 24 digits, but I think 16 should do. I did a short test on a WP34s in SP mode and it looks like the error stays below 3,8 E–13.

Edit: I did a final iteration which, for c = 5,081 and sufficient working precision, gave a max. rel. error of 3,7491 E–13 over the mentioned interval.

Dieter