Changeset 4b332e


Ignore:
Timestamp:
03/22/12 04:53:37 (3 years ago)
Author:
Raymond Toy <toy.raymond@…>
Branches:
master
Children:
1d9ec0, b1d9be
Parents:
9d3daf
Message:

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

Files:
2 edited

Legend:

Unmodified
Added
Removed
  • qd-gamma.lisp

    r0d5870 r4b332e  
    708708        (ci z))))
    709709
     710;; Array of values of the Bernoulli numbers.  We only have enough for
     711;; the evaluation of the psi function.
    710712(defconstant bern-values
    711713  (make-array 55
     
    794796
    795797  (cond ((minusp (realpart z))
    796          (- (psi (- 1 z))
    797             (* +pi+ (/ (tan (* +pi+ z))))))
     798         (let ((p (float +pi+ (realpart z))))
     799           (flet ((cot-pi (z)
     800                    ;; cot(%pi*z), car
     801                    (handler-case
     802                        (/ (tan (* p z)))
     803                      (division-by-zero ()
     804                        (* 0 z)))))
     805             (- (psi (- 1 z))
     806                (* p (cot-pi z))))))
    798807        (t
    799808         (let* ((k (* 2 (1+ (floor (* .41 (- (log (epsilon z) 10)))))))
  • rt-tests.lisp

    r9d3daf r4b332e  
    13691369  nil)
    13701370
     1371(rt:deftest psi.1d
     1372    (let* ((z 1d0)
     1373           (p (psi z))
     1374           (true (float (- +%gamma+) 1d0)))
     1375      (check-accuracy 52 p true))
     1376  nil)
     1377
     1378(rt:deftest psi.1q
     1379    (let* ((z #q1)
     1380           (p (psi z))
     1381           (true (- +%gamma+)))
     1382      (check-accuracy 208 p true))
     1383  nil)
     1384
     1385(rt:deftest psi.2d
     1386    (let* ((z (float 4/3 1d0))
     1387           (p (psi z))
     1388           (true (- 3
     1389                    +%gamma+
     1390                    (/ +pi+ (* 2 (sqrt #q3)))
     1391                    (* 1.5 (log #q3)))))
     1392      (check-accuracy 49.8 p true))
     1393  nil)
     1394
     1395(rt:deftest psi.2d
     1396    (let* ((z (float 4/3 #q1))
     1397           (p (psi z))
     1398           (true (- 3
     1399                    +%gamma+
     1400                    (/ +pi+ (* 2 (sqrt #q3)))
     1401                    (* 1.5 (log #q3)))))
     1402      (check-accuracy 205 p true))
     1403  nil)
     1404
     1405(rt:deftest psi.3d
     1406    (let* ((z (float -1/2 1d0))
     1407           (p (psi z))
     1408           (true (- 2
     1409                    +%gamma+
     1410                    (log #q4))))
     1411      (check-accuracy 48 p true))
     1412  nil)
     1413
     1414(rt:deftest psi.3q
     1415    (let* ((z (float -1/2 #q1))
     1416           (p (psi z))
     1417           (true (- 2
     1418                    +%gamma+
     1419                    (log #q4))))
     1420      (check-accuracy 204.1 p true))
     1421  nil)
     1422
    13711423(rt:deftest expintegral-e.1d
    13721424    (let* ((z 1d0)
Note: See TracChangeset for help on using the changeset viewer.