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

[Raymond Toy]

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  af9fc6fcdb824757c5aeeadc6ccae2098d243313 (commit)
      from  08d254374188f97026d114451bba5d839dbdad11 (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 af9fc6fcdb824757c5aeeadc6ccae2098d243313
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(-)


hooks/post-receive
-- 
OCT:  A portable Lisp implementation for quad-double precision floats