This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "OCT: A portable Lisp implementation for quad-double precision floats". The branch, master has been updated via b1d9be467f5f44383a356a4a11a04a64f0f29402 (commit) from bcb0aded33f57a2ea901fdcc43fdfd88e8dff17b (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit b1d9be467f5f44383a356a4a11a04a64f0f29402 Author: Raymond Toy <toy.raymond@gmail.com> Date: Wed Nov 27 09:09:50 2013 -0800 Add implementation for tan(x) using duplication. Includes Lentz' algorithm for evaluating continued fractions and updates the timing test to include tan(x) using duplication. Tests show this is much slower than existing algorithms. diff --git a/qd-extra.lisp b/qd-extra.lisp index 61ce7e4..3670be5 100644 --- a/qd-extra.lisp +++ b/qd-extra.lisp @@ -993,6 +993,66 @@ y (add-qd (mul-qd y ct) (mul-qd x st))) (div-qd y x)))) +(defvar *debug-cf-eval* nil) +(defvar *max-cf-iterations* 10000) + +(defun lentz-qd (bf af) + (let ((tiny-value-count 0)) + (flet ((value-or-tiny (v) + (cond ((zerop-qd v) + (incf tiny-value-count) + (make-qd-d (sqrt least-positive-normalized-double-float))) + (t + v)))) + (let* ((f (value-or-tiny (funcall bf 0))) + (c f) + (d (make-qd-d 0d0)) + (eps +qd-eps+)) + (loop + for j from 1 upto *max-cf-iterations* + for an = (funcall af j) + for bn = (funcall bf j) + do (progn + (setf d (value-or-tiny (add-qd bn (mul-qd an d)))) + (setf c (value-or-tiny (add-qd bn (div-qd an c)))) + (when *debug-cf-eval* + (format t "~&j = ~d~%" j) + (format t " an = ~s~%" an) + (format t " bn = ~s~%" bn) + (format t " c = ~s~%" c) + (format t " d = ~s~%" d)) + (let ((delta (div-qd c d))) + (setf d (div-qd +qd-one+ d)) + (setf f (mul-qd f delta)) + (when *debug-cf-eval* + (format t " dl= ~S (|dl - 1| = ~S)~%" delta (abs (1- delta))) + (format t " f = ~S~%" f)) + (when (qd-<= (abs-qd (sub-qd-d delta 1d0)) eps) + (return-from lentz-qd (values f j tiny-value-count))))) + finally + (error 'simple-error + :format-control "~<Continued fraction failed to converge after ~D iterations.~% Delta = ~S~>" + :format-arguments (list *max-cf-iterations* (/ c d)))))))) + +(defun tan-qd/duplication (x) + ;; tan(x) = 2*tan(x/2)/(1-tan(x/2)^2) + (cond ((< (abs (qd-0 x)) 1d-4) + ;; Continued fraction for tan(x) + ;; + ;; x + ;; tan(x) = ----------------------------- + ;; 1 - K(2*k+1, -z^2, k, 0, inf) + (let ((x2 (neg-qd (sqr-qd x)))) + (div-qd x + (lentz-qd #'(lambda (n) + (octi::make-qd-d (float (+ (* 2 n) 1) 1d0))) + #'(lambda (n) + (declare (ignore n)) + x2))))) + (t + (let* ((tanx/2 (tan-qd/duplication (scale-float-qd x -1)))) + (/ (mul-qd-d tanx/2 2d0) + (sub-d-qd 1d0 (sqr-qd tanx/2))))))) (defun sin-qd/cordic (r) (declare (type %quad-double r)) @@ -1374,17 +1434,22 @@ (format t "atan2-qd/newton~%") (time (dotimes (k n) (declare (fixnum k)) - (setf y (atan2-qd/newton x one)))) + (setf y (octi::atan2-qd/newton x one)))) #+cmu (ext:gc :full t) (format t "atan2-qd/cordic~%") (time (dotimes (k n) (declare (fixnum k)) - (setf y (atan2-qd/cordic x one)))) + (setf y (octi::atan2-qd/cordic x one)))) #+cmu (ext:gc :full t) (format t "atan-qd/duplication~%") (time (dotimes (k n) (declare (fixnum k)) - (setf y (atan-qd/duplication x)))) + (setf y (octi::atan-qd/duplication x)))) + #+cmu (ext:gc :full t) + (format t "atan-qd/taylor~%") + (time (dotimes (k n) + (declare (fixnum k)) + (setf y (octi::atan-qd/taylor x)))) )) ;; (time-tan #c(10w0 0) 10000) ----------------------------------------------------------------------- Summary of changes: qd-extra.lisp | 71 ++++++++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 68 insertions(+), 3 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats