Changeset 97bc71a13186735474eba77043d685b0bd0a00d6

Timestamp:

03/18/11 06:06:10 (2 years ago)

Author:

Raymond Toy <toy.raymond@…>

Children:

6b8dc51d96004406780ec1faa2cd4ce3b127c287

Parents:

50fc6412a3effd9bea8f3e481f7c35c3c911f856

git-committer:

Raymond Toy <toy.raymond@…> (03/18/11 06:06:10)

Message:

Add series for incomplete-gamma for when the fraction is slow.

o Add series for incomplete gamma function for small a and z. Needed

because the continued fraction is slow in this range.

o In INCOMPLETE-GAMMA-TAIL, call INCOMPLETE-GAMMA instead of

CF-INCOMPLETE-GAMMA just in case a and z are small.

o In INCOMPLETE-GAMMA, use the series for small a and z.
o Simplify evaluation of Si(z) when z is real.

Files:

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

qd-gamma.lisp

@@ -266 +266 @@
                           (- (* z (+ a n))))))))))
 
+;; Series expansion for incomplete gamma.  Intended for |a|<1 and
+;; |z|<1.  The series is
+;;
+;; g(a,z) = z^a * sum((-z)^k/k!/(a+k), k, 0, inf)
+(defun s-incomplete-gamma (a z)
+  (let ((-z (- z))
+        (eps (epsilon z)))
+    (loop for k from 0
+       for term = 1 then (* term (/ -z k))
+       for sum = (/ a) then (+ sum (/ term (+ a k)))
+       when (< (abs term) (* (abs sum) eps))
+       return (* sum (expt z a)))))
+
+  
+
 ;; Tail of the incomplete gamma function.
 (defun incomplete-gamma-tail (a z)
@@ -277 +292 @@
         ;; For real values, we split the result to compute either the
         ;; tail directly or from gamma(a) - incomplete-gamma
-        (if (> z (- a 1))
+        (if (> (abs z) (abs (- a 1)))
             (cf-incomplete-gamma-tail a z)
-            (- (gamma a) (cf-incomplete-gamma a z)))
+            (- (gamma a) (incomplete-gamma a z)))
         (cf-incomplete-gamma-tail a z))))
 
@@ -289 +304 @@
          (a (apply-contagion a prec))
          (z (apply-contagion z prec)))
-    (if (and (realp a) (realp z))
-        (if (< z (- a 1))
-            (cf-incomplete-gamma a z)
-            (- (gamma a) (cf-incomplete-gamma-tail a z)))
-        (if (< (abs z) (abs a))
-            (cf-incomplete-gamma a z)
-            (- (gamma a) (cf-incomplete-gamma-tail a z))))))
+    (if (and (< (abs a) 1) (< (abs z) 1))
+        (s-incomplete-gamma a z)
+        (if (and (realp a) (realp z))
+            (if (< z (- a 1))
+                (cf-incomplete-gamma a z)
+                (- (gamma a) (cf-incomplete-gamma-tail a z)))
+            ;; The continued fraction doesn't converge very fast if a
+            ;; and z are small.  In this case, use the series
+            ;; expansion instead, which converges quite rapidly.
+            (if (< (abs z) (abs a))
+                (cf-incomplete-gamma a z)
+                (- (gamma a) (cf-incomplete-gamma-tail a z)))))))
 
 (defun erf (z)
@@ -406 +426 @@
                       (log iz)))))))
     (if (realp z)
-        (if (< z 0)
-            (- (sin-integral (- z)))
-            (si z))
+        ;; Si is odd and real for real z.  In this case, we have
+        ;;
+        ;; Si(x) = %i/2*(gamma_inc_tail(0, -%i*x) - gamma_inc_tail(0, %i*x) - %i*%pi)
+        ;;       = %pi/2 + %i/2*(gamma_inc_tail(0, -%i*x) - gamma_inc_tail(0, %i*x))
+        ;; But gamma_inc_tail(0, conjugate(z)) = conjugate(gamma_inc_tail(0, z)), so
+        ;;
+        ;; Si(x) = %pi/2 + imagpart(gamma_inc_tail(0, %i*x))
+        (cond ((< z 0)
+               (- (sin-integral (- z))))
+              ((= z 0)
+               (* 0 z))
+              (t
+               (+ (* 1/2 (float-pi z))
+                  (imagpart (incomplete-gamma-tail 0 (complex 0 z))))))
         (si z))))