Commit f77890744f8e35f9c38ba19d435c6541f92d020b

Authored by Raymond Toy
2 parents 06e250ea c388f817

Merge branch 'master' of git://common-lisp.net/projects/oct/oct

Showing 1 changed file with 41 additions and 11 deletions
qd-gamma.lisp
... ... @@ -396,7 +396,7 @@
396 396 ;; Use the incomplete gamma function to evaluate in this
397 397 ;; region. (Arbitrarily selected the region to be a sector.
398 398 ;; But what is the correct size of this sector?)
399   - (if (<= (phase z) 3.1)
  399 + (if (<= (abs (phase z)) 3.1)
400 400 (cf-incomplete-gamma-tail a z)
401 401 (- (gamma a) (cf-incomplete-gamma a z)))))))
402 402  
... ... @@ -460,15 +460,7 @@
460 460 (/ (incomplete-gamma 1/2 (* z z))
461 461 (sqrt (float-pi z))))))
462 462  
463   -(defun exp-integral-e (v z)
464   - "Exponential integral E:
465   -
466   - E(v,z) = integrate(exp(-t)/t^v, t, 1, inf)"
467   - ;; E(v,z) = z^(v-1) * integrate(t^(-v)*exp(-t), t, z, inf);
468   - ;;
469   - ;; for |arg(z)| < pi.
470   - ;;
471   - ;;
  463 +(defun cf-exp-integral-e (v z)
472 464 ;; We use the continued fraction
473 465 ;;
474 466 ;; E(v,z) = exp(-z)/cf(z)
... ... @@ -485,7 +477,45 @@
485 477 (+ z+v (* 2 k)))
486 478 #'(lambda (k)
487 479 (* (- k)
488   - (1- (+ k v))))))))
  480 + (+ k v -1)))))))
  481 +
  482 +(defun s-exp-integral-e (v z)
  483 + ;; E(v,z) = gamma(1-v)*z^(v-1) - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf)
  484 + (let ((-z (- z))
  485 + (-v (- v))
  486 + (eps (epsilon z)))
  487 + (loop for k from 0
  488 + for term = 1 then (* term (/ -z k))
  489 + for sum = (/ (- 1 v)) then (+ sum (/ term (+ k 1 -v)))
  490 + when (< (abs term) (* (abs sum) eps))
  491 + return (- (* (gamma (- 1 v)) (expt z (- v 1)))
  492 + sum))))
  493 +
  494 +(defun exp-integral-e (v z)
  495 + "Exponential integral E:
  496 +
  497 + E(v,z) = integrate(exp(-t)/t^v, t, 1, inf)"
  498 + ;; E(v,z) = z^(v-1) * integrate(t^(-v)*exp(-t), t, z, inf);
  499 + ;;
  500 + ;; for |arg(z)| < pi.
  501 + ;;
  502 + ;;
  503 + (cond ((< (abs z) 1)
  504 + ;; Use series for small z
  505 + (s-exp-integral-e v z))
  506 + ((>= (abs (phase z)) 3.1)
  507 + ;; The continued fraction doesn't converge on the negative
  508 + ;; real axis, and converges very slowly near the negative
  509 + ;; real axis, so use the incomplete-gamma-tail function in
  510 + ;; this region. "Closeness" to the negative real axis is
  511 + ;; teken to mean that z is in a sector near the axis.
  512 + ;;
  513 + ;; E(v,z) = z^(v-1)*incomplete_gamma_tail(1-v,z)
  514 + (* (expt z (- v 1))
  515 + (incomplete-gamma-tail (- 1 v) z)))
  516 + (t
  517 + ;; Use continued fraction for everything else.
  518 + (cf-exp-integral-e v z))))
489 519  
490 520 ;; Series for Fresnel S
491 521 ;;
... ...