Update of /project/oct/cvsroot/oct
In directory clnet:/tmp/cvs-serv8050
Modified Files:
qd-dd.lisp qd-package.lisp qd.lisp
Log Message:
Make TWO-SUM a macro, just like we did for QUICK-TWO-SUM.
All RT tests pass on CMUCL and Allegro.
qd-package.lisp:
o Don't import C::TWO-SUM anymore.
qd-dd.lisp:
o Make TWO-SUM a macro.
qd.lisp
o Add TWO-SUM macro for CMUCL (which just calls C::TWO-SUM).
o Update all uses of TWO-SUM to use the macro appropriately.
--- /project/oct/cvsroot/oct/qd-dd.lisp 2007/09/16 02:46:24 1.9
+++ /project/oct/cvsroot/oct/qd-dd.lisp 2007/09/16 05:04:04 1.10
@@ -52,6 +52,7 @@
(setf ,s (+ ,a ,b))
(setf ,e (- ,b (- ,s ,a))))))
+#||
(declaim (inline two-sum))
(defun two-sum (a b)
"Computes fl(a+b) and err(a+b)"
@@ -63,21 +64,21 @@
(- b v))))
(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))
+ (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))))))))
+ (declare (double-float ,v))
+ (setf ,e (+ (- ,a (- ,s ,v))
+ (- ,b ,v)))))))
(declaim (inline two-prod))
(declaim (inline split))
--- /project/oct/cvsroot/oct/qd-package.lisp 2007/09/16 02:46:24 1.38
+++ /project/oct/cvsroot/oct/qd-package.lisp 2007/09/16 05:04:05 1.39
@@ -97,7 +97,6 @@
#:make-qd-dd)
#+cmu
(:import-from #:c
- #:two-sum
#:two-prod
#:two-sqr))
--- /project/oct/cvsroot/oct/qd.lisp 2007/09/16 02:46:24 1.47
+++ /project/oct/cvsroot/oct/qd.lisp 2007/09/16 05:04:05 1.48
@@ -46,6 +46,14 @@
(declaim (inline three-sum three-sum2))
;; Internal routines for implementing quad-double.
+
+#+cmu
+(defmacro two-sum (s e x y)
+ `(multiple-value-setq (,s ,e)
+ (c::two-sum ,x ,y)))
+
+
+#+nil
(defun three-sum (a b c)
(declare (double-float a b c)
(optimize (speed 3)))
@@ -57,6 +65,19 @@
(two-sum t2 t3)
(values a b c)))))
+(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)))
+
+#+nil
(defun three-sum2 (a b c)
(declare (double-float a b c)
(optimize (speed 3)))
@@ -66,6 +87,17 @@
(two-sum c t1)
(values a (cl:+ t2 t3) c))))
+(defun three-sum2 (a b c)
+ (declare (double-float a b c)
+ (optimize (speed 3)))
+ (let* ((t1 0d0)
+ (t2 t1)
+ (t3 t1))
+ (two-sum t1 t2 a b)
+ (two-sum a t3 c t1)
+ (values a (cl:+ t2 t3) c)))
+
+
;; Not needed????
#+nil
(declaim (inline quick-three-accum))
@@ -137,11 +169,6 @@
`(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
@@ -267,22 +294,24 @@
(double-float b)
(optimize (speed 3)
(space 0)))
- (multiple-value-bind (c0 e)
- (two-sum (qd-0 a) b)
+ (let* ((c0 0d0)
+ (e c0)
+ (c1 c0)
+ (c2 c0)
+ (c3 c0))
+ (declare (double-float e c0 c1 c2 c3))
+ (two-sum c0 e (qd-0 a) b)
#+cmu
(when (float-infinity-p c0)
(return-from add-qd-d (%make-qd-d c0 0d0 0d0 0d0)))
- (multiple-value-bind (c1 e)
- (two-sum (qd-1 a) e)
- (multiple-value-bind (c2 e)
- (two-sum (qd-2 a) e)
- (multiple-value-bind (c3 e)
- (two-sum (qd-3 a) e)
- (multiple-value-bind (r0 r1 r2 r3)
- (renorm-5 c0 c1 c2 c3 e)
- (if (and (zerop (qd-0 a)) (zerop b))
- (%make-qd-d c0 0d0 0d0 0d0)
- (%make-qd-d r0 r1 r2 r3))))))))
+ (two-sum c1 e (qd-1 a) e)
+ (two-sum c2 e (qd-2 a) e)
+ (two-sum c3 e (qd-3 a) e)
+ (multiple-value-bind (r0 r1 r2 r3)
+ (renorm-5 c0 c1 c2 c3 e)
+ (if (and (zerop (qd-0 a)) (zerop b))
+ (%make-qd-d c0 0d0 0d0 0d0)
+ (%make-qd-d r0 r1 r2 r3)))))
(defun add-d-qd (a b)
(declare (double-float a)
@@ -297,19 +326,21 @@
(double-double-float b)
(optimize (speed 3)
(space 0)))
- (multiple-value-bind (s0 t0)
- (two-sum (qd-0 a) (kernel:double-double-hi b))
- (multiple-value-bind (s1 t1)
- (two-sum (qd-1 a) (kernel:double-double-lo b))
- (multiple-value-bind (s1 t0)
- (two-sum s1 t0)
- (multiple-value-bind (s2 t0 t1)
- (three-sum (qd-2 a) t0 t1)
- (multiple-value-bind (s3 t0)
- (two-sum t0 (qd-3 a))
- (let ((t0 (cl:+ t0 t1)))
- (multiple-value-call #'%make-qd-d
- (renorm-5 s0 s1 s2 s3 t0)))))))))
+ (let* ((s0 0d0)
+ (t0 s0)
+ (s1 s0)
+ (t1 s0)
+ (s3 s0))
+ (declare (double-float s0 s1 s3 t0 t1))
+ (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))))))
#+cmu
(defun add-dd-qd (a b)
@@ -386,20 +417,19 @@
(t1 (cl:+ w1 u1))
(t2 (cl:+ w2 u2))
(t3 (cl:+ w3 u3)))
- (multiple-value-bind (s1 t0)
- (two-sum 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))))))))))))))
+ (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)))))))))))))
(defun neg-qd (a)
(declare (type %quad-double a))
@@ -451,50 +481,23 @@
(two-prod (qd-1 a) b)
(multiple-value-bind (p2 q2)
(two-prod (qd-2 a) b)
- (let ((p3 (cl:* (qd-3 a) b))
- (s0 p0))
- #+nil
- (format t "q0 p1 = ~A ~A~%" q0 p1)
- (multiple-value-bind (s1 s2)
- (two-sum q0 p1)
- #+nil
- (format t "s1 s2 = ~A ~A~%" s1 s2)
- #+nil
- (format t "s2 q1 p2 = ~A ~A ~A~%"
- s2 q1 p2)
- (multiple-value-bind (s2 q1 p2)
- (three-sum s2 q1 p2)
- #+nil
- (format t "s2,q1,p2 = ~A ~A ~A~%"
- s2 q1 p2)
- #+nil
- (format t "q1 q2 p3 = ~A ~A ~A~%"
- q1 q2 p3)
- (multiple-value-bind (q1 q2)
- (three-sum2 q1 q2 p3)
- #+nil
- (format t "q1 q2, p3 = ~A ~A ~A~%" q1 q2 p3)
- (let ((s3 q1)
- (s4 (cl:+ q2 p2)))
- #+nil
- (progn
- (format t "~a~%" s0)
- (format t "~a~%" s1)
- (format t "~a~%" s2)
- (format t "~a~%" s3)
- (format t "~a~%" s4))
- (multiple-value-bind (s0 s1 s2 s3)
- (renorm-5 s0 s1 s2 s3 s4)
- #+nil
- (progn
- (format t "~a~%" s0)
- (format t "~a~%" s1)
- (format t "~a~%" s2)
- (format t "~a~%" s3)
- (format t "~a~%" s4))
- (if (zerop s0)
- (%make-qd-d (float-sign p0 0d0) 0d0 0d0 0d0)
- (%make-qd-d s0 s1 s2 s3))))))))))))
+ (let* ((p3 (cl:* (qd-3 a) b))
+ (s0 p0)
+ (s1 p0)
+ (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)))))))))))
;; a0 * b0 0
;; a0 * b1 1
@@ -617,34 +620,35 @@
(multiple-value-setq (p3 p4 p5)
(three-sum p3 p4 p5))
;; Compute (s0,s1,s2) = (p2,q1,q2) + (p3,p4,p5)
- (multiple-value-bind (s0 t0)
- (two-sum p2 p3)
- (multiple-value-bind (s1 t1)
- (two-sum q1 p4)
- (let ((s2 (cl:+ q2 p5)))
- (declare (double-float s2))
- (multiple-value-bind (s1 t0)
- (two-sum s1 t0)
- (declare (double-float s1))
- (incf s2 (cl:+ t0 t1))
- ;; O(eps^3) order terms. This is the sloppy
- ;; multiplication version. Should we use
- ;; the precise version? It's significantly
- ;; more complex.
+ (let* ((s0 0d0)
+ (s1 s0)
+ (t0 s0)
+ (t1 s0))
+ (declare (double-float s0 s1 t0 t1))
+ (two-sum s0 t0 p2 p3)
+ (two-sum s1 t1 q1 p4)
+ (let ((s2 (cl:+ q2 p5)))
+ (declare (double-float s2))
+ (two-sum s1 t0 s1 t0)
+ (incf s2 (cl:+ t0 t1))
+ ;; O(eps^3) order terms. This is the sloppy
+ ;; multiplication version. Should we use
+ ;; the precise version? It's significantly
+ ;; more complex.
- (incf s1 (cl:+ (cl:* a0 b3)
- (cl:* a1 b2)
- (cl:* a2 b1)
- (cl:* a3 b0)
- q0 q3 q4 q5))
- #+nil
- (format t "p0,p1,s0,s1,s2 = ~a ~a ~a ~a ~a~%"
- p0 p1 s0 s1 s2)
- (multiple-value-bind (r0 r1 s0 s1)
- (renorm-5 p0 p1 s0 s1 s2)
- (if (zerop r0)
- (%make-qd-d p0 0d0 0d0 0d0)
- (%make-qd-d r0 r1 s0 s1))))))))))))))))
+ (incf s1 (cl:+ (cl:* a0 b3)
+ (cl:* a1 b2)
+ (cl:* a2 b1)
+ (cl:* a3 b0)
+ q0 q3 q4 q5))
+ #+nil
+ (format t "p0,p1,s0,s1,s2 = ~a ~a ~a ~a ~a~%"
+ p0 p1 s0 s1 s2)
+ (multiple-value-bind (r0 r1 s0 s1)
+ (renorm-5 p0 p1 s0 s1 s2)
+ (if (zerop r0)
+ (%make-qd-d p0 0d0 0d0 0d0)
+ (%make-qd-d r0 r1 s0 s1))))))))))))))
;; This is the non-sloppy version. I think this works just fine, but
;; since qd defaults to the sloppy multiplication version, we do the
@@ -760,45 +764,37 @@
(two-prod (cl:* 2 (qd-0 a)) (qd-2 a))
(multiple-value-bind (p3 q3)
(two-sqr (qd-1 a))
- (multiple-value-setq (p1 q0)
- (two-sum q0 p1))
- (multiple-value-setq (q0 q1)
- (two-sum q0 q1))
- (multiple-value-setq (p2 p3)
- (two-sum p2 p3))
-
- (multiple-value-bind (s0 t0)
- (two-sum q0 p2)
- (declare (double-float t0))
- (multiple-value-bind (s1 t1)
- (two-sum q1 p3)
- (declare (double-float s1 t1))
- (multiple-value-setq (s1 t0)
- (two-sum s1 t0))
- (incf t0 t1)
- (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))))
- (declare (double-float p4))
- (multiple-value-setq (p4 p5)
- (two-sum p4 p5))
- (multiple-value-setq (q2 q3)
- (two-sum q2 q3))
-
- (multiple-value-setq (t0 t1)
- (two-sum p4 q2))
-
- (incf t1 (cl:+ p5 q3))
-
- (multiple-value-setq (p3 p4)
- (two-sum p3 t0))
- (incf p4 (cl:+ q0 t1))
+ (two-sum p1 q0 q0 p1)
+ (two-sum q0 q1 q0 q1)
+ (two-sum p2 p3 p2 p3)
+
+ (let* ((s0 0d0)
+ (t0 s0)
+ (s1 s0)
+ (t1 s0))
+ (declare (double-float s0 s1 t0 t1))
+ (two-sum s0 t0 q0 p2)
+ (two-sum s1 t1 q1 p3)
+ (two-sum s1 t0 s1 t0)
+ (incf t0 t1)
+ (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))))
+ (declare (double-float p4 p5))
+ (two-sum p4 p5 p4 p5)
+ (two-sum q2 q3 q2 q3)
+ (two-sum t0 t1 p4 q2)
+
+ (incf t1 (cl:+ p5 q3))
- (multiple-value-call #'%make-qd-d
- (renorm-5 p0 p1 p2 p3 p4))))))))))
+ (two-sum p3 p4 p3 t0)
+ (incf p4 (cl:+ q0 t1))
+
+ (multiple-value-call #'%make-qd-d
+ (renorm-5 p0 p1 p2 p3 p4)))))))))
(defun div-qd (a b)
