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

[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  78801d6381aaaf4f21967582680f26889582db60 (commit)
      from  5566bc397937c2e8979ee6847e8f59f279f1f643 (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 78801d6381aaaf4f21967582680f26889582db60
Author: Raymond Toy <toy.raymond at gmail.com>
Date:   Sun Apr 15 10:50:16 2012 -0700

     * Add comments for {{{integer-bessel-j-exp-arc}}}.
     * Simplify {{{sum-an}}} so we stop the sum when the terms no longer
       contribute to the sum.
     * Change {{{big-n}}}.  This still needs work.

diff --git a/qd-bessel.lisp b/qd-bessel.lisp
index 0c06f8a..48a0dac 100644
--- a/qd-bessel.lisp
+++ b/qd-bessel.lisp
@@ -179,7 +179,13 @@
 	(format t " sum   - ~S~%" sum)))))
 
 
+;; Not really just for Bessel J for integer orders, but in that case,
+;; this is all that's needed to compute Bessel J.  For other values,
+;; this is just part of the computation needed.
 ;;
+;; Compute
+;;
+;;  1/(2*%pi) * (exp(-%i*v*%pi/2) * I(%i*z, v) + exp(%i*v*%pi/2) * I(-%i*z, v))
 (defun integer-bessel-j-exp-arc (v z)
   (let* ((iz (* #c(0 1) z))
 	 (i+ (exp-arc-i-2 iz v)))
@@ -257,6 +263,7 @@
 ;;
 ;; sum(exp(-k*z)*a[n](k,v), k, 1, N)
 ;;
+#+nil
 (defun sum-an (big-n n v z)
   (let ((sum 0))
     (loop for k from 1 upto big-n
@@ -265,6 +272,20 @@
 			  (an n k v))))
     sum))
 
+;; Like above, but we just stop when the terms no longer contribute to
+;; the sum.
+(defun sum-an (big-n n v z)
+  (let ((eps (epsilon (realpart z))))
+    (do* ((k 1 (+ 1 k))
+	  (term (* (exp (- (* k z)))
+		   (an n k v))
+		(* (exp (- (* k z)))
+		   (an n k v)))
+	  (sum term (+ sum term)))
+	 ((or (<= (abs term) (* eps (abs sum)))
+	      (> k big-n))
+	  sum))))
+
 ;; SUM-AB computes the series
 ;;
 ;; sum(alpha[n](z)*a[n](0,v) + beta[n](z)*sum_an(N, n, v, z), n, 0, inf)
@@ -399,7 +420,7 @@
 	   (integer-bessel-j-exp-arc v z))
 	  (t
 	   ;; Need to fine-tune the value of big-n.
-	   (let ((big-n 100)
+	   (let ((big-n 10)
 		 (vpi (* v (float-pi (realpart z)))))
 	     (+ (integer-bessel-j-exp-arc v z)
 		(if (= vv v)

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

Summary of changes:
 qd-bessel.lisp |   23 ++++++++++++++++++++++-
 1 files changed, 22 insertions(+), 1 deletions(-)


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