[oct-cvs] Oct commit: oct qd-fun.lisp
rtoy
rtoy at common-lisp.net
Sat Oct 13 20:14:46 UTC 2007
Update of /project/oct/cvsroot/oct
In directory clnet:/tmp/cvs-serv3017
Modified Files:
qd-fun.lisp
Log Message:
o New routine to compute x mod pi/2 accurately, using many bits of
2/pi.
o Implement accurate sin and cos routines to use this new routine.
(Not used yet.)
--- /project/oct/cvsroot/oct/qd-fun.lisp 2007/10/10 15:21:47 1.81
+++ /project/oct/cvsroot/oct/qd-fun.lisp 2007/10/13 20:14:45 1.82
@@ -612,7 +612,205 @@
(values (neg-qd s)
(neg-qd c))))))))))))
-
+;; A more accurate implementation of sin and cos. We use a more
+;; accurate argument reduction using 1584 bits of 2/pi. This makes
+;; sin and cos more expensive, of course.
+
+(defun rem-pi/2-int (qd)
+ ;; Compute qd rem pi/2 = k*pi/2+y. So we compute k + y*2/pi =
+ ;; qd*2/pi.
+ ;;
+ ;; First convert qd to 2^e*I. We already have 2/pi in the form
+ ;; 2^-1584*J. Then qd*2/pi = 2^(e-1584)*I*J. Extract out the
+ ;; integer and fractional parts of this. For the integer part we
+ ;; only need it modulo 4, because of the periodicity. For the
+ ;; fractional part, we only need 212 (or so bits of fraction).
+ ;;
+ ;; FIXME: But we don't really need to compute all the bits of I*J.
+ ;; In the product, we really only need the 2 bits to the left of the
+ ;; binary point, and then 212 bits to the right. This doesn't
+ ;; require doing the full multiplication.
+ (multiple-value-bind (i e s)
+ (integer-decode-qd qd)
+ (let* ((exp (- e 1584))
+ (prod (* (* s i) +2/pi-bits+))
+ (int (ash prod exp))
+ (mod (logand int 3))
+ (frac (ldb (byte 212 (- (- exp) 212)) prod))
+ (f (mul-qd (scale-float-qd (rational-to-qd frac) -212)
+ +qd-pi/2+)))
+ ;; We want the remainder part to be <= pi/4 because the trig
+ ;; functions want that. So if the fraction is too big, adjust
+ ;; it, and mod value.
+ (if (<= (abs (qd-0 f)) (/ pi 4))
+ (values mod f)
+ (values (mod (1+ mod) 4) (sub-qd f +qd-pi/2+))))))
+
+(defun accurate-sin-qd (a)
+ "Sin(a)"
+ (declare (type %quad-double a))
+ ;; To compute sin(x), choose integers a, b so that
+ ;;
+ ;; x = s + a * (pi/2) + b*(pi/1024)
+ ;;
+ ;; with |x| <= pi/2048. Using a precomputed table of sin(k*pi/1024)
+ ;; and cos(k*pi/1024), we can compute sin(x) from sin(s) and cos(s).
+ ;;
+ ;; sin(x) = sin(s+k*(pi/1024) + j*pi/2)
+ ;; = sin(s+k*(pi/1024))*cos(j*pi/2)
+ ;; + cos(s+k*(pi/1024))*sin(j*pi/2)
+ ;;
+ ;; sin(s+k*pi/1024) = sin(s)*cos(k*pi/1024)
+ ;; + cos(s)*sin(k*pi/1024)
+ ;;
+ ;; cos(s+k*pi/1024) = cos(s)*cos(k*pi/1024)
+ ;; - sin(s)*sin(k*pi/1024)
+ (when (zerop-qd a)
+ (return-from sin-qd a))
+
+ (multiple-value-bind (j r)
+ (rem-pi/2-int a)
+ (let* ((j (if (= j 3)
+ -1
+ j))
+ (abs-j (abs j)))
+ (multiple-value-bind (k tmp)
+ (divrem-qd r +qd-pi/1024+)
+ (let* ((k (truncate (qd-0 k)))
+ (abs-k (abs k)))
+ (assert (<= abs-j 2))
+ (assert (<= abs-k 256))
+ ;; Compute sin(s) and cos(s)
+ (multiple-value-bind (sin-t cos-t)
+ (sincos-taylor tmp)
+ (multiple-value-bind (s c)
+ (cond ((zerop abs-k)
+ (values sin-t cos-t))
+ (t
+ ;; Compute sin(s+k*pi/1024), cos(s+k*pi/1024)
+ (let ((u (aref +qd-cos-table+ (cl:1- abs-k)))
+ (v (aref +qd-sin-table+ (cl:1- abs-k))))
+ (cond ((plusp k)
+ ;; sin(s) * cos(k*pi/1024)
+ ;; + cos(s) * sin(k*pi/1024)
+ ;;
+ ;; cos(s) * cos(k*pi/1024)
+ ;; - sin(s) * sin(k*pi/1024)
+ (values (add-qd (mul-qd u sin-t)
+ (mul-qd v cos-t))
+ (sub-qd (mul-qd u cos-t)
+ (mul-qd v sin-t))))
+ (t
+ ;; sin(s) * cos(k*pi/1024)
+ ;; - cos(s) * sin(|k|*pi/1024)
+ ;;
+ ;; cos(s) * cos(k*pi/1024)
+ ;; + sin(s) * sin(|k|*pi/1024)
+ (values (sub-qd (mul-qd u sin-t)
+ (mul-qd v cos-t))
+ (add-qd (mul-qd u cos-t)
+ (mul-qd v sin-t))))))))
+ ;;(format t "s = ~/qd::qd-format/~%" s)
+ ;;(format t "c = ~/qd::qd-format/~%" c)
+ ;; sin(x) = sin(s+k*pi/1024) * cos(j*pi/2)
+ ;; + cos(s+k*pi/1024) * sin(j*pi/2)
+ (cond ((zerop abs-j)
+ ;; cos(j*pi/2) = 1, sin(j*pi/2) = 0
+ s)
+ ((= j 1)
+ ;; cos(j*pi/2) = 0, sin(j*pi/2) = 1
+ c)
+ ((= j -1)
+ ;; cos(j*pi/2) = 0, sin(j*pi/2) = -1
+ (neg-qd c))
+ (t
+ ;; cos(j*pi/2) = -1, sin(j*pi/2) = 0
+ (neg-qd s))))))))))
+
+(defun accurate-cos-qd (a)
+ "Cos(a)"
+ ;; Just like sin-qd, but for cos.
+ (declare (type %quad-double a))
+ ;; To compute sin(x), choose integers a, b so that
+ ;;
+ ;; x = s + a * (pi/2) + b*(pi/1024)
+ ;;
+ ;; with |x| <= pi/2048. Using a precomputed table of sin(k*pi/1024)
+ ;; and cos(k*pi/1024), we can compute sin(x) from sin(s) and cos(s).
+ ;;
+ ;; sin(x) = sin(s+k*(pi/1024) + j*pi/2)
+ ;; = sin(s+k*(pi/1024))*cos(j*pi/2)
+ ;; + cos(s+k*(pi/1024))*sin(j*pi/2)
+ ;;
+ ;; sin(s+k*pi/1024) = sin(s)*cos(k*pi/1024)
+ ;; + cos(s)*sin(k*pi/1024)
+ ;;
+ ;; cos(s+k*pi/1024) = cos(s)*cos(k*pi/1024)
+ ;; - sin(s)*sin(k*pi/1024)
+ (when (zerop-qd a)
+ (return-from cos-qd +qd-one+))
+
+ (multiple-value-bind (j r)
+ (rem-pi/2-int a)
+ (let* ((j (if (= j 3)
+ -1
+ j))
+ (abs-j (abs j)))
+ (multiple-value-bind (k tmp)
+ (divrem-qd r +qd-pi/1024+)
+ (let* ((k (truncate (qd-0 k)))
+ (abs-k (abs k)))
+ (assert (<= abs-j 2))
+ (assert (<= abs-k 256))
+ ;; Compute sin(s) and cos(s)
+ (multiple-value-bind (sin-t cos-t)
+ (sincos-taylor tmp)
+ (multiple-value-bind (s c)
+ (cond ((zerop abs-k)
+ (values sin-t cos-t))
+ (t
+ ;; Compute sin(s+k*pi/1024), cos(s+k*pi/1024)
+ (let ((u (aref +qd-cos-table+ (cl:1- abs-k)))
+ (v (aref +qd-sin-table+ (cl:1- abs-k))))
+ (cond ((plusp k)
+ ;; sin(s) * cos(k*pi/1024)
+ ;; + cos(s) * sin(k*pi/1024)
+ ;;
+ ;; cos(s) * cos(k*pi/1024)
+ ;; - sin(s) * sin(k*pi/1024)
+ (values (add-qd (mul-qd u sin-t)
+ (mul-qd v cos-t))
+ (sub-qd (mul-qd u cos-t)
+ (mul-qd v sin-t))))
+ (t
+ ;; sin(s) * cos(k*pi/1024)
+ ;; - cos(s) * sin(|k|*pi/1024)
+ ;;
+ ;; cos(s) * cos(k*pi/1024)
+ ;; + sin(s) * sin(|k|*pi/1024)
+ (values (sub-qd (mul-qd u sin-t)
+ (mul-qd v cos-t))
+ (add-qd (mul-qd u cos-t)
+ (mul-qd v sin-t))))))))
+ #+nil
+ (progn
+ (format t "s = ~/qd::qd-format/~%" s)
+ (format t "c = ~/qd::qd-format/~%" c))
+ ;; sin(x) = sin(s+k*pi/1024) * cos(j*pi/2)
+ ;; + cos(s+k*pi/1024) * sin(j*pi/2)
+ (cond ((zerop abs-j)
+ ;; cos(j*pi/2) = 1, sin(j*pi/2) = 0
+ c)
+ ((= j 1)
+ ;; cos(j*pi/2) = 0, sin(j*pi/2) = 1
+ (neg-qd s))
+ ((= j -1)
+ ;; cos(j*pi/2) = 0, sin(j*pi/2) = -1
+ s)
+ (t
+ ;; cos(j*pi/2) = -1, sin(j*pi/2) = 0
+ (neg-qd c))))))))))
+
(defun atan2-qd/newton (y x)
(declare (type %quad-double y x)
#+nil
More information about the oct-cvs
mailing list