Changeset 4b332ed2140e56c6fcaa689e6b8a48be36c988f6

Show
Ignore:
Timestamp:
03/21/12 21:53:37 (2 years 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:
2 modified

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  
    13641364  nil) 
    13651365 
     1366(rt:deftest psi.1d 
     1367    (let* ((z 1d0) 
     1368           (p (psi z)) 
     1369           (true (float (- +%gamma+) 1d0))) 
     1370      (check-accuracy 52 p true)) 
     1371  nil) 
     1372 
     1373(rt:deftest psi.1q 
     1374    (let* ((z #q1) 
     1375           (p (psi z)) 
     1376           (true (- +%gamma+))) 
     1377      (check-accuracy 208 p true)) 
     1378  nil) 
     1379 
     1380(rt:deftest psi.2d 
     1381    (let* ((z (float 4/3 1d0)) 
     1382           (p (psi z)) 
     1383           (true (- 3 
     1384                    +%gamma+ 
     1385                    (/ +pi+ (* 2 (sqrt #q3))) 
     1386                    (* 1.5 (log #q3))))) 
     1387      (check-accuracy 49.8 p true)) 
     1388  nil) 
     1389 
     1390(rt:deftest psi.2d 
     1391    (let* ((z (float 4/3 #q1)) 
     1392           (p (psi z)) 
     1393           (true (- 3 
     1394                    +%gamma+ 
     1395                    (/ +pi+ (* 2 (sqrt #q3))) 
     1396                    (* 1.5 (log #q3))))) 
     1397      (check-accuracy 205 p true)) 
     1398  nil) 
     1399 
     1400(rt:deftest psi.3d 
     1401    (let* ((z (float -1/2 1d0)) 
     1402           (p (psi z)) 
     1403           (true (- 2 
     1404                    +%gamma+ 
     1405                    (log #q4)))) 
     1406      (check-accuracy 48 p true)) 
     1407  nil) 
     1408 
     1409(rt:deftest psi.3q 
     1410    (let* ((z (float -1/2 #q1)) 
     1411           (p (psi z)) 
     1412           (true (- 2 
     1413                    +%gamma+ 
     1414                    (log #q4)))) 
     1415      (check-accuracy 204.1 p true)) 
     1416  nil) 
     1417 
    13661418(rt:deftest expintegral-e.1d 
    13671419    (let* ((z 1d0)