Gamma Function by Stieltjes Continued Fraction

[Dieter]

(09-13-2015 02:06 AM)lcwright1964 Wrote:  Hey, Dieter, for HP41 work with 10 digits I think I am sticking with the refinement mentioned above. Expressing the coefficient that is around 5 to a few more digits doesn't seem to confer extra benefit on the HP41, and indeed when I plot the error curve in an arbitrary precision environment I really see no big improvement at all.

The additional digits only help squeeze the error a tiny bit below the 2 E–12 mark. ;-) You won't notice much of a difference even in a higher precision environment. For the '41 and other 10-digit machines the earlier methods are perfectly OK.

BTW, with a=4,01 and b=5,00008 you can lower the threshold to 8 and still get similar accuracy (error < 5 E–12), while for small x one loop less is required.

(09-13-2015 02:06 AM)lcwright1964 Wrote:  I should play around with your four term versions, but I suspect that the advantages of them are lost in a 10 digit environment and, if anything, the extra computations due to the extra term produce their own round off issues.

Sure, all this was done with higher-precision environments in mind. The best I got so far is an error < 3,7 E–14 (threshold = 8), and I tried this on a WP34s in SP mode (16 digits). It's nice to see the 12 digit display show the same result as the built-in Gamma function. ;-)

(09-13-2015 02:06 AM)lcwright1964 Wrote:  Of course, all of these great refinements are good candidates for M-CODE programming, but I think none of them can top some of the optimized Lanczos approximations we were discussing elsewhere. Getting rid of the shift-divide step for smaller arguments really goes far to speeding things up across the board.

Right, this shift-and-divide thing slows down everything, especially on older 10-digit calculators like the '41 or 67/97. But in MCode speed should be no big issue, and with some careful, accuracy-preserving coding one may expect something like 11 valid digits.

Finally I would like to sum up some results for three-term Stieltjes approximations in the attached graphics. The y-axis is the EDD figure, i.e. the base-10 log of the relative error, or roughly the number of exact digits (±1 unit). The peaks are the points where the error passes through zero and accuracy becomes infinite. The sawtooth pattern for small x is due to the shift-and-divide method applied there. You can see the original Stieltjes method in the black graph, while the red one shows the improvement due to the a=4 simplification. I especially like the blue version as it reaches an error level slightly below 1 E–11 down to t=8. Yes, you may omit the ...04 at the end of the b coefficient. ;-)

Dieter