Increasing the accuracy of Gamma for double float

Dieter Kaiser wrote:
> Ray wrote:
>> I am looking into this again, finally.  I was wondering why you made
>> the following change to gamma-lanczos:
>>
>> 	(let* ((z (- z 1))
>> 	       (zh (+ z 1/2))
>> 	       (zgh (+ zh 607/128))
>>               ;; Calculate log(zp) to avoid overflow at this point.
>> 	       (lnzp (* (/ zh 2) (log zgh)))
>>
>> You say it overflows, but I don't see it.  It used to be zp =
>> zgh^(zh/2).
>>
>> For z up to at least 200, this would not overflow.  Yes, the final
>> result will overflow when we compute zp*exp(-zgh)*zp, but that's to be
>> expected.  zp was done this way so we wouldn't prematurely overflow.
> 
> Yes, zp overflows not as fast as the result. I think I have tried to avoid an
> overflow of an intermediate result. Another possible was to to catch an error
> for the whole calculation. When we get an overflow of an intermediate result I
> expect again a Lisp error with CLISP. GCL might continue the calculation and the
> test and the end of the function might be enough.

Ah, ok. I understand the problem now.  I think we can move the
calculation of zp into the ignore-errors part, and that should take care
of it.  Using log loses accuracy, so we should avoid that if possible.

Ray