Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv12645
Modified Files: qd-dd.lisp qd-package.lisp qd.lisp Log Message: To speed up Allegro (and other Lisp's that don't support inline functions), change QUICK-TWO-SUM from a function to a macro. Note that macro has different calling convention than the function. This is needed because Allegro apparently doesn't handle VALUES without boxing.
All rt tests pass.
qd-package.lisp: o For CMUCL, don't import C::QUICK-TWO-SUM into the QDI package anymore.
qd-dd.lisp: o New QUICK-TWO-SUM macro.
qd.lisp: o Add CMUCL version of QUICK-TWO-SUM macro, which just calls C::QUICK-TWO-SUM. o Update all users of QUICK-TWO-SUM appropriately.
--- /project/oct/cvsroot/oct/qd-dd.lisp 2007/09/13 17:28:30 1.8 +++ /project/oct/cvsroot/oct/qd-dd.lisp 2007/09/16 02:46:24 1.9 @@ -31,6 +31,7 @@ ;;; ;;; These routines were taken directly from CMUCL.
+#|| (declaim (inline quick-two-sum)) (defun quick-two-sum (a b) "Computes fl(a+b) and err(a+b), assuming |a| >= |b|" @@ -40,6 +41,16 @@ (e (- b (- s a)))) (declare (double-float s e)) (values s e))) +||# + +(defmacro quick-two-sum (s e x y) + (let ((a (gensym)) + (b (gensym))) + `(let* ((,a ,x) + (,b ,y)) + (declare (double-float ,a ,b ,s ,e)) + (setf ,s (+ ,a ,b)) + (setf ,e (- ,b (- ,s ,a))))))
(declaim (inline two-sum)) (defun two-sum (a b) @@ -53,6 +64,21 @@ (declare (double-float s v e)) (values s e)))
+#+nil +(defmacro two-sum (s e x y) + "Computes fl(a+b) and err(a+b)" + (let ((a (gensym)) + (b (gensym)) + (v (gensym)) + `(let ((,a ,x) + (,b ,y)) + (declare (double-float ,a ,b)) + (setf ,s (+ ,a ,b)) + (let ((,v (- ,s ,a))) + (declare (double-float v)) + (setf e (+ (- ,a (- ,s ,v)) + (- ,b ,v)))))))) + (declaim (inline two-prod)) (declaim (inline split)) ;; This algorithm is the version given by Yozo Hida. It has problems --- /project/oct/cvsroot/oct/qd-package.lisp 2007/09/12 21:01:13 1.37 +++ /project/oct/cvsroot/oct/qd-package.lisp 2007/09/16 02:46:24 1.38 @@ -98,7 +98,6 @@ #+cmu (:import-from #:c #:two-sum - #:quick-two-sum #:two-prod #:two-sqr))
--- /project/oct/cvsroot/oct/qd.lisp 2007/09/12 02:31:14 1.46 +++ /project/oct/cvsroot/oct/qd.lisp 2007/09/16 02:46:24 1.47 @@ -132,6 +132,16 @@ div-qd-d div-qd-dd))
+#+cmu +(defmacro quick-two-sum (s e x y) + `(multiple-value-setq (,s ,e) + (c::quick-two-sum ,x ,y))) + +#+(and nil cmu) +(defmacro two-sum (s e x y) + `(multiple-value-setq (s e) + (c::two-sum x y))) + #-(or qd-inline (not cmu)) (declaim (ext:start-block renorm-4 renorm-5 make-qd-d @@ -172,91 +182,67 @@
(defun renorm-4 (c0 c1 c2 c3) (declare (double-float c0 c1 c2 c3) - (optimize (speed 3) (safety 0))) - (let ((s2 0d0) + (optimize (speed 3) (safety 0) (debug 0))) + (let ((s0 0d0) + (s1 0d0) + (s2 0d0) (s3 0d0)) - (multiple-value-bind (s0 c3) - (quick-two-sum c2 c3) - (multiple-value-bind (s0 c2) - (quick-two-sum c1 s0) - (multiple-value-bind (c0 c1) - (quick-two-sum c0 s0) - (let ((s0 c0) - (s1 c1)) - (cond ((/= s1 0) - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c2)) - (if (/= s2 0) - (multiple-value-setq (s2 s3) - (quick-two-sum s2 c3)) - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c3)))) - (t - (multiple-value-setq (s0 s1) - (quick-two-sum s0 c2)) - (if (/= s1 0) - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c3)) - (multiple-value-setq (s0 s1) - (quick-two-sum s0 c3))))) - (values s0 s1 s2 s3))))))) - + (declare (double-float s0 s1 s2 s3)) + (quick-two-sum s0 c3 c2 c3) + (quick-two-sum s0 c2 c1 s0) + (quick-two-sum c0 c1 c0 s0) + (setf s0 c0) + (setf s1 c1) + (cond ((/= s1 0) + (quick-two-sum s1 s2 s1 c2) + (if (/= s2 0) + (quick-two-sum s2 s3 s2 c3) + (quick-two-sum s1 s2 s1 c3))) + (t + (quick-two-sum s0 s1 s0 c2) + (if (/= s1 0) + (quick-two-sum s1 s2 s1 c3) + (quick-two-sum s0 s1 s0 c3)))) + (values s0 s1 s2 s3))) + (defun renorm-5 (c0 c1 c2 c3 c4) - (declare (double-float c0 c1 c2 c3) + (declare (double-float c0 c1 c2 c3 c4) (optimize (speed 3) (safety 0))) - (let ((s2 0d0) + (let ((s0 0d0) + (s1 0d0) + (s2 0d0) (s3 0d0)) - (declare (double-float s2 s3)) - (multiple-value-bind (s0 c4) - (quick-two-sum c3 c4) - (multiple-value-bind (s0 c3) - (quick-two-sum c2 s0) - (multiple-value-bind (s0 c2) - (quick-two-sum c1 s0) - (multiple-value-bind (c0 c1) - (quick-two-sum c0 s0) - (let ((s0 c0) - (s1 c1)) - (declare (double-float s0 s1)) - (multiple-value-setq (s0 s1) - (quick-two-sum c0 c1)) - (cond ((/= s1 0) - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c2)) - (cond ((/= s2 0) - (multiple-value-setq (s2 s3) - (quick-two-sum s2 c3)) - (if (/= s3 0) - (incf s3 c4) - (incf s2 c4))) - (t - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c3)) - (if (/= s2 0) - (multiple-value-setq (s2 s3) - (quick-two-sum s2 c4)) - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c4)))))) - (t - (multiple-value-setq (s0 s1) - (quick-two-sum s0 c2)) - (cond ((/= s1 0) - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c3)) - (if (/= s2 0) - (multiple-value-setq (s2 s3) - (quick-two-sum s2 c4)) - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c4)))) - (t - (multiple-value-setq (s0 s1) - (quick-two-sum s0 c3)) - (if (/= s1 0) - (multiple-value-setq (s1 s2) - (quick-two-sum s1 c4)) - (multiple-value-setq (s0 s1) - (quick-two-sum s0 c4))))))) - (values s0 s1 s2 s3)))))))) + (declare (double-float s0 s1 s2 s3)) + (quick-two-sum s0 c4 c3 c4) + (quick-two-sum s0 c3 c2 s0) + (quick-two-sum s0 c2 c1 s0) + (quick-two-sum c0 c1 c0 s0) + (quick-two-sum s0 s1 c0 c1) + (cond ((/= s1 0) + (quick-two-sum s1 s2 s1 c2) + (cond ((/= s2 0) + (quick-two-sum s2 s3 s2 c3) + (if (/= s3 0) + (incf s3 c4) + (incf s2 c4))) + (t + (quick-two-sum s1 s2 s1 c3) + (if (/= s2 0) + (quick-two-sum s2 s3 s2 c4) + (quick-two-sum s1 s2 s1 c4))))) + (t + (quick-two-sum s0 s1 s0 c2) + (cond ((/= s1 0) + (quick-two-sum s1 s2 s1 c3) + (if (/= s2 0) + (quick-two-sum s2 s3 s2 c4) + (quick-two-sum s1 s2 s1 c4))) + (t + (quick-two-sum s0 s1 s0 c3) + (if (/= s1 0) + (quick-two-sum s1 s2 s1 c4) + (quick-two-sum s0 s1 s0 c4)))))) + (values s0 s1 s2 s3)))
(defun make-qd-d (a0 &optional (a1 0d0 a1-p) (a2 0d0) (a3 0d0)) "Create a %quad-double from four double-floats, appropriately @@ -786,17 +772,13 @@ (declare (double-float t0)) (multiple-value-bind (s1 t1) (two-sum q1 p3) - (declare (double-float t1)) + (declare (double-float s1 t1)) (multiple-value-setq (s1 t0) (two-sum s1 t0)) (incf t0 t1) - - (multiple-value-setq (s1 t0) - (quick-two-sum s1 t0)) - (multiple-value-setq (p2 t1) - (quick-two-sum s0 s1)) - (multiple-value-setq (p3 q0) - (quick-two-sum t1 t0)) + (quick-two-sum s1 t0 s1 t0) + (quick-two-sum p2 t1 s0 s1) + (quick-two-sum p3 q0 t1 t0)
(let ((p4 (cl:* 2 (qd-0 a) (qd-3 a))) (p5 (cl:* 2 (qd-1 a) (qd-2 a))))