Commit 0d5870201359817c679921a2d740fdd1697469b2

Authored by Raymond Toy
1 parent 405df618

Implement psi and fix exp-integral-e for integral values of v. Needs

some more work.
Showing 1 changed file with 139 additions and 6 deletions
qd-gamma.lisp
... ... @@ -479,17 +479,44 @@
479 479 (* (- k)
480 480 (+ k v -1)))))))
481 481  
  482 +
  483 +;; For v not an integer:
  484 +;;
  485 +;; E(v,z) = gamma(1-v)*z^(v-1) - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf)
  486 +;;
  487 +;; For v an integer:
  488 +;;
  489 +;; E(v,z) = (-z)^(v-1)/(v-1)!*(psi(v)-log(z))
  490 +;; - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf, k != n-1)
  491 +;;
482 492 (defun s-exp-integral-e (v z)
483 493 ;; E(v,z) = gamma(1-v)*z^(v-1) - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf)
484 494 (let ((-z (- z))
485 495 (-v (- v))
486 496 (eps (epsilon z)))
487   - (loop for k from 0
488   - for term = 1 then (* term (/ -z k))
489   - for sum = (/ (- 1 v)) then (+ sum (/ term (+ k 1 -v)))
490   - when (< (abs term) (* (abs sum) eps))
491   - return (- (* (gamma (- 1 v)) (expt z (- v 1)))
492   - sum))))
  497 + (if (and (realp v)
  498 + (= v (ftruncate v)))
  499 + ;; v is an integer
  500 + (let ((n (truncate v)))
  501 + (- (* (/ (expt -z (- v 1))
  502 + (gamma v))
  503 + (- (psi v) (log z)))
  504 + (loop for k from 0 below n
  505 + for term = 1 then (* term (/ -z k))
  506 + for sum = (/ (- 1 v)) then (+ sum (/ term (+ k 1 -v)))
  507 + when (< (abs term) (* (abs sum) eps))
  508 + return sum)
  509 + (loop for k from n
  510 + for term = 1 then (* term (/ -z k))
  511 + for sum = 0 then (+ sum (/ term (+ k 1 -v)))
  512 + when (< (abs term) (* (abs sum) eps))
  513 + return sum)))
  514 + (loop for k from 0
  515 + for term = 1 then (* term (/ -z k))
  516 + for sum = (/ (- 1 v)) then (+ sum (/ term (+ k 1 -v)))
  517 + when (< (abs term) (* (abs sum) eps))
  518 + return (- (* (gamma (- 1 v)) (expt z (- v 1)))
  519 + sum)))))
493 520  
494 521 (defun exp-integral-e (v z)
495 522 "Exponential integral E:
... ... @@ -679,3 +706,109 @@
679 706 (if (and (realp z) (plusp z))
680 707 (realpart (ci z))
681 708 (ci z))))
  709 +
  710 +(defconstant bern-values
  711 + (make-array 55
  712 + :initial-contents
  713 + '(1
  714 + -1/2
  715 + 1/6
  716 + 0
  717 + -1/30
  718 + 0
  719 + 1/42
  720 + 0
  721 + -1/30
  722 + 0
  723 + 5/66
  724 + 0
  725 + -691/2730
  726 + 0
  727 + 7/6
  728 + 0
  729 + -3617/510
  730 + 0
  731 + 43867/798
  732 + 0
  733 + -174611/330
  734 + 0
  735 + 854513/138
  736 + 0
  737 + -236364091/2730
  738 + 0
  739 + 8553103/6
  740 + 0
  741 + -23749461029/870
  742 + 0
  743 + 8615841276005/14322
  744 + 0
  745 + -7709321041217/510
  746 + 0
  747 + 2577687858367/6
  748 + 0
  749 + -26315271553053477373/1919190
  750 + 0
  751 + 2929993913841559/6
  752 + 0
  753 + -261082718496449122051/13530
  754 + 0
  755 + 1520097643918070802691/1806
  756 + 0
  757 + -27833269579301024235023/690
  758 + 0
  759 + 596451111593912163277961/282
  760 + 0
  761 + -5609403368997817686249127547/46410
  762 + 0
  763 + 495057205241079648212477525/66
  764 + 0
  765 + -801165718135489957347924991853/1590
  766 + 0
  767 + 29149963634884862421418123812691/798
  768 + )))
  769 +
  770 +(defun bern (k)
  771 + (aref bern-values k))
  772 +
  773 +(defun psi (z)
  774 + "Digamma function defined by
  775 +
  776 + - %gamma + sum(1/k-1/(k+z-1), k, 1, inf)
  777 +
  778 + where %gamma is Euler's constant"
  779 +
  780 + ;; A&S 6.3.7: Reflection formula
  781 + ;;
  782 + ;; psi(1-z) = psi(z) + %pi*cot(%pi*z)
  783 + ;;
  784 + ;; A&S 6.3.6: Recurrence formula
  785 + ;;
  786 + ;; psi(n+z) = 1/(z+n-1)+1/(z+n-2)+...+1/(z+2)+1/(1+z)+psi(1+z)
  787 + ;;
  788 + ;; A&S 6.3.8: Asymptotic formula
  789 + ;;
  790 + ;; psi(z) ~ log(z) - sum(bern(2*n)/(2*n*z^(2*n)), n, 1, inf)
  791 + ;;
  792 + ;; So use reflection formula if Re(z) < 0. For z > 0, use the recurrence
  793 + ;; formula to increase the argument and then apply the asymptotic formula.
  794 +
  795 + (cond ((minusp (realpart z))
  796 + (- (psi (- 1 z))
  797 + (* +pi+ (/ (tan (* +pi+ z))))))
  798 + (t
  799 + (let* ((k (* 2 (1+ (floor (* .41 (- (log (epsilon z) 10)))))))
  800 + (m 0)
  801 + (y (expt (+ z k) 2))
  802 + (x 0))
  803 + (loop for i from 1 upto (floor k 2) do
  804 + (progn
  805 + (incf m (+ (/ (+ z i i -1))
  806 + (/ (+ z i i -2))))
  807 + (setf x (/ (+ x (/ (bern (+ k 2 (* -2 i)))
  808 + (- k i i -2)))
  809 + y))))
  810 + (- (log (+ z k))
  811 + (/ (* 2 (+ z k)))
  812 + x
  813 + m)))))
  814 +
... ...