Changeset 0d5870


Ignore:
Timestamp:
03/22/12 01:45:44 (3 years ago)
Author:
Raymond Toy <toy.raymond@…>
Branches:
master
Children:
9d3daf
Parents:
405df6
Message:

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

File:
1 edited

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
Note: See TracChangeset for help on using the changeset viewer.