# [oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. af9fc6fcdb824757c5aeeadc6ccae2098d243313

Raymond Toy rtoy at common-lisp.net
Tue Mar 29 13:59:14 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

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 -----------------------------------------------------------------
Author: Raymond Toy <toy.raymond at gmail.com>
Date:   Tue Mar 29 09:58:45 2011 -0400

Make gamma accurate for integers; add precision test for
incomplete-gamma-tail.

qd-gamma.lisp:
o For integer values, just compute the gamma value directly by
multiplication.  This works around the problem that the current
algorithm is not as accurate as we would like.

rt-test.lisp:
o Reduce required accuracy in gamma-incomplete-tail.3.d.
o Add precision test for gamm incomplete tail near the negative real
axis.

diff --git a/qd-gamma.lisp b/qd-gamma.lisp
index 7048b80..8cab1bb 100644
--- a/qd-gamma.lisp
+++ b/qd-gamma.lisp
@@ -142,6 +142,18 @@
(/ (float-pi z)
(sin (* (float-pi z) z))
(gamma-aux (- 1 z) limit nterms)))
+	  ((and (zerop (imagpart z))
+		(= z (truncate z)))
+	   ;; We have gamma(n) where an integer value n and is small
+	   ;; enough.  In this case, just compute the product
+	   ;; directly.  We do this because our current implementation
+	   ;; has some round-off for these values, and that's annoying
+	   ;; and unexpected.
+	   (let ((n (truncate z)))
+	     (loop
+		for prod = (apply-contagion 1 precision) then (* prod k)
+		for k from 2 below n
+		finally (return (apply-contagion prod precision)))))
(t
(let ((absz (abs z)))
(cond ((>= absz limit)
@@ -168,7 +180,6 @@
(defmethod gamma ((z qd-complex))
(gamma-aux z 39 32))

-
;; Lentz's algorithm for evaluating continued fractions.
;;
;; Let the continued fraction be:
diff --git a/rt-tests.lisp b/rt-tests.lisp
index 66aab0c..acc8b74 100644
--- a/rt-tests.lisp
+++ b/rt-tests.lisp
@@ -1230,7 +1230,7 @@
(let* ((z -5d0)
(gi (incomplete-gamma-tail 2 z))
(true (* (+ z 1) (exp (- z)))))
-      (check-accuracy 52 gi true))
+      (check-accuracy 50 gi true))
nil)

(rt:deftest gamma-incomplete-tail.3.q
@@ -1240,6 +1240,16 @@
(check-accuracy 206 gi true))
nil)

+;; See http://www.wolframalpha.com/input/?i=Gamma[1%2F2%2C-100%2Bi%2F%2810^10%29]
+
+(rt:deftest gamma-incomplete-tail.4.q
+    (let* ((z #q(#q-100 #q1q-10))
+	   (gi (incomplete-gamma-tail 1/2 z))
+	   (true #q(#q-2.68811714181613544840818982228135651231579313476267430888499241497530341422025007816745898370049200133136q32
+		    #q-2.70176456134384383878883307528351227886457379834795655467745609829086928772079968479767583764284583465328q42)))
+      (check-accuracy 205 gi true))
+  nil)
+

;; Fresnel integrals.
;;

-----------------------------------------------------------------------

Summary of changes:
qd-gamma.lisp |   13 ++++++++++++-
rt-tests.lisp |   12 +++++++++++-
2 files changed, 23 insertions(+), 2 deletions(-)