Changeset bc851c79e6fffe06c4383142bc2dfeb154fbde14 for qd-gamma.lisp

Timestamp:

03/23/12 08:38:53 (14 months 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:

• qd-gamma.lisp (modified) (2 diffs)

qd-gamma.lisp

@@ -498 +498 @@
              (= v (ftruncate v)))
         ;; v is an integer
-        (let ((n (truncate v)))
-          (- (* (/ (expt -z (- v 1))
+        (let* ((n (truncate v))
+               (n-1 (1- n)))
+          (- (* (/ (expt -z n-1)
                    (gamma v))
                 (- (psi v) (log z)))
-             (loop for k from 0 below n
+             (loop for k from 0
                    for term = 1 then (* term (/ -z k))
-                   for sum = (/ (- 1 v)) then (+ sum (/ term (+ k 1 -v)))
-                   when (< (abs term) (* (abs sum) eps))
-                     return sum)
-             (loop for k from n
-                   for term = 1 then (* term (/ -z k))
-                   for sum = 0 then (+ sum (/ term (+ k 1 -v)))
+                   for sum = (/ (- 1 v)) then (+ sum (let ((denom (- k n-1)))
+                                                       (if (zerop denom)
+                                                           0
+                                                           (/ term denom))))
                    when (< (abs term) (* (abs sum) eps))
                      return sum)))
@@ -528 +527 @@
   ;;
   ;;
-  (cond ((< (abs z) 1)
+  (cond ((and (realp v) (minusp v))
+         ;; E(-v, z) = z^(-v-1)*incomplete_gamma_tail(v+1,z)
+         (let ((-v (- v)))
+           (* (expt z (- v 1))
+              (incomplete-gamma-tail (+ -v 1) z))))
+        ((< (abs z) 1)
          ;; Use series for small z
          (s-exp-integral-e v z))