Changeset 97bc71a13186735474eba77043d685b0bd0a00d6
 Timestamp:
 03/18/11 06:06:10 (3 years ago)
 Author:
 Raymond Toy <toy.raymond@…>
 Children:
 6b8dc51d96004406780ec1faa2cd4ce3b127c287
 Parents:
 50fc6412a3effd9bea8f3e481f7c35c3c911f856
 gitcommitter:
 Raymond Toy <toy.raymond@…> (03/18/11 06:06:10)
 Message:

Add series for incompletegamma 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 INCOMPLETEGAMMATAIL, call INCOMPLETEGAMMA instead of
CFINCOMPLETEGAMMA just in case a and z are small.
o In INCOMPLETEGAMMA, use the series for small a and z.
o Simplify evaluation of Si(z) when z is real.
 Files:

Legend:
 Unmodified
 Added
 Removed

r50fc64

r97bc71


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 sincompletegamma (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 incompletegammatail (a z) 
… 
… 

277  292  ;; For real values, we split the result to compute either the 
278  293  ;; tail directly or from gamma(a)  incompletegamma 
279   (if (> z ( a 1)) 
 294  (if (> (abs z) (abs ( a 1))) 
280  295  (cfincompletegammatail a z) 
281   ( (gamma a) (cfincompletegamma a z))) 
 296  ( (gamma a) (incompletegamma a z))) 
282  297  (cfincompletegammatail a z)))) 
283  298  
… 
… 

289  304  (a (applycontagion a prec)) 
290  305  (z (applycontagion z prec))) 
291   (if (and (realp a) (realp z)) 
292   (if (< z ( a 1)) 
293   (cfincompletegamma a z) 
294   ( (gamma a) (cfincompletegammatail a z))) 
295   (if (< (abs z) (abs a)) 
296   (cfincompletegamma a z) 
297   ( (gamma a) (cfincompletegammatail a z)))))) 
 306  (if (and (< (abs a) 1) (< (abs z) 1)) 
 307  (sincompletegamma a z) 
 308  (if (and (realp a) (realp z)) 
 309  (if (< z ( a 1)) 
 310  (cfincompletegamma a z) 
 311  ( (gamma a) (cfincompletegammatail 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  (cfincompletegamma a z) 
 317  ( (gamma a) (cfincompletegammatail a z))))))) 
298  318  
299  319  (defun erf (z) 
… 
… 

406  426  (log iz))))))) 
407  427  (if (realp z) 
408   (if (< z 0) 
409   ( (sinintegral ( 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  ( (sinintegral ( z)))) 
 437  ((= z 0) 
 438  (* 0 z)) 
 439  (t 
 440  (+ (* 1/2 (floatpi z)) 
 441  (imagpart (incompletegammatail 0 (complex 0 z)))))) 
411  442  (si z)))) 
412  443  