--- maxima-pre59/src/csimp2.lisp Mon May 8 02:09:41 2000
+++ maxima-cvs-1025/src/csimp2.lisp Wed Nov 28 12:01:50 2001
@@ -241,8 +241,10 @@
((lessp q 1) ans)
(setq q (sub1 q) m (*dif m n) ans (times m ans))))
+#+nil
(declare-top (flonum a r))
+#+nil
(defun gammafloat (a)
(cond ((= a 1.0) 1.0)
((= a 0.0) (merror "GAMMA(0.0) has been generated."))
@@ -261,6 +263,7 @@
; (dbg:floating-exponent-overflow
; (merror "GAMMA(~A) - arithmetic overflow" a))))))
+#+nil
(defun do-gammafloat (a)
(do ((r 1.0 (*$ z r))
(s (minusp a)) (z (abs a)))
@@ -270,5 +273,75 @@
(t r)))
(setq z (1-$ z))))
+;; This implementation is based on Lanczos convergent formula for the
+;; gamma function for Re(z) > 0. We can use the reflection formula
+;;
+;; -z*Gamma(z)*Gamma(-z) = pi/sin(pi*z)
+;;
+;; to handle the case of Re(z) <= 0.
+;;
+;; See http://winnie.fit.edu/~gabdo/gamma.m for some matlab code to
+;; compute this and http://winnie.fit.edu/~gabdo/gamma.txt for a nice
+;; discussion of Lanczos method and an improvement of Lanczos method.
+;;
+;;
+;; The document says this should give about 15 digits of accuracy for
+;; double-precision IEEE floats. The document also indicates how to
+;; compute a new set of coefficients if you need more range or
+;; accuracy.
+
+(defun gamma-lanczos (z)
+ (declare (type (complex double-float) z)
+ (optimize (safety 3)))
+ (let ((g 607/128)
+ (c (make-array 15 :element-type 'double-float
+ :initial-contents
+ '(0.99999999999999709182d0
+ 57.156235665862923517d0
+ -59.597960355475491248d0
+ 14.136097974741747174d0
+ -0.49191381609762019978d0
+ .33994649984811888699d-4
+ .46523628927048575665d-4
+ -.98374475304879564677d-4
+ .15808870322491248884d-3
+ -.21026444172410488319d-3
+ .21743961811521264320d-3
+ -.16431810653676389022d-3
+ .84418223983852743293d-4
+ -.26190838401581408670d-4
+ .36899182659531622704d-5))))
+ (declare (type (rational 607/128 607/128) g)
+ (type (simple-array double-float (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? What
+ (/ (float pi 1d0)
+ (* (- z) (sin (* (float pi 1d0)
+ z))
+ (gamma-lanczos (- z))))
+ (let* ((z (- z 1))
+ (zh (+ z 1/2))
+ (zgh (+ zh 607/128))
+ (zp (expt zgh (/ zh 2)))
+ (ss
+ (do ((sum 0d0)
+ (pp (1- (length c)) (1- pp)))
+ ((< pp 1)
+ sum)
+ (incf sum (/ (aref c pp) (+ z pp))))
+ ))
+ (* (sqrt (float (* 2 pi) 1d0))
+ (+ ss (aref c 0))
+ (* zp (exp (- zgh)) zp))))))
+
+(defun gammafloat (a)
+ (realpart (gamma-lanczos (complex a 0d0))))
+
(declare-top (notype a r))
(declare-top (SPLITFILE ERF) (SPECIAL $NUMER $TRIGSIGN))
@@ -285,6 +358,7 @@
((and $trigsign (mminusp* y)) (neg (list '(%erf simp) (neg y))))
(t (eqtest (list '(%erf) y) x))))
+#+nil
(defmfun erf (y)
(cond ((> (abs y) 4.0) (cond ((> y 0.0) 1.0) (t -1.0)))
(t ((lambda (t1 xf)
@@ -308,6 +382,216 @@
t1))))
(cond ((> y 0.0) t1) (t (minus t1))))
0.0
- (abs y)))))
+ (abs y)))))
+
+;;; This is a hand-modified version of the code generated by f2cl
+;;; applied to the routine calerf from TOMS Algorithm 715. The changes are:
+;;;
+;;; o added some comments
+;;; o reindented some parts of the code
+;;; o changed the type integer4 to (signed-byte 32)
+;;; o changed AINT to floor
+;;; o removed the usage of the fref macro
+;;; o removed the fdo macro.
+;;; o Compute the constants instead of having the approximations given
+;;; in the Fortran code
+;;; o removed the arg result that was used to return the result. The
+;;; function value is the result.
+
+(let ((four 4.0d0)
+ (one 1.0d0)
+ (half 0.5d0)
+ (two 2.0d0)
+ (zero 0.0d0)
+ (sqrpi (coerce (/ (sqrt pi)) 'double-float))
+ (thresh 0.46875d0)
+ (sixten 16.0d0)
+ (xinf most-positive-double-float)
+ ;; XNEG is the negative of the solution of 2*exp(x*x) = XINF.
+ ;; Thus XNEG = -sqrt(log(XINF/2))
+ (xneg (- (sqrt (log (/ most-positive-double-float 2)))))
+ ;; argument below which erf(x) may be represented by
+ ;; 2*x/sqrt(pi) and above which x*x will not underflow.
+ ;; Conservatively, X such that 1+x=1.
+ (xsmall double-float-epsilon)
+ ;; largest argument acceptable to erfc; solution to the
+ ;; equation: W(x) * (1-0.5/x**2) = XMIN, where W(x) =
+ ;; exp(-x*x)/[x*sqrt(pi)].
+ ;;
+ ;; There's no analytic solution, and I'm too lazy to compute
+ ;; this more accurately and erfc would underflow in this case.
+ (xbig 26.543d0)
+ ;; Number for which 1-1/(2*x*x) = 1. That is, 1/(2*x*x) is
+ ;; double-float-negative-epsilon.
+ (xhuge (/ (sqrt (* 2 double-float-negative-epsilon))))
+ ;; Largest acceptable arg to erfcx; the minimum of XINF and
+ ;; 1/(sqrt(pi)*XMIN), where XMIN is the smallest positive
+ ;; floating-point number (normalized)
+ (xmax (min most-positive-double-float
+ (/ (* (coerce (sqrt pi) 'double-float)
+ #-gcl least-positive-normalized-double-float
+ #+gcl least-positive-double-float))))
+ (a
+ (make-array 5
+ :element-type
+ 'double-float
+ :initial-contents
+ '(3.1611237438705655d0 113.86415415105016d0
+ 377.485237685302d0 3209.3775891384694d0
+ 0.18577770618460318d0)))
+ (b
+ (make-array 4
+ :element-type
+ 'double-float
+ :initial-contents
+ '(23.601290952344122d0 244.02463793444417d0
+ 1282.6165260773723d0 2844.2368334391704d0)))
+ (c
+ (make-array 9
+ :element-type
+ 'double-float
+ :initial-contents
+ '(0.5641884969886701d0 8.883149794388377d0
+ 66.11919063714163d0 298.6351381974001d0
+ 881.9522212417692d0 1712.0476126340707d0
+ 2051.078377826071d0 1230.3393547979972d0
+ 2.1531153547440382d-8)))
+ (d
+ (make-array 8
+ :element-type
+ 'double-float
+ :initial-contents
+ '(15.744926110709834d0 117.6939508913125d0
+ 537.1811018620099d0 1621.3895745666903d0
+ 3290.7992357334597d0 4362.619090143247d0
+ 3439.3676741437216d0 1230.3393548037493d0)))
+ (p
+ (make-array 6
+ :element-type
+ 'double-float
+ :initial-contents
+ '(0.30532663496123236d0 0.36034489994980445d0
+ 0.12578172611122926d0 0.016083785148742275d0
+ 6.587491615298379d-4 0.016315387137302097d0)))
+ (q
+ (make-array 5
+ :element-type
+ 'double-float
+ :initial-contents
+ '(2.568520192289822d0 1.8729528499234604d0
+ 0.5279051029514285d0 0.06051834131244132d0
+ 0.0023352049762686918d0))))
+ (declare (type (simple-array double-float (6)) p)
+ (type (simple-array double-float (8)) d)
+ (type (simple-array double-float (9)) c)
+ (type (simple-array double-float (4)) b)
+ (type (simple-array double-float (5)) q a)
+ (type double-float xmax xhuge xbig xsmall xneg xinf sixten thresh
+ sqrpi zero two half one four))
+ (defun calerf (arg jint)
+ (declare (type (integer 0 2) jint)
+ (type double-float arg)
+ (optimize (speed 3)))
+ (prog ((del 0.0d0) (x 0.0d0) (xden 0.0d0) (xnum 0.0d0) (y 0.0d0)
+ (ysq 0.0d0) (result 0d0))
+ (declare (type double-float ysq y xnum xden x del result))
+ (setf x arg)
+ (setf y (abs x))
+ (cond
+ ((<= y thresh)
+ ;; Compute erf(x) for |x| < 0.46875
+ (setf ysq zero)
+ (if (> y xsmall)
+ (setf ysq (* y y)))
+ (setf xnum (* (aref a (- 5 1)) ysq))
+ (setf xden ysq)
+ #+nil
+ (loop for i of-type (integer 1 4) from 1 upto 3 do
+ (tagbody
+ (setf xnum (* (+ xnum (aref a (- i 1))) ysq))
+ (setf xden (* (+ xden (aref b (- i 1))) ysq))
+ label20))
+ (do ((i 1 (1+ i)))
+ ((> i 3))
+ (setf xnum (* (+ xnum (aref a (- i 1))) ysq))
+ (setf xden (* (+ xden (aref b (- i 1))) ysq)))
+
+ (setf result
+ (/ (* x (+ xnum (aref a (- 4 1))))
+ (+ xden (aref b (- 4 1)))))
+ (if (/= jint 0) (setf result (- one result)))
+ (if (= jint 2) (setf result (* (exp ysq) result))) (go label800))
+ ((<= y four)
+ ;; Compute erfc for 0.46785 <= |x| <= 4
+ (setf xnum (* (aref c (- 9 1)) y)) (setf xden y)
+ #+nil
+ (loop for i of-type (integer 1 8) from 1 upto 7 do
+ (tagbody
+ (setf xnum (* (+ xnum (aref c (- i 1))) y))
+ (setf xden (* (+ xden (aref d (- i 1))) y))
+ label120))
+ (do ((i 1 (1+ i)))
+ ((> i 7))
+ (setf xnum (* (+ xnum (aref c (- i 1))) y))
+ (setf xden (* (+ xden (aref d (- i 1))) y)))
+
+ (setf result
+ (/ (+ xnum (aref c (- 8 1)))
+ (+ xden (aref d (- 8 1)))))
+ (cond
+ ((/= jint 2) (setf ysq (/ (the (signed-byte 32) (floor (* y sixten))) sixten))
+ (setf del (* (- y ysq) (+ y ysq)))
+ (setf result (* (exp (* (- ysq) ysq)) (exp (- del)) result)))))
+ (t
+ ;; Compute erfc for |x| > 4
+ (setf result zero)
+ (cond
+ ((>= y xbig) (if (or (/= jint 2) (>= y xmax)) (go label300))
+ (cond ((>= y xhuge) (setf result (/ sqrpi y)) (go label300)))))
+ (setf ysq (/ one (* y y))) (setf xnum (* (aref p (- 6 1)) ysq))
+ (setf xden ysq)
+ #+nil
+ (loop for i of-type (integer 1 5) from 1 upto 4 do
+ (tagbody
+ (setf xnum (* (+ xnum (aref p (- i 1))) ysq))
+ (setf xden (* (+ xden (aref q (- i 1))) ysq))
+ label240))
+ (do ((i 1 (1+ i)))
+ ((> i 4))
+ (setf xnum (* (+ xnum (aref p (- i 1))) ysq))
+ (setf xden (* (+ xden (aref q (- i 1))) ysq)))
+
+ (setf result
+ (/ (* ysq (+ xnum (aref p (- 5 1))))
+ (+ xden (aref q (- 5 1)))))
+ (setf result (/ (- sqrpi result) y))
+ (cond
+ ((/= jint 2) (setf ysq (/ (the (signed-byte 32) (floor (* y sixten))) sixten))
+ (setf del (* (- y ysq) (+ y ysq)))
+ (setf result (* (exp (* (- ysq) ysq)) (exp (- del)) result))))))
+ label300
+ (cond
+ ((= jint 0) (setf result (+ (- half result) half))
+ (if (< x zero) (setf result (- result))))
+ ((= jint 1) (if (< x zero) (setf result (- two result))))
+ (t
+ (cond
+ ((< x zero)
+ (cond ((< x xneg) (setf result xinf))
+ (t
+ (setf ysq (/ (the (signed-byte 32) (floor (* x sixten))) sixten))
+ (setf del (* (- x ysq) (+ x ysq)))
+ (setf y (* (exp (* ysq ysq)) (exp del)))
+ (setf result (- (+ y y) result))))))))
+ label800
+ (go end_label)
+ end_label
+ (return result))))
+
+(defmfun erf (y)
+ (calerf (float y 1d0) 0))
+(defmfun erfc (y)
+ (calerf (float y 1d0) 1))
+
(declare-top (SPLITFILE EMATRIX))