Increasing the accuracy of Gamma for double float

Raymond Toy wrote:
> Dieter Kaiser wrote:
>> -----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.
>>
> 
> Our conclusion about the accuracy of gamma is wrong.  According to this
> <http://www.netlib.org/cephes/doubldoc.html#gamma>;, the peak relative
> error of 2.3d-15 can be achieved:
> 
>  *                      Relative error:
>  * arithmetic   domain     # trials      peak         rms
>  *    DEC      -34, 34      10000       1.3e-16     2.5e-17
>  *    IEEE    -170,-33      20000       2.3e-15     3.3e-16
>  *    IEEE     -33,  33     20000       9.4e-16     2.2e-16
>  *    IEEE      33, 171.6   20000       2.3e-15     3.2e-16
> 
> We should examine the algorithm used there.

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.

Ray