Changeset 0d5870201359817c679921a2d740fdd1697469b2

Show
Ignore:
Timestamp:
03/21/12 18:45:44 (2 years ago)
Author:
Raymond Toy <toy.raymond@…>
Children:
9d3daf46c3e396941c8eb43209a45c0105217840
Parents:
405df618a38d3b8ddaae691f865bbf068e931ac5
git-committer:
Raymond Toy <toy.raymond@…> (03/21/12 18:45:44)
Message:

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

Files:
1 modified

Legend:

Unmodified
Added
Removed
  • qd-gamma.lisp

    rc388f8 r0d5870  
    480480                     (+ k v -1))))))) 
    481481 
     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;; 
    482492(defun s-exp-integral-e (v z) 
    483493  ;; E(v,z) = gamma(1-v)*z^(v-1) - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf) 
     
    485495        (-v (- v)) 
    486496        (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))))) 
    493520 
    494521(defun exp-integral-e (v z) 
     
    680707        (realpart (ci z)) 
    681708        (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