Changeset 97bc71
 Timestamp:
 03/18/11 13:06:10 (4 years ago)
 Branches:
 master
 Children:
 6b8dc5
 Parents:
 50fc64
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

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