diff --git a/code/irrat-dd.lisp b/code/irrat-dd.lisp index 24532471f35c0377da34648fcd4e20e7b7580465..6ff711d3f7f9102b0024810f90e332358bbe543a 100644 --- a/code/irrat-dd.lisp +++ b/code/irrat-dd.lisp @@ -5,7 +5,7 @@ ;;; Carnegie Mellon University, and has been placed in the public domain. ;;; (ext:file-comment - "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/irrat-dd.lisp,v 1.13 2007/05/29 16:28:36 rtoy Exp $") + "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/irrat-dd.lisp,v 1.14 2007/06/11 19:01:23 rtoy Exp $") ;;; ;;; ********************************************************************** ;;; @@ -1703,15 +1703,6 @@ Z may be any number, but the result is always a complex." ;; ATANH is continuous with quadrant III in this case. (dd-complex-atanh (complex z -0d0))) (t - (flet ((careful-mul (a b) - ;; Carefully multiply a and b, taking care to handle - ;; signed zeroes. Only need to handle the case of b - ;; being zero. - (if (zerop b) - (if (minusp (* (float-sign a) (float-sign b))) - -0w0 - 0w0) - (* a b)))) (let* ( ;; Constants (theta (/ (sqrt most-positive-double-float) 4.0w0)) (rho (/ 4.0w0 (sqrt most-positive-double-float))) @@ -1719,7 +1710,7 @@ Z may be any number, but the result is always a complex." (rp (float (realpart z) 1.0w0)) (beta (float-sign rp 1.0w0)) (x (* beta rp)) - (y (careful-mul beta (- (float (imagpart z) 1.0w0)))) + (y (* beta (- (float (imagpart z) 1.0w0)))) (eta 0.0w0) (nu 0.0w0)) ;; Shouldn't need this declare. @@ -1758,7 +1749,7 @@ Z may be any number, but the result is always a complex." (+ (square (- 1.0d0 x)) (square t1)))))) (setf nu (* 0.5d0 - (dd-%atan2 (careful-mul 2.0d0 y) + (dd-%atan2 (* 2.0d0 y) (- (* (- 1.0d0 x) (+ 1.0d0 x)) (square t1)))))))) diff --git a/compiler/float-tran.lisp b/compiler/float-tran.lisp index b75357ecfdcd00ac9187fc2844594b35f1a910ae..4f6cde14a61d8a466f025e8328428954fd838b0c 100644 --- a/compiler/float-tran.lisp +++ b/compiler/float-tran.lisp @@ -5,7 +5,7 @@ ;;; Carnegie Mellon University, and has been placed in the public domain. ;;; (ext:file-comment - "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/compiler/float-tran.lisp,v 1.111 2007/05/09 03:18:09 rtoy Exp $") + "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/compiler/float-tran.lisp,v 1.112 2007/06/11 19:01:23 rtoy Exp $") ;;; ;;; ********************************************************************** ;;; @@ -1630,6 +1630,8 @@ (multiple-value-bind (s1 s2) (two-sum a0 b0) (declare (double-float s1 s2)) + (when (float-infinity-p s1) + (return-from add-dd (values s1 0d0))) (multiple-value-bind (t1 t2) (two-sum a1 b1) (declare (double-float t1 t2)) @@ -1638,7 +1640,13 @@ (quick-two-sum s1 s2) (declare (double-float s1 s2)) (incf s2 t2) - (quick-two-sum s1 s2))))) + (multiple-value-bind (r1 r2) + (quick-two-sum s1 s2) + (if (and (zerop a0) (zerop b0)) + ;; Handle sum of signed zeroes here. + (values (float-sign (+ a0 b0) 0d0) + 0d0) + (values r1 r2))))))) (deftransform + ((a b) (vm::double-double-float vm::double-double-float) *) `(multiple-value-bind (hi lo) @@ -1671,6 +1679,8 @@ (multiple-value-bind (s1 s2) (two-diff a0 b0) (declare (double-float s2)) + (when (float-infinity-p s1) + (return-from add-dd (values s1 0d0))) (multiple-value-bind (t1 t2) (two-diff a1 b1) (incf s2 t1) @@ -1678,7 +1688,12 @@ (quick-two-sum s1 s2) (declare (double-float s2)) (incf s2 t2) - (quick-two-sum s1 s2))))) + (multiple-value-bind (r1 r2) + (quick-two-sum s1 s2) + (if (and (zerop a0) (zerop b0)) + (values (float-sign (- a0 b0) 0d0) + 0d0) + (values r1 r2))))))) (declaim (inline sub-d-dd)) (defun sub-d-dd (a b0 b1) @@ -1688,8 +1703,14 @@ (multiple-value-bind (s1 s2) (two-diff a b0) (declare (double-float s2)) + (when (float-infinity-p s1) + (return-from add-dd (values s1 0d0))) (decf s2 b1) - (quick-two-sum s1 s2))) + (multiple-value-bind (r1 r2) + (quick-two-sum s1 s2) + (if (and (zerop a) (zerop b0)) + (values (float-sign (- a b0) 0d0) 0d0) + (values r1 r2))))) (declaim (inline sub-dd-d)) (defun sub-dd-d (a0 a1 b) @@ -1699,8 +1720,14 @@ (multiple-value-bind (s1 s2) (two-diff a0 b) (declare (double-float s2)) + (when (float-infinity-p s1) + (return-from add-dd (values s1 0d0))) (incf s2 a1) - (quick-two-sum s1 s2))) + (multiple-value-bind (r1 r2) + (quick-two-sum s1 s2) + (if (and (zerop a0) (zerop b)) + (values (float-sign (- a0 b) 0d0) 0d0) + (values r1 r2))))) (deftransform - ((a b) (vm::double-double-float vm::double-double-float) *) `(multiple-value-bind (hi lo) @@ -1854,10 +1881,17 @@ (multiple-value-bind (p1 p2) (two-prod a0 b) (declare (double-float p2)) + (when (float-infinity-p p1) + (return-from mul-dd (values p1 0d0))) ;;(format t "mul-dd-d p1,p2 = ~A ~A~%" p1 p2) (incf p2 (* a1 b)) ;;(format t "mul-dd-d p2 = ~A~%" p2) - (quick-two-sum p1 p2))) + (multiple-value-bind (r1 r2) + (quick-two-sum p1 p2) + (when (zerop r1) + (setf r1 (float-sign p1 0d0)) + (setf r2 p1)) + (values r1 r2)))) (declaim (maybe-inline mul-dd)) (defun mul-dd (a0 a1 b0 b1) @@ -1867,11 +1901,15 @@ (multiple-value-bind (p1 p2) (two-prod a0 b0) (declare (double-float p1 p2)) + (when (float-infinity-p p1) + (return-from mul-dd (values p1 0d0))) (incf p2 (* a0 b1)) (incf p2 (* a1 b0)) - (multiple-value-bind (p1 p2) + (multiple-value-bind (r1 r2) (quick-two-sum p1 p2) - (values p1 p2)))) + (if (zerop r1) + (values (float-sign p1 0d0) 0d0) + (values r1 r2))))) (declaim (maybe-inline add-dd-d)) (defun add-dd-d (a0 a1 b) @@ -1881,8 +1919,14 @@ (multiple-value-bind (s1 s2) (two-sum a0 b) (declare (double-float s1 s2)) + (when (float-infinity-p s1) + (return-from add-dd (values s1 0d0))) (incf s2 a1) - (quick-two-sum s1 s2))) + (multiple-value-bind (r1 r2) + (quick-two-sum s1 s2) + (if (and (zerop a0) (zerop b)) + (values (float-sign (+ a0 b) 0d0) 0d0) + (values r1 r2))))) (declaim (maybe-inline sqr-dd)) (defun sqr-dd (a0 a1) @@ -1963,6 +2007,8 @@ (optimize (speed 3)) (inline sub-dd)) (let ((q1 (/ a0 b0))) + (when (float-infinity-p q1) + (return-from div-dd (values q1 0d0))) ;; (q1b0, q1b1) = q1*(b0,b1) ;;(format t "q1 = ~A~%" q1) (multiple-value-bind (q1b0 q1b1)