Update of /project/oct/cvsroot/oct In directory clnet:/tmp/cvs-serv365

Modified Files: qd.lisp Log Message: o Replace THREE-SUM with a macro to make sure it's inlined everywhere. o Update code for new THREE-SUM macro.

--- /project/oct/cvsroot/oct/qd.lisp 2007/09/16 14:23:24 1.49 +++ /project/oct/cvsroot/oct/qd.lisp 2007/09/17 03:08:25 1.50 @@ -53,17 +53,25 @@ (c::two-sum ,x ,y)))

-(defun three-sum (a b c) - (declare (double-float a b c) - (optimize (speed 3))) - (let* ((t1 0d0) - (t2 t1) - (t3 t1)) - (declare (double-float t1 t2 t3)) - (two-sum t1 t2 a b) - (two-sum a t3 c t1) - (two-sum b c t2 t3) - (values a b c))) +(defmacro three-sum (s0 s1 s2 a b c) + (let ((aa (gensym)) + (bb (gensym)) + (cc (gensym)) + (t1 (gensym)) + (t2 (gensym)) + (t3 (gensym))) + `(let* ((,aa ,a) + (,bb ,b) + (,cc ,c) + (,t1 0d0) + (,t2 ,t1) + (,t3 ,t1)) + (declare (double-float ,t1 ,t2 ,t3)) + (two-sum ,t1 ,t2 ,aa ,bb) + (two-sum ,s0 ,t3 ,cc ,t1) + (two-sum ,s1 ,s2 ,t2 ,t3)))) + +

(defun three-sum2 (a b c) (declare (double-float a b c) @@ -310,17 +318,20 @@ (t0 s0) (s1 s0) (t1 s0) + (s2 s0) (s3 s0)) - (declare (double-float s0 s1 s3 t0 t1)) + (declare (double-float s0 s1 s3 t0 t1 s2)) (two-sum s0 t1 (qd-0 a) (kernel:double-double-hi b)) (two-sum s1 t1 (qd-1 a) (kernel:double-double-lo b)) (two-sum s1 t0 s1 t0) - (multiple-value-bind (s2 t0 t1) - (three-sum (qd-2 a) t0 t1) - (two-sum s3 t0 t0 (qd-3 a)) - (let ((t0 (cl:+ t0 t1))) - (multiple-value-call #'%make-qd-d - (renorm-5 s0 s1 s2 s3 t0)))))) + (three-sum s2 t0 t1 t0 t0 (qd-3 a)) + (two-sum s3 t0 t0 (qd-3 a)) + (let ((t0 (cl:+ t0 t1))) + (declare (double-float t0)) + (multiple-value-bind (a0 a1 a2 a3) + (renorm-5 s0 s1 s2 s3 t0) + (declare (double-float a0 a1 a2 a3)) + (%make-qd-d a0 a1 a2 a3)))))

#+cmu (defun add-dd-qd (a b) @@ -400,18 +411,17 @@ (t2 (cl:+ w2 u2)) (t3 (cl:+ w3 u3))) (two-sum s1 t0 s1 t0) - (multiple-value-bind (s2 t0 t1) - (three-sum s2 t0 t1) - (multiple-value-bind (s3 t0) - (three-sum2 s3 t0 t2) - (declare (double-float t0)) - (setf t0 (cl:+ t0 t1 t3)) - ;; Renormalize - (multiple-value-setq (s0 s1 s2 s3) - (renorm-5 s0 s1 s2 s3 t0)) - (if (and (zerop a0) (zerop b0)) - (%make-qd-d (+ a0 b0) 0d0 0d0 0d0) - (%make-qd-d s0 s1 s2 s3))))))))))))) + (three-sum s2 t0 t1 s2 t0 t1) + (multiple-value-bind (s3 t0) + (three-sum2 s3 t0 t2) + (declare (double-float t0)) + (setf t0 (cl:+ t0 t1 t3)) + ;; Renormalize + (multiple-value-setq (s0 s1 s2 s3) + (renorm-5 s0 s1 s2 s3 t0)) + (if (and (zerop a0) (zerop b0)) + (%make-qd-d (+ a0 b0) 0d0 0d0 0d0) + (%make-qd-d s0 s1 s2 s3))))))))))))

(defun neg-qd (a) (declare (type %quad-double a)) @@ -471,17 +481,16 @@ (s2 p0)) (declare (double-float s0 s1 s2 p3)) (two-sum s1 s2 q0 p1) - (multiple-value-bind (s2 q1 p2) - (three-sum s2 q1 p2) - (multiple-value-bind (q1 q2) - (three-sum2 q1 q2 p3) - (let ((s3 q1) - (s4 (cl:+ q2 p2))) - (multiple-value-bind (s0 s1 s2 s3) - (renorm-5 s0 s1 s2 s3 s4) - (if (zerop s0) - (%make-qd-d (float-sign p0 0d0) 0d0 0d0 0d0) - (%make-qd-d s0 s1 s2 s3))))))))))) + (three-sum s2 q1 p2 s2 q1 p2) + (multiple-value-bind (q1 q2) + (three-sum2 q1 q2 p3) + (let ((s3 q1) + (s4 (cl:+ q2 p2))) + (multiple-value-bind (s0 s1 s2 s3) + (renorm-5 s0 s1 s2 s3 s4) + (if (zerop s0) + (%make-qd-d (float-sign p0 0d0) 0d0 0d0 0d0) + (%make-qd-d s0 s1 s2 s3))))))))))

;; a0 * b0 0 ;; a0 * b1 1 @@ -597,14 +606,11 @@ (multiple-value-bind (p5 q5) (two-prod a2 b0) ;; Start accumulation - (multiple-value-setq (p1 p2 q0) - (three-sum p1 p2 q0)) + (three-sum p1 p2 q0 p1 p2 q0)

;; six-three-sum of p2, q1, q2, p3, p4, p5 - (multiple-value-setq (p2 q1 q2) - (three-sum p2 q1 q2)) - (multiple-value-setq (p3 p4 p5) - (three-sum p3 p4 p5)) + (three-sum p2 q1 q2 p2 q1 q2) + (three-sum p3 p4 p5 p3 p4 p5) ;; Compute (s0,s1,s2) = (p2,q1,q2) + (p3,p4,p5) (let* ((s0 0d0) (s1 s0) @@ -779,8 +785,9 @@ (two-sum p3 p4 p3 t0) (incf p4 (cl:+ q0 t1))

- (multiple-value-call #'%make-qd-d - (renorm-5 p0 p1 p2 p3 p4))))))))) + (multiple-value-bind (a0 a1 a2 a3) + (renorm-5 p0 p1 p2 p3 p4) + (%make-qd-d a0 a1 a2 a3)))))))))

(defun div-qd (a b)