[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 6b8a8a6fe864050f9c6371150d8070f9b38fe76e

[Raymond Toy]

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  6b8a8a6fe864050f9c6371150d8070f9b38fe76e (commit)
      from  b609574d330c7b0c9a6de588d3884a00f9aad13b (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 6b8a8a6fe864050f9c6371150d8070f9b38fe76e
Author: Raymond Toy <toy.raymond at gmail.com>
Date:   Sun Dec 4 23:22:01 2011 -0800

    Use exp-integral-e for (incomplete-gamma-tail 0 z).
    
    We already have exp-integral-e function so move expintegral-e
    implementation to exp-integral-e.

diff --git a/qd-gamma.lisp b/qd-gamma.lisp
index 05f85d4..8e523ed 100644
--- a/qd-gamma.lisp
+++ b/qd-gamma.lisp
@@ -375,27 +375,30 @@
   (let* ((prec (float-contagion a z))
 	 (a (apply-contagion a prec))
 	 (z (apply-contagion z prec)))
-    (if (and (zerop (imagpart a))
-	     (zerop (imagpart z)))
-	;; For real values, we split the result to compute either the
-	;; tail directly from the continued fraction or from gamma(a)
-	;; - incomplete-gamma.  The continued fraction doesn't
-	;; converge on the negative real axis, so we can't use that
-	;; there.  And accuracy appears to be better if z is "small".
-	;; We take this to mean |z| < |a-1|.  Note that |a-1| is the
-	;; peak of the integrand.
-	(if (and (> (abs z) (abs (- a 1)))
-		(not (minusp (realpart z))))
-	    (cf-incomplete-gamma-tail a z)
-	    (- (gamma a) (cf-incomplete-gamma a z)))
-	;; If the argument is close enough to the negative real axis,
-	;; the continued fraction for the tail is not very accurate.
-	;; 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)
-	    (cf-incomplete-gamma-tail a z)
-	    (- (gamma a) (cf-incomplete-gamma a z))))))
+    (if (zerop a)
+	;; incomplete_gamma_tail(0, z) = exp_integral_e(1,z)
+	(exp-integral-e 1 z)
+	(if (and (zerop (imagpart a))
+		 (zerop (imagpart z)))
+	    ;; For real values, we split the result to compute either the
+	    ;; tail directly from the continued fraction or from gamma(a)
+	    ;; - incomplete-gamma.  The continued fraction doesn't
+	    ;; converge on the negative real axis, so we can't use that
+	    ;; there.  And accuracy appears to be better if z is "small".
+	    ;; We take this to mean |z| < |a-1|.  Note that |a-1| is the
+	    ;; peak of the integrand.
+	    (if (and (> (abs z) (abs (- a 1)))
+		     (not (minusp (realpart z))))
+		(cf-incomplete-gamma-tail a z)
+		(- (gamma a) (cf-incomplete-gamma a z)))
+	    ;; If the argument is close enough to the negative real axis,
+	    ;; the continued fraction for the tail is not very accurate.
+	    ;; 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)
+		(cf-incomplete-gamma-tail a z)
+		(- (gamma a) (cf-incomplete-gamma a z)))))))
 
 (defun incomplete-gamma (a z)
   "Incomplete gamma function defined by:
@@ -461,9 +464,28 @@
   "Exponential integral E:
 
    E(v,z) = integrate(exp(-t)/t^v, t, 1, inf)"
-  ;; Wolfram gives E(v,z) = z^(v-1)*gamma_incomplete_tail(1-v,z)
-  (* (expt z (- v 1))
-     (incomplete-gamma-tail (- 1 v) z)))
+  ;; E(v,z) = z^(v-1) * integrate(t^(-v)*exp(-t), t, z, inf);
+  ;;
+  ;; for |arg(z)| < pi.
+  ;;
+  ;;
+  ;; We use the continued fraction
+  ;;
+  ;; E(v,z) = exp(-z)/cf(z)
+  ;;
+  ;; where the continued fraction cf(z) is
+  ;;
+  ;; a[k] = -k*(k+v-1)
+  ;; b[k] = v + 2*k + z
+  ;;
+  ;; for k = 1, inf
+  (let ((z+v (+ z v)))
+    (/ (exp (- z))
+       (lentz #'(lambda (k)
+		  (+ z+v (* 2 k)))
+	      #'(lambda (k)
+		  (* (- k)
+		     (1- (+ k v))))))))
 
 ;; Series for Fresnel S
 ;;
@@ -647,10 +669,3 @@
 ;; for k = 1, inf
 
 (defun expintegral-e (v z)
-  (let ((z+v (+ z v)))
-    (/ (exp (- z))
-       (lentz #'(lambda (k)
-		  (+ z+v (* 2 k)))
-	      #'(lambda (k)
-		  (* (- k)
-		     (1- (+ k v))))))))

-----------------------------------------------------------------------

Summary of changes:
 qd-gamma.lisp |   77 ++++++++++++++++++++++++++++++++++-----------------------
 1 files changed, 46 insertions(+), 31 deletions(-)


hooks/post-receive
-- 
OCT:  A portable Lisp implementation for quad-double precision floats