Show
Ignore:
Timestamp:
03/23/12 08:38:53 (2 years ago)
Author:
Raymond Toy <rtoy@…>
Children:
efbf11a6f5cd7e15f51cdb52fefd1584ab58c3d9
Parents:
6162b30bff1a78bcf7757476e8b05346af85f0e2, fe8cffb5ee8e7161addf7586ecaee00682c6bf1b
git-committer:
Raymond Toy <rtoy@…> (03/23/12 08:38:53)
Message:

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

Files:
1 modified

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))