[octscm] [octgit]OCT: A portable Lisp implementation for quaddouble precision floats branch master updated. c388f81713d7b2a483000d3cee1af030ed2c1cac
Raymond Toy
rtoy at commonlisp.net
Tue Dec 6 04:03:03 UTC 2011
This is an automated email from the git hooks/postreceive script. It was
generated because a ref change was pushed to the repository containing
the project "OCT: A portable Lisp implementation for quaddouble 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 at gmail.com>
Date: Mon Dec 5 19:58:44 2011 0800
Better expintegrale computation and fix for incompletegammatail.
For expintegrale, use the series for small z and the
incompletegammatail for near the negative real axis. Otherwise, use
the continued fraction.
In incompletegammatail, we were using the continued fraction instead
of the incompletegamma function for the region just below the negative
real axis. We should use the cf except in that region.
diff git a/qdgamma.lisp b/qdgamma.lisp
index f13b9a8..48f8345 100644
 a/qdgamma.lisp
+++ b/qdgamma.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)
(cfincompletegammatail a z)
( (gamma a) (cfincompletegamma a z)))))))
@@ 460,15 +460,7 @@
(/ (incompletegamma 1/2 (* z z))
(sqrt (floatpi z))))))
(defun expintegrale (v z)
 "Exponential integral E:

 E(v,z) = integrate(exp(t)/t^v, t, 1, inf)"
 ;; E(v,z) = z^(v1) * integrate(t^(v)*exp(t), t, z, inf);
 ;;
 ;; for arg(z) < pi.
 ;;
 ;;
+(defun cfexpintegrale (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 sexpintegrale (v z)
+ ;; E(v,z) = gamma(1v)*z^(v1)  sum((1)^k*z^k/(kv+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 expintegrale (v z)
+ "Exponential integral E:
+
+ E(v,z) = integrate(exp(t)/t^v, t, 1, inf)"
+ ;; E(v,z) = z^(v1) * integrate(t^(v)*exp(t), t, z, inf);
+ ;;
+ ;; for arg(z) < pi.
+ ;;
+ ;;
+ (cond ((< (abs z) 1)
+ ;; Use series for small z
+ (sexpintegrale 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 incompletegammatail 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^(v1)*incomplete_gamma_tail(1v,z)
+ (* (expt z ( v 1))
+ (incompletegammatail ( 1 v) z)))
+ (t
+ ;; Use continued fraction for everything else.
+ (cfexpintegrale v z))))
;; Series for Fresnel S
;;

Summary of changes:
qdgamma.lisp  52 +++++++++++++++++++++++++++++++++++++++++
1 files changed, 41 insertions(+), 11 deletions()
hooks/postreceive

OCT: A portable Lisp implementation for quaddouble precision floats
More information about the octcvs
mailing list