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 c31b1cdd112f26334b7762014d3afb781917ebda (commit) via 88cff63cfb4996e2e90e499afd34e1fe16ecc179 (commit) from a44327a2532c5afc3cbcb1194aa1f15b6dccd2cb (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 c31b1cdd112f26334b7762014d3afb781917ebda Author: Raymond Toy <toy.raymond@gmail.com> Date: Mon Mar 28 22:37:32 2011 -0400 More accurate incomplete-gamma function, add debugging to lentz, and some random clean ups. o Add *DEBUG-CF-EVAL* to enable debugging prints in LENTZ. o Modify LENTZ to terminate with an error if *MAX-CF-ITERATIONS* is reached. o Modify LENTZ to return the function value, the number of iterations, and the number of times a zero value had to be replaced. o Adjust cf-incomplete-gamma and cf-incomplete-gamma-tail not to signal overflow prematurely when calculating z^a*exp(-z). o Fix doc bug in reference for continued fraction for (original) cf-incomplete-gamma. o Add new version of cf-incomplete-gamma using a different continued fraction. This appears to converge faster and to be more accurate than the original, especially for points near the negative real axis. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index 78b5a99..e020a91 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -177,32 +177,62 @@ ;; b0 + ---- ---- ---- ;; b1 + b2 + b3 + ;; + +(defvar *debug-cf-eval* + nil + "When true, enable some debugging prints when evaluating a + continued fraction.") + +;; Max number of iterations allowed when evaluating the continued +;; fraction. When this is reached, we assume that the continued +;; fraction did not converge. +(defvar *max-cf-iterations* + 10000 + "Max number of iterations allowed when evaluating the continued + fraction. When this is reached, we assume that the continued + fraction did not converge.") + (defun lentz (bf af) - (flet ((value-or-tiny (v) - (if (zerop v) - (etypecase v - ((or double-float cl:complex) - least-positive-normalized-double-float) - ((or qd-real qd-complex) - (make-qd least-positive-normalized-double-float))) - v))) - (let* ((f (value-or-tiny (funcall bf 0))) - (c f) - (d 0) - (eps (epsilon f))) - (loop - for j from 1 - for an = (funcall af j) - for bn = (funcall bf j) - do (progn - (setf d (value-or-tiny (+ bn (* an d)))) - (setf c (value-or-tiny (+ bn (/ an c)))) - (setf d (/ d)) - (let ((delta (* c d))) - (setf f (* f delta)) - (when (<= (abs (- delta 1)) eps) - (return))))) - f))) + (let ((tiny-value-count 0)) + (flet ((value-or-tiny (v) + (if (zerop v) + (progn + (incf tiny-value-count) + (etypecase v + ((or double-float cl:complex) + least-positive-normalized-double-float) + ((or qd-real qd-complex) + (make-qd least-positive-normalized-double-float)))) + v))) + (let* ((f (value-or-tiny (funcall bf 0))) + (c f) + (d 0) + (eps (epsilon f))) + (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 (+ bn (* an d)))) + (setf c (value-or-tiny (+ bn (/ 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 (/ c d))) + (setf d (/ d)) + (setf f (* f delta)) + (when *debug-cf-eval* + (format t " dl= ~S~%" delta) + (format t " f = ~S~%" f)) + (when (<= (abs (- delta 1)) eps) + (return-from lentz (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)))))))) ;; Continued fraction for erf(b): ;; @@ -240,8 +270,13 @@ :function-name 'cf-incomplete-gamma-tail :format-arguments (list 'z z) :format-control "Argument ~S should not be on the negative real axis: ~S")) - (/ (* (expt z a) - (exp (- z))) + (/ (handler-case (* (expt z a) + (exp (- z))) + (arithmetic-error () + ;; z^a*exp(-z) can overflow prematurely. In this case, use + ;; the equivalent exp(a*log(z)-z). We don't use this latter + ;; form because it has more roundoff error than the former. + (exp (- (* a (log z)) z)))) (let ((z-a (- z a))) (lentz #'(lambda (n) (+ n n 1 z-a)) @@ -251,18 +286,23 @@ ;; Incomplete gamma function: ;; integrate(x^(a-1)*exp(-x), x, 0, z) ;; -;; The continued fraction, valid for all z except the negative real -;; axis: +;; The continued fraction, valid for all z: ;; ;; b[n] = n - 1 + z + a ;; a[n] = -z*(a + n) ;; -;; See http://functions.wolfram.com/06.06.10.0005.01. We modified the +;; See http://functions.wolfram.com/06.06.10.0007.01. We modified the ;; continued fraction slightly and discarded the first quotient from ;; the fraction. +#+nil (defun cf-incomplete-gamma (a z) - (/ (* (expt z a) - (exp (- z))) + (/ (handler-case (* (expt z a) + (exp (- z))) + (arithmetic-error () + ;; z^a*exp(-z) can overflow prematurely. In this case, use + ;; the equivalent exp(a*log(z)-z). We don't use this latter + ;; form because it has more roundoff error than the former. + (exp (- (* a (log z)) z)))) (let ((za1 (+ z a 1))) (- a (/ (* a z) (lentz #'(lambda (n) @@ -270,6 +310,35 @@ #'(lambda (n) (- (* z (+ a n)))))))))) +;; Incomplete gamma function: +;; integrate(x^(a-1)*exp(-x), x, 0, z) +;; +;; The continued fraction, valid for all z: +;; +;; b[n] = a + n +;; a[n] = -(a+n/2)*z if n odd +;; n/2*z if n even +;; +;; See http://functions.wolfram.com/06.06.10.0009.01. +;; +;; Some experiments indicate that this converges faster than the above +;; and is actually quite a bit more accurate, expecially near the +;; negative real axis. +(defun cf-incomplete-gamma (a z) + (/ (handler-case (* (expt z a) + (exp (- z))) + (arithmetic-error () + ;; z^a*exp(-z) can overflow prematurely. In this case, use + ;; the equivalent exp(a*log(z)-z). We don't use this latter + ;; form because it has more roundoff error than the former. + (exp (- (* a (log z)) z)))) + (lentz #'(lambda (n) + (+ n a)) + #'(lambda (n) + (if (evenp n) + (* (ash n -1) z) + (- (* (+ a (ash n -1)) z))))))) + ;; Series expansion for incomplete gamma. Intended for |a|<1 and ;; |z|<1. The series is ;; @@ -283,8 +352,6 @@ when (< (abs term) (* (abs sum) eps)) return (* sum (expt z a))))) - - ;; Tail of the incomplete gamma function. (defun incomplete-gamma-tail (a z) "Tail of the incomplete gamma function defined by: commit 88cff63cfb4996e2e90e499afd34e1fe16ecc179 Author: Raymond Toy <toy.raymond@gmail.com> Date: Mon Mar 28 09:08:37 2011 -0400 Fix typo in #q() number in fresnel-s.2q test. diff --git a/rt-tests.lisp b/rt-tests.lisp index c0b59fb..66aab0c 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1279,7 +1279,7 @@ nil) (rt:deftest fresnel-s.2q - (let* ((z #q(1q-3 1q-3)) + (let* ((z #q(#q1q-3 #q1q-3)) (s (fresnel-s z)) (true (fresnel-s-series z))) (check-accuracy 212 s true)) ----------------------------------------------------------------------- Summary of changes: qd-gamma.lisp | 135 ++++++++++++++++++++++++++++++++++++++++++-------------- rt-tests.lisp | 2 +- 2 files changed, 102 insertions(+), 35 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats