Commit 47e6f99e authored by Raymond Toy's avatar Raymond Toy

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.
parent 92630dff
......@@ -47,6 +47,7 @@ that we can always return an integer"
#+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 @@ that we can always return an integer"
(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))
......
......@@ -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)
......
......@@ -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)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment