[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. c31b1cdd112f26334b7762014d3afb781917ebda
Raymond Toy
rtoy at common-lisp.net
Tue Mar 29 02:37:47 UTC 2011
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 at 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 at 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
More information about the oct-cvs
mailing list