Changeset 7c5a3186070096ee93e16f2ddf51b2c84e7c5895 for qd-bessel.lisp

Timestamp:

04/11/12 09:18:31 (14 months ago)

Author:

Raymond Toy <rtoy@…>

Children:

c0e12ddf6b61f571555d46c6f168e6bebabc80b1, 25129063aaa6168763b447db6f033359a19b60e0

Parents:

bba9f8940c9f904bf14adc405d795a38ac333c24

git-committer:

Raymond Toy <rtoy@…> (04/11/12 09:18:31)

Message:

Correct some comments, remove unused code.

Files:

• qd-bessel.lisp (modified) (4 diffs)

qd-bessel.lisp

@@ -38 +38 @@
 ;; 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 +110 @@
 ;;
 ;; 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 +151 @@
 
 (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 +164 @@
           (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*