Hello Robert,
Thank you very much for your hint. I tried to implemented the code to detect the
overflow.
The following code use the macro ignore-errors which work for CLISP and the
tests float-nan-p and float-inf-p for GCL. Both tests are necessary. Numbers
which are bigger than 171.6 but near by give T with float-inf-p. Bigger numbers
e.g. gamma(250.0) give T with float-nan-p.
With this extension the last double float the algorithm works for is ~171.6243.
For bigger values a Maxima Error is generated.
I hope this code will work for other compilers too.
Do you think we should commit this extension?
Here is the modified routine gamma-lanczos:
(defun gamma-lanczos (z)
(declare (type (complex flonum) z)
(optimize (safety 3)))
(let ((c (make-array 15 :element-type 'flonum
:initial-contents
'(0.99999999999999709182
57.156235665862923517
-59.597960355475491248
14.136097974741747174
-0.49191381609762019978
.33994649984811888699e-4
.46523628927048575665e-4
-.98374475304879564677e-4
.15808870322491248884e-3
-.21026444172410488319e-3
.21743961811521264320e-3
-.16431810653676389022e-3
.84418223983852743293e-4
-.26190838401581408670e-4
.36899182659531622704e-5))))
(declare (type (simple-array flonum (15)) c))
(if (minusp (realpart z))
;; Use the reflection formula
;; -z*Gamma(z)*Gamma(-z) = pi/sin(pi*z)
;; or
;; Gamma(z) = pi/z/sin(pi*z)/Gamma(-z)
;;
;; If z is a negative integer, Gamma(z) is infinity. Should
;; we test for this? Throw an error?
;; The test will be done by the calling routine.
(/ (float pi)
(* (- z) (sin (* (float pi) z))
(gamma-lanczos (- z))))
(let* ((z (- z 1))
(zh (+ z 1/2))
(zgh (+ zh 607/128))
;; Calculate log(zp) to avoid overflow at this point.
(lnzp (* (/ zh 2) (log zgh)))
(ss
(do ((sum 0.0)
(pp (1- (length c)) (1- pp)))
((< pp 1)
sum)
(incf sum (/ (aref c pp) (+ z pp))))))
(let ((result
;; We check for an overflow. The last positive value in
;; double-float precicsion for which Maxima can calculate
;; gamma is ~171.6243 (CLISP 2.46 and GCL 2.6.8)
(ignore-errors
(* (sqrt (float (* 2 pi)))
(+ ss (aref c 0))
(exp (+ (- zgh) (* 2 lnzp)))))))
(cond ((null result)
;; No result. Overflow.
(merror "Overflow in `gamma-lanczos'."))
((or (float-nan-p (realpart result))
(float-inf-p (realpart result)))
;; Result, but beyond extreme values. Overflow.
(merror "Overflow in `gamma-lanczos'."))
(t result)))))))
Dieter Kaiser
-----Urspr?ngliche Nachricht-----
Von: robert.dodier at gmail.com [mailto:robert.dodier at gmail.com]
Gesendet: Samstag, 13. September 2008 16:43
An: Dieter Kaiser
Cc: maxima at math.utexas.edu
Betreff: Re: [Maxima] Problems with the Gamma function
On 9/8/08, Dieter Kaiser <drdieterkaiser at web.de> wrote:
> For big values we get further Lisp errors (this is reported as a bug
> SF[2013650]):
>
> (%i70) gamma(250.0);
> (%o70)
> Maxima encountered a Lisp error:
>
> Error in PROGN [or a callee]: Can't print a non-number.
This is a bug in GCL: it computes a non-numeric floating point
result (infinity or not-a-number) but it fails to print it.
My advice is to call FLOAT-NAN-P and FLOAT-INF-P
(see below) and take some appropriate action if either is true.
Incidentally, I've concluded that Maxima cannot expect the Lisp
implementation to provide floating point infinity and not-a-number;
the CL spec seems to exclude them. However, I think Maxima
should allow non-numeric floats if the underlying Lisp allows them.
I think that is the case at present.
Hope this helps,
Robert Dodier
PS. Here is the code, as it stands, to detect floating point infinity
and not-a-number. These functions are now in src/float.lisp.
I wouldn't be surprised if they can be improved.
(defun float-nan-p (x)
(and (floatp x) (not (= x x))))
(defun float-inf-p (x)
(and (floatp x) (not (float-nan-p x)) (beyond-extreme-values x)))
(defun beyond-extreme-values (x)
(multiple-value-bind (most-negative most-positive) (extreme-float-values x)
(cond
((< x 0) (< x most-negative))
((> x 0) (> x most-positive))
(t nil))))
(defun extreme-float-values (x)
;; BLECHH, I HATE ENUMERATING CASES. IS THERE A BETTER WAY ??
(case (type-of x)
(short-float (values most-negative-short-float most-positive-short-float))
(single-float (values most-negative-single-float
most-positive-single-float))
(double-float (values most-negative-double-float
most-positive-double-float))
(long-float (values most-negative-long-float most-positive-long-float))
;; NOT SURE THE FOLLOWING REALLY WORKS
;; #+(and cmu double-double)
;; (kernel:double-double-float
;; (values most-negative-double-double-float
most-positive-double-double-float))
))