[oct-cvs] Oct commit: oct qd-fun.lisp qd-rep.lisp qd.lisp

[rtoy]

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

Modified Files:
      Tag: THREE-ARG-BRANCH
	qd-fun.lisp qd-rep.lisp qd.lisp 
Log Message:
qd-rep.lisp:
o Fix typo in compiler macro for sub-d-qd

qd.lisp:
o Use 3-arg versions in div-qd-t to speed things up.  Approximately
  doubles the speed with clisp.

qd-fun.lisp:
o Use 3-arg versions in sqrt-qd to speed things up.  Approximately
  doubles the speed with clisp.


--- /project/oct/cvsroot/oct/qd-fun.lisp	2007/10/18 14:38:56	1.90
+++ /project/oct/cvsroot/oct/qd-fun.lisp	2007/11/07 03:45:41	1.90.2.1
@@ -47,6 +47,7 @@
 
 #+cmu
 (declaim (ext:maybe-inline sqrt-qd))
+#+nil
 (defun sqrt-qd (a)
   "Square root of the (non-negative) quad-float"
   (declare (type %quad-double a)
@@ -79,6 +80,47 @@
 	(setf r (add-qd r (mul-qd r (sub-d-qd half (mul-qd h (sqr-qd r)))))))
       (scale-float-qd (mul-qd r new-a) (ash k -1)))))
 
+(defun sqrt-qd (a)
+  "Square root of the (non-negative) quad-float"
+  (declare (type %quad-double a)
+	   (optimize (speed 3) (space 0)))
+  ;; Perform the following Newton iteration:
+  ;;
+  ;;  x' = x + (1 - a * x^2) * x / 2
+  ;;
+  ;; which converges to 1/sqrt(a).
+  ;;
+  ;; However, there appear to be round-off errors when x is either
+  ;; very large or very small.  So let x = f*2^(2*k).  Then sqrt(x) =
+  ;; 2^k*sqrt(f), and sqrt(f) doesn't have round-off problems.
+  (when (zerop-qd a)
+    (return-from sqrt-qd a))
+  (when (float-infinity-p (qd-0 a))
+    (return-from sqrt-qd a))
+
+  (let* ((k (logandc2 (logb-finite (qd-0 a)) 1))
+	 (new-a (scale-float-qd a (- k))))
+    (assert (evenp k))
+    (let* ((r (make-qd-d (cl:/ (sqrt (the (double-float (0d0))
+				       (qd-0 new-a))))))
+	   (temp (%make-qd-d 0d0 0d0 0d0 0d0))
+	   (half 0.5d0)
+	   (h (mul-qd-d new-a half)))
+      (declare (type %quad-double r))
+      ;; Since we start with double-float precision, three more
+      ;; iterations should give us full accuracy.
+      (dotimes (k 3)
+	#+nil
+	(setf r (add-qd r (mul-qd r (sub-d-qd half (mul-qd h (sqr-qd r))))))
+	(sqr-qd r temp)
+	(mul-qd h temp temp)
+	(sub-d-qd half temp temp)
+	(mul-qd r temp temp)
+	(add-qd r temp r)
+	)
+      (mul-qd r new-a r)
+      (scale-float-qd r (ash k -1)))))
+
 (defun hypot-aux-qd (x y)
   (declare (type %quad-double x y))
   (let ((k (- (logb-finite (max (cl:abs (qd-0 x))
--- /project/oct/cvsroot/oct/qd-rep.lisp	2007/11/07 03:08:26	1.10.2.5
+++ /project/oct/cvsroot/oct/qd-rep.lisp	2007/11/07 03:45:41	1.10.2.6
@@ -299,7 +299,7 @@
 
 #-cmu
 (define-compiler-macro sub-d-qd (a b &optional (c #-cmu (%make-qd-d 0d0 0d0 0d0 0d0)))
-  `(add-d-qd a (neg-qd ,b) ,c))
+  `(add-d-qd ,a (neg-qd ,b) ,c))
 
 #+cmu
 (define-compiler-macro neg-qd (a &optional c)
--- /project/oct/cvsroot/oct/qd.lisp	2007/11/07 03:08:26	1.60.2.6
+++ /project/oct/cvsroot/oct/qd.lisp	2007/11/07 03:45:41	1.60.2.7
@@ -863,6 +863,7 @@
 (defun div-qd (a b &optional (target #-cmu (%make-qd-d 0d0 0d0 0d0 0d0)))
   (div-qd-t a b target))
 
+#+nil
 (defun div-qd-t (a b target)
   (declare (type %quad-double a b)
 	   (optimize (speed 3)
@@ -886,6 +887,33 @@
 	      (renorm-4 q0 q1 q2 q3)
 	    (%store-qd-d target q0 q1 q2 q3)))))))
 
+(defun div-qd-t (a b target)
+  (declare (type %quad-double a b)
+	   (optimize (speed 3)
+		     (space 0))
+	   (inline float-infinity-p)
+	   #+cmu
+	   (ignore target))
+  (let ((a0 (qd-0 a))
+	(b0 (qd-0 b)))
+    (let* ((q0 (cl:/ a0 b0))
+	   (r (%make-qd-d 0d0 0d0 0d0 0d0)))
+      (mul-qd-d b q0 r)
+      (sub-qd a r r)
+      (let* ((q1 (cl:/ (qd-0 r) b0))
+	     (temp (mul-qd-d b q1)))
+	(when (float-infinity-p q0)
+	  (return-from div-qd-t (%store-qd-d target q0 0d0 0d0 0d0)))
+	
+	(sub-qd r temp r)
+	(let ((q2 (cl:/ (qd-0 r) b0)))
+	  (mul-qd-d b q2 temp)
+	  (sub-qd r temp r)
+	  (let ((q3 (cl:/ (qd-0 r) b0)))
+	    (multiple-value-bind (q0 q1 q2 q3)
+		(renorm-4 q0 q1 q2 q3)
+	      (%store-qd-d target q0 q1 q2 q3))))))))
+
 (declaim (inline invert-qd))
 
 (defun invert-qd (v)