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)