Increasing the accuracy of Gamma for double float

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.

I have not done further tests to check it. Today I have no more time to do more
tests.

Dieter Kaiser