Increasing the accuracy of Gamma for double float

Dieter Kaiser wrote:
> 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.

How did you arrive at this conclusion?  Why does the computation require
3 to 5 extra digits?

> 
> The reason might be that the algorithm for gamma are based on expansions with
> alternating coefficients.

I don't think that's the problem.  The coefficients in the series have
31 digits.  The weighted error is 1.28e-31.  Don't know what the
weighting function is, but that seems reasonable for something that only
has 31 digits.

Ray