Commit 4b332ed2140e56c6fcaa689e6b8a48be36c988f6

Authored by Raymond Toy
1 parent 9d3daf46

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

tests too.
Showing 2 changed files with 63 additions and 2 deletions
qd-gamma.lisp
... ... @@ -707,6 +707,8 @@
707 707 (realpart (ci z))
708 708 (ci z))))
709 709  
  710 +;; Array of values of the Bernoulli numbers. We only have enough for
  711 +;; the evaluation of the psi function.
710 712 (defconstant bern-values
711 713 (make-array 55
712 714 :initial-contents
... ... @@ -793,8 +795,15 @@
793 795 ;; formula to increase the argument and then apply the asymptotic formula.
794 796  
795 797 (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))))))
798 807 (t
799 808 (let* ((k (* 2 (1+ (floor (* .41 (- (log (epsilon z) 10)))))))
800 809 (m 0)
... ...
rt-tests.lisp
... ... @@ -1363,6 +1363,58 @@
1363 1363 (check-accuracy 212 s true))
1364 1364 nil)
1365 1365  
  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 +
1366 1418 (rt:deftest expintegral-e.1d
1367 1419 (let* ((z 1d0)
1368 1420 (e (exp-integral-e 0 z))
... ...