Subject: Increasing the accuracy of Gamma for double float
From: Dieter Kaiser
Date: Wed, 22 Oct 2008 23:23:52 +0200
-----Urspr?ngliche Nachricht-----
Von: raymond.toy at ericsson.com [mailto:raymond.toy at ericsson.com]
>Ah, we don't have to. The following works pretty well:
>(defun gammafloat (a)
> (let ((a (float a)))
> (cond ((< (abs a) 10)
> (slatec::dgamma a))
> (t
> ;; Stirling approximation
> (* (expt a (- a 1/2))
> (exp (- a))
> (sqrt (* 2 (float pi)))
> (exp (slatec::d9lgmc a)))))))
>For z from 10 to 140 step 1/4 the peak relative error is about 4d-16.
>The code above needs some work because it overflows around 142, and
>doesn't handle negative values of a, but we know how to handle that.
>My conclusion from this is that computing gamma from log gamma
>needlessly loses digits because log gamma itself wastes digits. A
>better way is to use Stirling's approximation with an accurate value for
>the "correction" term.
>This idea could also be applied to the computation of the beta function.
> It's more coding work to get this correct, but with some care, we won't
>lose accuracy, and we can probably not prematurely overflow when the
>beta function itself would fit in a double-float.
Thus it seems to be more useful to use different algorithm for different regions
of the Gamma function.
I had a look at functions.wolfram.com. There I found under the topic "Asymptotic
series expansions" for the Gamma function first the Stirling formula (Ref.
06.05.06.0010.01). What is interessting is the generalization of the Stirling
formula to complex numbers (formula 06.05.06.0011.01) and (formula
06.05.06.0041.01). These formulas I think represent the algorithm of bffac.
Perhaps the algorithm of d9lgmc is nothing else than the sum of the first terms
of the series of Bernoulli numbers. This could be a hint to extend the
alogorithm to complex numbers too.
Dieter Kaiser