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:
> 

Here is a suggested replacement for gammafloat:

(defun gammafloat (a)
  (let ((a (float a)))
    (cond ((minusp a)
	   (/ (float (- pi))
	      (* a (sin (* (float pi) a)))
	      (gammafloat (- a))))
	  ((< a 10)
	   (slatec::dgamma a))
	  (t
	   (let ((c (* (sqrt (* 2 (float pi)))
		       (exp (slatec::d9lgmc a)))))
	     (let ((v (expt a (- (* .5d0 a) 0.25d0))))
	       (* v
		  (/ v (exp a))
		  c)))))))

I think this handles all the cases we need.  A test for z from 1/2 to
150 by 1/4 compared against the bfloat version with 64 digits shows that
the peak relative error is 5.23e-16, with a mean error of 1.45e-16.

I think that's pretty good.

Barring any bugs (I haven't run the testsuite yet), I would like to
check this in.

Ray