Changeset bc851c for qd-gamma.lisp


Ignore:
Timestamp:
03/23/12 15:38:53 (3 years ago)
Author:
Raymond Toy <rtoy@…>
Branches:
master
Children:
efbf11
Parents:
6162b30bff1a78bcf7757476e8b05346af85f0e2 (diff), fe8cff (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • qd-gamma.lisp

    r4b332e rfe8cff  
    498498             (= v (ftruncate v)))
    499499        ;; v is an integer
    500         (let ((n (truncate v)))
    501           (- (* (/ (expt -z (- v 1))
     500        (let* ((n (truncate v))
     501               (n-1 (1- n)))
     502          (- (* (/ (expt -z n-1)
    502503                   (gamma v))
    503504                (- (psi v) (log z)))
    504              (loop for k from 0 below n
     505             (loop for k from 0
    505506                   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)))
     507                   for sum = (/ (- 1 v)) then (+ sum (let ((denom (- k n-1)))
     508                                                       (if (zerop denom)
     509                                                           0
     510                                                           (/ term denom))))
    512511                   when (< (abs term) (* (abs sum) eps))
    513512                     return sum)))
     
    528527  ;;
    529528  ;;
    530   (cond ((< (abs z) 1)
     529  (cond ((and (realp v) (minusp v))
     530         ;; E(-v, z) = z^(-v-1)*incomplete_gamma_tail(v+1,z)
     531         (let ((-v (- v)))
     532           (* (expt z (- v 1))
     533              (incomplete-gamma-tail (+ -v 1) z))))
     534        ((< (abs z) 1)
    531535         ;; Use series for small z
    532536         (s-exp-integral-e v z))
Note: See TracChangeset for help on using the changeset viewer.