[oct-cvs] Oct commit: oct qd-fun.lisp
rtoy
rtoy at common-lisp.net
Mon Oct 15 03:26:38 UTC 2007
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)
More information about the oct-cvs
mailing list