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