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

Modified Files: qd-fun.lisp Log Message: o Fix some typos accurate-sin-qd and accurate-cos-qd. o Adjust code in accurate-sin-qd and accurate-cos-qd to handle values of 0 <= j <= 3, instead of -1 <= j <= 2. o Add accurate-sincos-qd.

--- /project/oct/cvsroot/oct/qd-fun.lisp 2007/10/13 20:14:45 1.82 +++ /project/oct/cvsroot/oct/qd-fun.lisp 2007/10/15 03:26:38 1.83 @@ -666,66 +666,62 @@ ;; 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)) + (return-from accurate-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)) + (multiple-value-bind (k tmp) + (divrem-qd r +qd-pi/1024+) + (let* ((k (truncate (qd-0 k))) + (abs-k (abs k))) + (assert (<= 0 j 3)) + (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 - ;; cos(j*pi/2) = -1, sin(j*pi/2) = 0 - (neg-qd s)))))))))) + ;; 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 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 2) + ;; cos(j*pi/2) = -1, sin(j*pi/2) = 0 + (neg-qd s)) + ((= j 3) + ;; cos(j*pi/2) = 0, sin(j*pi/2) = -1 + (neg-qd c)))))))))

(defun accurate-cos-qd (a) "Cos(a)" @@ -748,68 +744,130 @@ ;; 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+)) + (return-from accurate-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) + (multiple-value-bind (k tmp) + (divrem-qd r +qd-pi/1024+) + (let* ((k (truncate (qd-0 k))) + (abs-k (abs k))) + (assert (<= 0 j 3)) + (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 - ;; cos(j*pi/2) = -1, sin(j*pi/2) = 0 - (neg-qd c)))))))))) + ;; 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 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 2) + ;; cos(j*pi/2) = -1, sin(j*pi/2) = 0 + (neg-qd c)) + ((= j 3) + ;; cos(j*pi/2) = 0, sin(j*pi/2) = -1 + s)))))))) + +(defun accurate-sincos-qd (a) + (declare (type %quad-double a)) + (when (zerop-qd a) + (return-from accurate-sincos-qd + (values +qd-zero+ + +qd-one+))) + + (multiple-value-bind (j r) + (rem-pi/2-int a) + (multiple-value-bind (k tmp) + (divrem-qd tmp +qd-pi/1024+) + (let* ((k (truncate (qd-0 k))) + (abs-j (abs j)) + (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 + (values s c)) + ((= j 1) + ;; cos(j*pi/2) = 0, sin(j*pi/2) = 1 + (values c (neg-qd s))) + ((= j 2) + ;; cos(j*pi/2) = -1, sin(j*pi/2) = 0 + (values (neg-qd s) + (neg-qd c))) + ((= j 3) + ;; cos(j*pi/2) = 0, sin(j*pi/2) = -1 + (values (neg-qd c) s))))))))))

(defun atan2-qd/newton (y x) (declare (type %quad-double y x)