Index: qd-bessel.lisp
===================================================================
--- qd-bessel.lisp	(revision e9cd1a46fcf2a5c0101b3473648cb242b55987e8)
+++ qd-bessel.lisp	(revision 7c5a3186070096ee93e16f2ddf51b2c84e7c5895)
@@ -38,5 +38,5 @@
 ;; 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.
@@ -110,4 +110,10 @@
 ;;
 ;; 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))))
@@ -145,12 +151,7 @@
 
 (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)
@@ -163,4 +164,10 @@
 	  (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*