Changeset 4b332ed2140e56c6fcaa689e6b8a48be36c988f6

Timestamp:

03/21/12 21:53:37 (14 months ago)

Author:

Raymond Toy <toy.raymond@…>

Children:

6162b30bff1a78bcf7757476e8b05346af85f0e2, 1d9ec007bb4172eaccb8ca1db543c1218192cdb5

Parents:

9d3daf46c3e396941c8eb43209a45c0105217840

git-committer:

Raymond Toy <toy.raymond@…> (03/21/12 21:53:37)

Message:

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

Files:

• qd-gamma.lisp (modified) (2 diffs)
• rt-tests.lisp (modified) (1 diff)

qd-gamma.lisp

@@ -708 +708 @@
         (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
@@ -794 +796 @@
 
   (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)))))))

rt-tests.lisp

@@ -1364 +1364 @@
   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)