This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "OCT: A portable Lisp implementation for quad-double precision floats".
The branch, master has been updated via c388f81713d7b2a483000d3cee1af030ed2c1cac (commit) from 26ef83a4543cdca6f752abca6f68e9d7f586c68e (commit)
Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below.
- Log ----------------------------------------------------------------- commit c388f81713d7b2a483000d3cee1af030ed2c1cac Author: Raymond Toy toy.raymond@gmail.com Date: Mon Dec 5 19:58:44 2011 -0800
Better exp-integral-e computation and fix for incomplete-gamma-tail.
For exp-integral-e, use the series for small z and the incomplete-gamma-tail for near the negative real axis. Otherwise, use the continued fraction.
In incomplete-gamma-tail, we were using the continued fraction instead of the incomplete-gamma function for the region just below the negative real axis. We should use the cf except in that region.
diff --git a/qd-gamma.lisp b/qd-gamma.lisp index f13b9a8..48f8345 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -396,7 +396,7 @@ ;; Use the incomplete gamma function to evaluate in this ;; region. (Arbitrarily selected the region to be a sector. ;; But what is the correct size of this sector?) - (if (<= (phase z) 3.1) + (if (<= (abs (phase z)) 3.1) (cf-incomplete-gamma-tail a z) (- (gamma a) (cf-incomplete-gamma a z)))))))
@@ -460,15 +460,7 @@ (/ (incomplete-gamma 1/2 (* z z)) (sqrt (float-pi z))))))
-(defun exp-integral-e (v z) - "Exponential integral E: - - E(v,z) = integrate(exp(-t)/t^v, t, 1, inf)" - ;; E(v,z) = z^(v-1) * integrate(t^(-v)*exp(-t), t, z, inf); - ;; - ;; for |arg(z)| < pi. - ;; - ;; +(defun cf-exp-integral-e (v z) ;; We use the continued fraction ;; ;; E(v,z) = exp(-z)/cf(z) @@ -485,7 +477,45 @@ (+ z+v (* 2 k))) #'(lambda (k) (* (- k) - (1- (+ k v)))))))) + (+ k v -1))))))) + +(defun s-exp-integral-e (v z) + ;; E(v,z) = gamma(1-v)*z^(v-1) - sum((-1)^k*z^k/(k-v+1)/k!, k, 0, inf) + (let ((-z (- z)) + (-v (- v)) + (eps (epsilon z))) + (loop for k from 0 + for term = 1 then (* term (/ -z k)) + for sum = (/ (- 1 v)) then (+ sum (/ term (+ k 1 -v))) + when (< (abs term) (* (abs sum) eps)) + return (- (* (gamma (- 1 v)) (expt z (- v 1))) + sum)))) + +(defun exp-integral-e (v z) + "Exponential integral E: + + E(v,z) = integrate(exp(-t)/t^v, t, 1, inf)" + ;; E(v,z) = z^(v-1) * integrate(t^(-v)*exp(-t), t, z, inf); + ;; + ;; for |arg(z)| < pi. + ;; + ;; + (cond ((< (abs z) 1) + ;; Use series for small z + (s-exp-integral-e v z)) + ((>= (abs (phase z)) 3.1) + ;; The continued fraction doesn't converge on the negative + ;; real axis, and converges very slowly near the negative + ;; real axis, so use the incomplete-gamma-tail function in + ;; this region. "Closeness" to the negative real axis is + ;; teken to mean that z is in a sector near the axis. + ;; + ;; E(v,z) = z^(v-1)*incomplete_gamma_tail(1-v,z) + (* (expt z (- v 1)) + (incomplete-gamma-tail (- 1 v) z))) + (t + ;; Use continued fraction for everything else. + (cf-exp-integral-e v z))))
;; Series for Fresnel S ;;
-----------------------------------------------------------------------
Summary of changes: qd-gamma.lisp | 52 +++++++++++++++++++++++++++++++++++++++++----------- 1 files changed, 41 insertions(+), 11 deletions(-)
hooks/post-receive