Commit 6162b30bff1a78bcf7757476e8b05346af85f0e2

Authored by Raymond Toy
2 parents 7eb8fdd7 4b332ed2

Merge branch 'master' of git://common-lisp.net/projects/oct/oct

Conflicts:
	qd-gamma.lisp
Showing 2 changed files with 98 additions and 3 deletions
qd-gamma.lisp
... ... @@ -504,7 +504,8 @@
504 504 (loop for k from 0 below n
505 505 for term = 1 then (* term (/ -z k))
506 506 for sum = (/ (- 1 v)) then (+ sum (/ term (+ k 1 -v)))
507   - finally (return sum))
  507 + when (< (abs term) (* (abs sum) eps))
  508 + return sum)
508 509 (loop for k from n
509 510 for term = 1 then (* term (/ -z k))
510 511 for sum = 0 then (+ sum (/ term (+ k 1 -v)))
... ... @@ -706,6 +707,8 @@
706 707 (realpart (ci z))
707 708 (ci z))))
708 709  
  710 +;; Array of values of the Bernoulli numbers. We only have enough for
  711 +;; the evaluation of the psi function.
709 712 (defconstant bern-values
710 713 (make-array 55
711 714 :initial-contents
... ... @@ -792,8 +795,15 @@
792 795 ;; formula to increase the argument and then apply the asymptotic formula.
793 796  
794 797 (cond ((minusp (realpart z))
795   - (- (psi (- 1 z))
796   - (* +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))))))
797 807 (t
798 808 (let* ((k (* 2 (1+ (floor (* .41 (- (log (epsilon z) 10)))))))
799 809 (m 0)
... ...
rt-tests.lisp
... ... @@ -1363,3 +1363,88 @@
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 +
  1418 +(rt:deftest expintegral-e.1d
  1419 + (let* ((z 1d0)
  1420 + (e (exp-integral-e 0 z))
  1421 + (true (/ (exp (- z)) z)))
  1422 + (check-accuracy 53 e true))
  1423 + nil)
  1424 +
  1425 +(rt:deftest expintegral-e.1q
  1426 + (let* ((z #q1)
  1427 + (e (exp-integral-e 0 z))
  1428 + (true (/ (exp (- z)) z)))
  1429 + (check-accuracy 212 e true))
  1430 + nil)
  1431 +
  1432 +(rt:deftest expintegral-e.2d
  1433 + (let* ((z 15d0)
  1434 + (e (exp-integral-e 0 z))
  1435 + (true (/ (exp (- z)) z)))
  1436 + (check-accuracy 53 e true))
  1437 + nil)
  1438 +
  1439 +(rt:deftest expintegral-e.2q
  1440 + (let* ((z #q15)
  1441 + (e (exp-integral-e 0 z))
  1442 + (true (/ (exp (- z)) z)))
  1443 + (check-accuracy 212 e true))
  1444 + nil)
  1445 +
  1446 +(rt:deftest expintegral-e.3d
  1447 + (let* ((e (exp-integral-e 2 1d0))
  1448 + (true 0.14849550677592204791835999d0))
  1449 + (check-accuracy 47.5 e true))
  1450 + nil)
1366 1451 \ No newline at end of file
... ...