Subject: Increasing the accuracy of Gamma for double float
From: Dieter Kaiser
Date: Tue, 21 Oct 2008 22:27:41 +0200
-----Urspr?ngliche Nachricht-----
Von: raymond.toy at ericsson.com [mailto:raymond.toy at ericsson.com]
>FWIW, I played around with the dlngam function in slatec. For z from 10
>to 150, dlngam differs from log_gamma by about 2.75d-16 (relative
>error). Not too bad compared to double-float-epsilon of about 1.1d-16.
>If we compare gamma(z) to exp(dlngam(z)), the max rel error is about
>1.3d-13. Assuming gamma(z) with fpprec = 32 is correct, the error comes
>from computing the exp. Not much we can do about that, I think.
Because of this result concerning the accuracy of log_gamma, gamma and the fact
that we have to do an exponentiation I have concluded some time ago that it
could be better for the beta function to implement an algorithm using the
log_gamma function for double float precision.
Now we are using the gamma function directly. A second problem with this
algorithm is, that we get very fast an overflow even if the result would be
small enough to fit into a double float value.
This is the code we use to calculate the beta function in double float
precision:
((or (and (floatp u) (floatp v))
(and $numer (numberp u) (numberp v)))
($makegamma (list '($beta) u v)))
Dieter Kaiser