# [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)"
+  ;; 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))
+						(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.
+  ;; 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))
+						(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)