Changeset 97bc71a13186735474eba77043d685b0bd0a00d6

Show
Ignore:
Timestamp:
03/18/11 06:06:10 (3 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:
1 modified

Legend:

Unmodified
Added
Removed
  • qd-gamma.lisp

    r50fc64 r97bc71  
    266266                          (- (* z (+ a n)))))))))) 
    267267 
     268;; Series expansion for incomplete gamma.  Intended for |a|<1 and 
     269;; |z|<1.  The series is 
     270;; 
     271;; g(a,z) = z^a * sum((-z)^k/k!/(a+k), k, 0, inf) 
     272(defun s-incomplete-gamma (a z) 
     273  (let ((-z (- z)) 
     274        (eps (epsilon z))) 
     275    (loop for k from 0 
     276       for term = 1 then (* term (/ -z k)) 
     277       for sum = (/ a) then (+ sum (/ term (+ a k))) 
     278       when (< (abs term) (* (abs sum) eps)) 
     279       return (* sum (expt z a))))) 
     280 
     281   
     282 
    268283;; Tail of the incomplete gamma function. 
    269284(defun incomplete-gamma-tail (a z) 
     
    277292        ;; For real values, we split the result to compute either the 
    278293        ;; tail directly or from gamma(a) - incomplete-gamma 
    279         (if (> z (- a 1)) 
     294        (if (> (abs z) (abs (- a 1))) 
    280295            (cf-incomplete-gamma-tail a z) 
    281             (- (gamma a) (cf-incomplete-gamma a z))) 
     296            (- (gamma a) (incomplete-gamma a z))) 
    282297        (cf-incomplete-gamma-tail a z)))) 
    283298 
     
    289304         (a (apply-contagion a prec)) 
    290305         (z (apply-contagion z prec))) 
    291     (if (and (realp a) (realp z)) 
    292         (if (< z (- a 1)) 
    293             (cf-incomplete-gamma a z) 
    294             (- (gamma a) (cf-incomplete-gamma-tail a z))) 
    295         (if (< (abs z) (abs a)) 
    296             (cf-incomplete-gamma a z) 
    297             (- (gamma a) (cf-incomplete-gamma-tail a z)))))) 
     306    (if (and (< (abs a) 1) (< (abs z) 1)) 
     307        (s-incomplete-gamma a z) 
     308        (if (and (realp a) (realp z)) 
     309            (if (< z (- a 1)) 
     310                (cf-incomplete-gamma a z) 
     311                (- (gamma a) (cf-incomplete-gamma-tail a z))) 
     312            ;; The continued fraction doesn't converge very fast if a 
     313            ;; and z are small.  In this case, use the series 
     314            ;; expansion instead, which converges quite rapidly. 
     315            (if (< (abs z) (abs a)) 
     316                (cf-incomplete-gamma a z) 
     317                (- (gamma a) (cf-incomplete-gamma-tail a z))))))) 
    298318 
    299319(defun erf (z) 
     
    406426                      (log iz))))))) 
    407427    (if (realp z) 
    408         (if (< z 0) 
    409             (- (sin-integral (- z))) 
    410             (si z)) 
     428        ;; Si is odd and real for real z.  In this case, we have 
     429        ;; 
     430        ;; Si(x) = %i/2*(gamma_inc_tail(0, -%i*x) - gamma_inc_tail(0, %i*x) - %i*%pi) 
     431        ;;       = %pi/2 + %i/2*(gamma_inc_tail(0, -%i*x) - gamma_inc_tail(0, %i*x)) 
     432        ;; But gamma_inc_tail(0, conjugate(z)) = conjugate(gamma_inc_tail(0, z)), so 
     433        ;; 
     434        ;; Si(x) = %pi/2 + imagpart(gamma_inc_tail(0, %i*x)) 
     435        (cond ((< z 0) 
     436               (- (sin-integral (- z)))) 
     437              ((= z 0) 
     438               (* 0 z)) 
     439              (t 
     440               (+ (* 1/2 (float-pi z)) 
     441                  (imagpart (incomplete-gamma-tail 0 (complex 0 z)))))) 
    411442        (si z)))) 
    412443