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

Modified Files: qd.lisp Log Message: MAKE-QD-DD is only defined for CMUCL, so don't unconditionally inline it.

--- /project/oct/cvsroot/oct/qd.lisp 2007/10/16 17:05:13 1.59 +++ /project/oct/cvsroot/oct/qd.lisp 2007/10/18 14:38:11 1.60 @@ -136,6 +136,7 @@ sub-qd-dd sub-qd-d sub-d-qd + #+cmu make-qd-dd abs-qd qd-> @@ -833,6 +834,32 @@ (renorm-4 q0 q1 q2 q3) (%make-qd-d q0 q1 q2 q3)))))))

+(declaim (inline invert-qd)) + +(defun invert-qd(v) ;; a quartic newton iteration for 1/v + ;; to invert v, start with a good guess, x. + ;; let h= 1-v*x ;; h is small + ;; return x+ x*(h+h^2+h^3) . compute h3 in double-float + ;; enough accuracy. + + (let* + ((x (%make-qd-d (cl:/ (qd-0 v)) 0d0 0d0 0d0)) + (h (add-qd-d (neg-qd (mul-qd v x)) + 1.0d0)) + (h2 (mul-qd h h)) ;also use h2 for target + (h3 (* (qd-0 h) (qd-0 h2)))) + (declare (type %quad-double v h h2) + (double-float h3)) + (add-qd x + (mul-qd x + (add-qd h (add-qd-d h2 h3)))))) + +(defun div-qd-i (a b) + (declare (type %quad-double a b) + (optimize (speed 3) + (space 0))) + (mul-qd a (invert-qd b))) + ;; Non-sloppy divide #+(or) (defun div-qd (a b)