Commit cb1a5d41baf9d4db12a2563a230d0a1e55c8adea

Authored by Raymond Toy
1 parent 6cfb0ac4

Fix some mistakes in exp-integral-e when changing algorithm to use

series.
Showing 1 changed file with 9 additions and 24 deletions
qd-gamma.lisp
... ... @@ -532,36 +532,21 @@
532 532 ;; for |arg(z)| < pi.
533 533 ;;
534 534 ;;
535   - (let* ((prec (float-contagion v z))
536   - (v (apply-contagion v prec))
537   - (z (apply-contagion z prec)))
  535 + (with-floating-point-contagion (v z)
538 536 (cond ((and (realp v) (minusp v))
539 537 ;; E(-v, z) = z^(-v-1)*incomplete_gamma_tail(v+1,z)
540 538 (let ((-v (- v)))
541 539 (* (expt z (- v 1))
542 540 (incomplete-gamma-tail (+ -v 1) z))))
543   - ((< (abs z) 1)
544   - ;; Use series for small z
  541 + ((or (< (abs z) 1) (>= (abs (phase z)) 3.1))
  542 + ;; Use series for small z or if z is near the negative real
  543 + ;; axis because the continued fraction does not converge on
  544 + ;; the negative axis and converges slowly near the negative
  545 + ;; axis.
545 546 (s-exp-integral-e v z))
546   - ((>= (abs (phase z)) 3.1)
547   - ;; The continued fraction doesn't converge on the negative
548   - ;; real axis, and converges very slowly near the negative
549   - ;; real axis, so use the incomplete-gamma-tail function in
550   - ;; this region. "Closeness" to the negative real axis is
551   - ;; teken to mean that z is in a sector near the axis.
552   - ;;
553   - ;; E(v,z) = z^(v-1)*incomplete_gamma_tail(1-v,z)
554   - (* (expt z (- v 1))
555   - (incomplete-gamma-tail (+ -v 1) z))))
556   - ((or (< (abs z) 1) (>= (abs (phase z)) 3.1))
557   - ;; Use series for small z or if z is near the negative real
558   - ;; axis because the continued fraction does not converge on
559   - ;; the negative axis and converges slowly near the negative
560   - ;; axis.
561   - (s-exp-integral-e v z))
562   - (t
563   - ;; Use continued fraction for everything else.
564   - (cf-exp-integral-e v z))))
  547 + (t
  548 + ;; Use continued fraction for everything else.
  549 + (cf-exp-integral-e v z)))))
565 550  
566 551 ;; Series for Fresnel S
567 552 ;;
... ...