Increasing the accuracy of Gamma for double float

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

Yes, I had guessed bffac used asymptotic series for log gamma, but I was
never really sure.

> 
> Perhaps the algorithm of d9lgmc is nothing else than the sum of the first terms
> of the series of Bernoulli numbers. 

Difficult to say.  The comment in d9lgmc says:

C Compute the log gamma correction factor for X .GE. 10. so that
C LOG (DGAMMA(X)) = LOG(SQRT(2*PI)) + (X-.5)*LOG(X) - X + D9lGMC(X)

It's not clear if d9lgmc is an approximation to the real correction term
or to the bernoulli series.  In any case the coefficients used there
don't match the bernoulli series.  The first few are close, but
definitely aren't the same.

I'm not sure if we can extend d9lgmc itself to the complex domain.
There shouldn't be any issue with using the asymptotic formula though.

Ray