Commit 97bc71a13186735474eba77043d685b0bd0a00d6

Authored by Raymond Toy
1 parent 50fc6412

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.
Showing 1 changed file with 43 additions and 12 deletions
qd-gamma.lisp
... ... @@ -265,6 +265,21 @@
265 265 #'(lambda (n)
266 266 (- (* z (+ a n))))))))))
267 267  
  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 +
268 283 ;; Tail of the incomplete gamma function.
269 284 (defun incomplete-gamma-tail (a z)
270 285 "Tail of the incomplete gamma function defined by:
... ... @@ -276,9 +291,9 @@
276 291 (if (and (realp a) (realp z))
277 292 ;; For real values, we split the result to compute either the
278 293 ;; tail directly or from gamma(a) - incomplete-gamma
279   - (if (> z (- a 1))
  294 + (if (> (abs z) (abs (- a 1)))
280 295 (cf-incomplete-gamma-tail a z)
281   - (- (gamma a) (cf-incomplete-gamma a z)))
  296 + (- (gamma a) (incomplete-gamma a z)))
282 297 (cf-incomplete-gamma-tail a z))))
283 298  
284 299 (defun incomplete-gamma (a z)
... ... @@ -288,13 +303,18 @@
288 303 (let* ((prec (float-contagion a z))
289 304 (a (apply-contagion a prec))
290 305 (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)))))))
298 318  
299 319 (defun erf (z)
300 320 "Error function:
... ... @@ -405,9 +425,20 @@
405 425 (- (log -iz)
406 426 (log iz)))))))
407 427 (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))))))
411 442 (si z))))
412 443  
413 444 (defun cos-integral (z)
... ...