Changeset 97bc71


Ignore:
Timestamp:
03/18/11 13:06:10 (4 years ago)
Author:
Raymond Toy <toy.raymond@…>
Branches:
master
Children:
6b8dc5
Parents:
50fc64
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.

File:
1 edited

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
Note: See TracChangeset for help on using the changeset viewer.