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

[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  7c5a3186070096ee93e16f2ddf51b2c84e7c5895 (commit)
      from  bba9f8940c9f904bf14adc405d795a38ac333c24 (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 7c5a3186070096ee93e16f2ddf51b2c84e7c5895
Author: Raymond Toy <rtoy at google.com>
Date:   Wed Apr 11 09:18:31 2012 -0700

    Correct some comments, remove unused code.

diff --git a/qd-bessel.lisp b/qd-bessel.lisp
index 367b847..0197e7f 100644
--- a/qd-bessel.lisp
+++ b/qd-bessel.lisp
@@ -37,7 +37,7 @@
 
 ;; B[k](p) = 1/2^(k+3/2)*integrate(exp(-p*u)*u^(k-1/2),u,0,1)
 ;;         = 1/2^(k+3/2)/p^(k+1/2)*integrate(t^(k-1/2)*exp(-t),t,0,p)
-;;         = 1/2^(k+3/2)/p^(k+1/2) * g(k+1/2, p)
+;;         = 1/2^(k+3/2)/p^(k+1/2) * G(k+1/2, p)
 ;;
 ;; where G(a,z) is the lower incomplete gamma function.
 ;;
@@ -109,6 +109,12 @@
 ;;                        = conj(B[k](%i*z).
 ;;
 ;; Hence I(-%i*z, v) = conj(I(%i*z, v)) when both z and v are real.
+;;
+;; Also note that when v is an integer of the form (2*m+1)/2, then
+;;   r[2*k+1](-2*%i*v) = r[2*k+1](-%i*(2*m+1))
+;;                     = -%i*(2*m+1)*product(-(2*m+1)^2+(2*j-1)^2, j, 1, k)
+;; so the product is zero when k >= m and the series I(p, q) is
+;; finite.
 (defun exp-arc-i (p q)
   (let* ((sqrt2 (sqrt (float 2 (realpart p))))
 	 (exp/p/sqrt2 (/ (exp (- p)) p sqrt2))
@@ -144,14 +150,9 @@
 	(format t " sum   - ~S~%" sum)))))
 
 (defun exp-arc-i-2 (p q)
-  (let* ((sqrt2 (sqrt (float 2 (realpart p))))
-	 (exp/p/sqrt2 (/ (exp (- p)) p sqrt2))
-	 (v (* #c(0 -2) q))
+  (let* ((v (* #c(0 -2) q))
 	 (v2 (expt v 2))
 	 (eps (epsilon (realpart p))))
-    (when *debug-exparc*
-      (format t "sqrt2 = ~S~%" sqrt2)
-      (format t "exp/p/sqrt2 = ~S~%" exp/p/sqrt2))
     (do* ((k 0 (1+ k))
 	  (bk (bk 0 p)
 	      (bk k p))
@@ -162,6 +163,12 @@
 		(* ratio bk))
 	  (sum term (+ sum term)))
 	 ((< (abs term) (* (abs sum) eps))
+	  (when *debug-exparc*
+	    (format t "Final k= ~D~%" k)
+	    (format t " bk    = ~S~%" bk)
+	    (format t " ratio = ~S~%" ratio)
+	    (format t " term  = ~S~%" term)
+	    (format t " sum   - ~S~%" sum))
 	  (* sum #c(0 2) (/ (exp p) q)))
       (when *debug-exparc*
 	(format t "k      = ~D~%" k)

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

Summary of changes:
 qd-bessel.lisp |   21 ++++++++++++++-------
 1 files changed, 14 insertions(+), 7 deletions(-)


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