## maxima-commits

 [Maxima-commits] CVS: maxima-pre59/src csimp2.lisp,1.1.1.1,1.2 From: Raymond Toy - 2002-01-09 23:04:49 ```Update of /cvsroot/maxima/maxima-pre59/src In directory usw-pr-cvs1:/tmp/cvs-serv25102 Modified Files: csimp2.lisp Log Message: o Added more accurate gamma routine that should be accurate to within 15 digits or so. Still needs a bit of work because gamma(10.0) returns 362880.0000000005 instead of 362880.0. o Added more accurate erf function. Index: csimp2.lisp =================================================================== RCS file: /cvsroot/maxima/maxima-pre59/src/csimp2.lisp,v retrieving revision 1.1.1.1 retrieving revision 1.2 diff -u -d -r1.1.1.1 -r1.2 --- csimp2.lisp 2000/05/08 06:09:41 1.1.1.1 +++ csimp2.lisp 2002/01/09 23:04:46 1.2 @@ -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,6 +273,76 @@ (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,7 +382,217 @@ 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)) ```