[Maxima] Increasing the accuracy of Gamma for double float

[Dieter Kaiser]

Hello Ray,

I would appreciate if we get the more accurate float precicison for gamma. I
think a lot of more special functions would get more accurate also.

We had implemented an overflow check for gamma-lanzcos which is now again
missing. Therefore one test (Problem 196) which check the overflow code failes
with GCL (not with CLISP).

Do you would like to implement the overflow check again?

Dieter Kaiser

-----Ursprüngliche Nachricht-----
Von: raymond.toy at ericsson.com [mailto:raymond.toy at ericsson.com] 
Gesendet: Donnerstag, 30. Oktober 2008 15:46
An: Dieter Kaiser
Cc: maxima at math.utexas.edu
Betreff: Re: AW: [Maxima] 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