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

[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  4b332ed2140e56c6fcaa689e6b8a48be36c988f6 (commit)
      from  9d3daf46c3e396941c8eb43209a45c0105217840 (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 4b332ed2140e56c6fcaa689e6b8a48be36c988f6
Author: Raymond Toy <toy.raymond at gmail.com>
Date:   Wed Mar 21 21:53:37 2012 -0700

    Fix bug in psi for -n/2 for n odd which was causing an overflow.  Add
    tests too.

diff --git a/qd-gamma.lisp b/qd-gamma.lisp
index 02b5547..d49acf8 100644
--- a/qd-gamma.lisp
+++ b/qd-gamma.lisp
@@ -707,6 +707,8 @@
 	(realpart (ci z))
 	(ci z))))
 
+;; Array of values of the Bernoulli numbers.  We only have enough for
+;; the evaluation of the psi function.
 (defconstant bern-values
   (make-array 55
 	      :initial-contents
@@ -793,8 +795,15 @@
   ;; formula to increase the argument and then apply the asymptotic formula.
 
   (cond ((minusp (realpart z))
-	 (- (psi (- 1 z))
-	    (* +pi+ (/ (tan (* +pi+ z))))))
+	 (let ((p (float +pi+ (realpart z))))
+	   (flet ((cot-pi (z)
+		    ;; cot(%pi*z), car
+		    (handler-case
+			(/ (tan (* p z)))
+		      (division-by-zero ()
+		        (* 0 z)))))
+	     (- (psi (- 1 z))
+		(* p (cot-pi z))))))
 	(t
 	 (let* ((k (* 2 (1+ (floor (* .41 (- (log (epsilon z) 10)))))))
 		(m 0)
diff --git a/rt-tests.lisp b/rt-tests.lisp
index 9426c98..c5b6743 100644
--- a/rt-tests.lisp
+++ b/rt-tests.lisp
@@ -1363,6 +1363,58 @@
       (check-accuracy 212 s true))
   nil)
 
+(rt:deftest psi.1d
+    (let* ((z 1d0)
+	   (p (psi z))
+	   (true (float (- +%gamma+) 1d0)))
+      (check-accuracy 52 p true))
+  nil)
+
+(rt:deftest psi.1q
+    (let* ((z #q1)
+	   (p (psi z))
+	   (true (- +%gamma+)))
+      (check-accuracy 208 p true))
+  nil)
+
+(rt:deftest psi.2d
+    (let* ((z (float 4/3 1d0))
+	   (p (psi z))
+	   (true (- 3
+		    +%gamma+
+		    (/ +pi+ (* 2 (sqrt #q3)))
+		    (* 1.5 (log #q3)))))
+      (check-accuracy 49.8 p true))
+  nil)
+
+(rt:deftest psi.2d
+    (let* ((z (float 4/3 #q1))
+	   (p (psi z))
+	   (true (- 3
+		    +%gamma+
+		    (/ +pi+ (* 2 (sqrt #q3)))
+		    (* 1.5 (log #q3)))))
+      (check-accuracy 205 p true))
+  nil)
+
+(rt:deftest psi.3d
+    (let* ((z (float -1/2 1d0))
+	   (p (psi z))
+	   (true (- 2
+		    +%gamma+
+		    (log #q4))))
+      (check-accuracy 48 p true))
+  nil)
+
+(rt:deftest psi.3q
+    (let* ((z (float -1/2 #q1))
+	   (p (psi z))
+	   (true (- 2
+		    +%gamma+
+		    (log #q4))))
+      (check-accuracy 204.1 p true))
+  nil)
+
 (rt:deftest expintegral-e.1d
     (let* ((z 1d0)
 	   (e (exp-integral-e 0 z))

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

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


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