Gamma Function by Stieltjes Continued Fraction

[Dieter]

(09-04-2015 01:54 AM)lcwright1964 Wrote:  Indeed, my approach was taken from JMB's GAM+ routine, where he does this. I wonder if your modification would improve that too?

The proposed modification divides e^(quite_a_small_value) by e^x, as opposed to e^(quite_a_small_value – x). Since x >> quite_a_small_value, 2 or 3 digits are lost after the subtraction. On the other hand the two individual antilogs are more or less exact (within working precision). That's why the results are more accurate. Since the mentioned small value is close to zero I even tried a method that uses e^x–1, but that does not make much of a difference.

(09-04-2015 01:54 AM)lcwright1964 Wrote:  With adequate extra precision EDD = 10.4 is plenty to get 10 full digits across the board,

Let me add a short remark on these EDD figures. EDD=n means a relative error of 10^–n, which in turn means the error can approach 1 unit in the n-th significant digit. If that is the last displayed digit, the error of the displayed result may be up to 1,5 ULP. So EDD=n means that the last digit can be off by 1 or slightly more. 10,4 EDD equals a relative error of 4 E–10 or up to 0,9 ULP including display roundoff. That's not exactly what I'd call "plenty". ;-) I prefer a somewhat stricter measure: for n valid digits the relative error should not exceed 10^–(n+1), i.e. 1 E–11 for 10 digits. This makes sure the n-digit rounded result is within 0,6 ULP, which matches HP's accuracy for most of the transcendental functions. But since we do all the calculations with merely the same 10 digits, we should not be too picky here... ;-)

(09-04-2015 01:54 AM)lcwright1964 Wrote:  I am figuring that tweaks can start to get out of hand after a bit without driving one crazy,

Exactly. ;-)

Dieter