Changeset a5a4c7acd95d1946e7cf426f9a927a918c9e2afc

Timestamp:

04/07/12 09:28:16 (14 months ago)

Author:

Raymond Toy <toy.raymond@…>

Children:

015558a3df8fae2686c7b6e2f049da3f4b1a2cd8

Parents:

1d404aca1e6652ad2456ba14e8bbdf5d329c06a0

git-committer:

Raymond Toy <toy.raymond@…> (04/07/12 09:28:16)

Message:

Add more parts of the exp-arc algorithm. Needs lots of work, but it
seems that bessel_j(n,z) mostly works.

Files:

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

qd-bessel.lisp

@@ -40 +40 @@
 ;;         = 1/2^(k+3/2)/p^(k+1/2) * g(k+1/2, p)
 ;;
-;; where g(a,z) is the lower incomplete gamma function.
-;;
-;; There is the continued fraction expansion for g(a,z) (see
+;; where G(a,z) is the lower incomplete gamma function.
+;;
+;; There is the continued fraction expansion for G(a,z) (see
 ;; cf-incomplete-gamma in qd-gamma.lisp):
 ;;
-;;  g(a,z) = z^a*exp(-z)/ CF
+;;  G(a,z) = z^a*exp(-z)/ CF
 ;;
 ;; So
@@ -184 +184 @@
        (float-pi i+)
        2)))
+
+;; alpha[n](z) = integrate(exp(-z*s)*s^n, s, 0, 1/2)
+;; beta[n](z)  = integrate(exp(-z*s)*s^n, s, -1/2, 1/2)
+;;
+;; The recurrence in [2] is
+;;
+;; alpha[n](z) = - exp(-z/2)/2^n/z + n/z*alpha[n-1](z)
+;; beta[n]z)   = ((-1)^n*exp(z/2)-exp(-z/2))/2^n/z + n/z*beta[n-1](z)
+;;
+;; We also note that
+;;
+;; alpha[n](z) = G(n+1,z/2)/z^(n+1)
+;; beta[n](z)  = G(n+1,z/2)/z^(n+1) - G(n+1,-z/2)/z^(n+1)
+
+(defun alpha (n z)
+  (let ((n (float n (realpart z))))
+    (/ (cf-incomplete-gamma (1+ n) (/ z 2))
+       (expt z (1+ n)))))
+
+(defun beta (n z)
+  (let ((n (float n (realpart z))))
+    (/ (- (cf-incomplete-gamma (1+ n) (/ z 2))
+          (cf-incomplete-gamma (1+ n) (/ z -2)))
+       (expt z (1+ n)))))
+
+;; a[0](k,v) := (k+sqrt(k^2+1))^(-v);
+;; a[1](k,v) := -v*a[0](k,v)/sqrt(k^2+1);
+;; a[n](k,v) := 1/(k^2+1)/(n-1)/n*((v^2-(n-2)^2)*a[n-2](k,v)-k*(n-1)*(2*n-3)*a[n-1](k,v));
+
+;; Convert this to iteration instead of using this quick-and-dirty
+;; memoization?
+(let ((hash (make-hash-table :test 'equal)))
+  (defun an-clrhash ()
+    (clrhash hash))
+  (defun an-dump-hash ()
+    (maphash #'(lambda (k v)
+                 (format t "~S -> ~S~%" k v))
+             hash))
+  (defun an (n k v)
+    (or (gethash (list n k v) hash)
+        (let ((result
+                (cond ((= n 0)
+                       (expt (+ k (sqrt (float (1+ (* k k)) (realpart v)))) (- v)))
+                      ((= n 1)
+                       (- (/ (* v (an 0 k v))
+                             (sqrt (float (1+ (* k k)) (realpart v))))))
+                      (t
+                       (/ (- (* (- (* v v) (expt (- n 2) 2)) (an (- n 2) k v))
+                             (* k (- n 1) (+ n n -3) (an (- n 1) k v)))
+                          (+ 1 (* k k))
+                          (- n 1)
+                          n)))))
+          (setf (gethash (list n k v) hash) result)
+          result))))
+
+;; SUM-AN computes the series
+;;
+;; sum(exp(-k*z)*a[n](k,v), k, 1, N)
+;;
+(defun sum-an (big-n n v z)
+  (let ((sum 0))
+    (loop for k from 1 upto big-n
+          do
+             (incf sum (* (exp (- (* k z)))
+                          (an n k v))))
+    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)
+(defun sum-ab (big-n v z)
+  (let ((eps (epsilon (realpart z))))
+    (an-clrhash)
+    (do* ((n 0 (+ 1 n))
+          (term (+ (* (alpha n z) (an n 0 v))
+                   (* (beta n z) (sum-an big-n n v z)))
+                (+ (* (alpha n z) (an n 0 v))
+                   (* (beta n z) (sum-an big-n n v z))))
+          (sum term (+ sum term)))
+         ((<= (abs term) (* eps (abs sum)))
+          sum)
+      (when nil
+        (format t "n = ~D~%" n)
+        (format t " term = ~S~%" term)
+        (format t " sum  = ~S~%" sum)))))
+
+;; Convert to iteration instead of this quick-and-dirty memoization?
+(let ((hash (make-hash-table :test 'equal)))
+  (defun %big-a-clrhash ()
+    (clrhash hash))
+  (defun %big-a-dump-hash ()
+    (maphash #'(lambda (k v)
+                 (format t "~S -> ~S~%" k v))
+             hash))
+  (defun %big-a (n v)
+    (or (gethash (list n v) hash)
+        (let ((result
+                (cond ((zerop n)
+                       (expt 2 (- v)))
+                      (t
+                       (* (%big-a (- n 1) v)
+                          (/ (* (+ v n n -2) (+ v n n -1))
+                             (* 4 n (+ n v))))))))
+          (setf (gethash (list n v) hash) result)
+          result))))
+
+;; Computes A[n](v) =
+;; (-1)^n*v*2^(-v)*pochhammer(v+n+1,n-1)/(2^(2*n)*n!)  If v is a
+;; negative integer -m, use A[n](-m) = (-1)^(m+1)*A[n-m](m) for n >=
+;; m.
+(defun big-a (n v)
+  (let ((m (ftruncate v)))
+    (cond ((and (= m v) (minusp m))
+           (if (< n m)
+               (%big-a n v)
+               (let ((result (%big-a (+ n m) v)))
+                 (if (oddp (truncate m))
+                     result
+                     (- result)))))
+          (t
+           (%big-a n v)))))
+
+;; I[n](t, z, v) = exp(-t*z)/t^(2*n+v-1) *
+;;                  integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf)
+;;
+;; Use the substitution u=1+s to get a new integral
+;;
+;; integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf)
+;;   = exp(t*z) * integrate(u^(-v-2*n)*exp(-t*u*z), u, 1, inf)
+;;   = exp(t*z)*t^(v+2*n-1)*z^(v+2*n-1)*incomplete_gamma_tail(1-v-2*n,t*z)
+;;
+;; The continued fraction for incomplete_gamma_tail(a,z) is
+;;
+;;   z^a*exp(-z)/CF
+;;
+;; So incomplete_gamma_tail(1-v-2*n, t*z) is
+;;
+;;   (t*z)^(1-v-2*n)*exp(-t*z)/CF
+;;
+;; which finally gives
+;;
+;; integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf)
+;;   = CF
+;;
+;; and I[n](t, z, v) = exp(-t*z)/t^(2*n+v-1)/CF
+(defun big-i (n t z v)
+  (/ (exp (- (* t z)))
+     (expt t (+ n n v -1))
+     (let* ((a (- 1 v n n))
+            (z-a (- z a)))
+       (lentz #'(lambda (n)
+                  (+ n n 1 z-a))
+              #'(lambda (n)
+                  (* n (- a n)))))))
+
+(defun sum-big-ia (big-n v z)
+  )
+
+(defun bessel-j (v z)
+  (let ((vv (ftruncate v)))
+    (cond ((= vv v)
+           ;; v is an integer
+           (integer-bessel-j-exp-arc v z))
+          (t
+           (let ((big-n 100)
+                 (vpi (* v (float-pi (realpart z)))))
+             (+ (integer-bessel-j-exp-arc v z)
+                (* z
+                   (/ (sin vpi) vpi)
+                   (+ (/ -1 z)
+                      (sum-ab big-n v z)))))))))
 
 (defun paris-series (v z n)