Gamma Function by Stieltjes Continued Fraction

[lcwright1964]

Hey Dieter, thank you for all that.

There is no doubt now that across the range that your choice to take the antilog of both results and divide the them as opposed to subtract the value from the CF computation and antilog the whole thing is the better choice. Indeed, my approach was taken from JMB's GAM+ routine, where he does this. I wonder if your modification would improve that too?

I had figured out the line addressing issue as a step saver before I read your post, and I was keen to try it out when I got the chance. In my original and in the streamlined version using LastX there is a saving of 3 steps due to loss of the LBL lines but a gain of one step because .5 takes two steps. I also figured that on the 33 and 25 step 49 doesn't have to be an R/S or RTN as it just stops there automatically.

I would like to comment on my choice of 7 as the shift and divide threshold. On Mathematica it was easy to plot out exact decimal digits against x (where EDD = -Log10[Abs[Relative Error]]). The curve sweeps up to about EDD = 10.4 when x = 7. Then there is a spike a little after that, indicating that the approximation briefly becomes really accurate around x = 8 for some reason. The curve then swoops back down to to about EDD = 10.4 at roughly x =8.5, from whence it gradually climbs to EDD = 11 at x = 13 or so as you have observed. With adequate extra precision EDD = 10.4 is plenty to get 10 full digits across the board, and I figured that taking the cut-off higher would slow things down for smaller values of x without much added benefit. So that was my logic there.

I am figuring that tweaks can start to get out of hand after a bit without driving one crazy, but I think we have taken this pretty far. Thanks!

Les