Subject: Increasing the accuracy of Gamma for double float
From: Dieter Kaiser
Date: Tue, 21 Oct 2008 18:14:52 +0200
Von: raymond.toy at ericsson.com [mailto:raymond.toy at ericsson.com]
>Dieter Kaiser wrote:
>> I have done some further work to improve the precision of the Gamma function
for
>> double float calculation.
>>
>> This is the actual result for gamma-lanczos:
>> ...................................................
>> Test of precision for gamma-lanczos
>> gamma((float(z)) against gamma(bfloat(z)), fpprec:32
>> z: 0.5 thru 150.0 step 0.25 (values about 170.0 give an overflow)
>I modified gammafloat in csimp2.lisp to:
>(defun gammafloat (a)
> (slatec::dgamma (float a)))
>A quick test says
>meanerr = 2.35b-14
>maxerr = 1.27b-13 at 142.75
>So gamma-lanczos is actually slightly more accurate than the SLATEC
>routine. I'm a bit surprised. I would have expected the SLATEC routine
>to do better.
>Also, the slatec routine does pretty well until x = 10 at which point
>the errors start to grow. x = 10 happens to be the point where the
>slatec routine switches from a series to something like Stirling's
>formula. The code claims the approximation there is accurate to
>1.28d-31 (which we can't actually achieve since we only have
>double-float precision). Perhaps the computation of (x-1/2*log(x)-x
>causes the loss of precision we are seeing.
I am not really surprised. It is no problem to get any desired precision. But
calculating gamma we always need at least 3 to 5 extra digits. In double float
precision we have no extra digits but lost the 3 to 5 digits we need to
calculate gamma.
The reason might be that the algorithm for gamma are based on expansions with
alternating coefficients.
Dieter Kaiser