source: qd-bessel.lisp @ 7c5a31

Last change on this file since 7c5a31 was 7c5a31, checked in by Raymond Toy <rtoy@…>, 3 years ago

Correct some comments, remove unused code.

  • Property mode set to 100644
File size: 12.1 KB
Line 
1;;;; -*- Mode: lisp -*-
2;;;;
3;;;; Copyright (c) 2011 Raymond Toy
4;;;; Permission is hereby granted, free of charge, to any person
5;;;; obtaining a copy of this software and associated documentation
6;;;; files (the "Software"), to deal in the Software without
7;;;; restriction, including without limitation the rights to use,
8;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
9;;;; copies of the Software, and to permit persons to whom the
10;;;; Software is furnished to do so, subject to the following
11;;;; conditions:
12;;;;
13;;;; The above copyright notice and this permission notice shall be
14;;;; included in all copies or substantial portions of the Software.
15;;;;
16;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
17;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
18;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
19;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
20;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
21;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
22;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
23;;;; OTHER DEALINGS IN THE SOFTWARE.
24
25(in-package #:oct)
26
27;;; References:
28;;;
29;;; [1] Borwein, Borwein, Crandall, "Effective Laguerre Asymptotics",
30;;; http://people.reed.edu/~crandall/papers/Laguerre-f.pdf
31;;;
32;;; [2] Borwein, Borwein, Chan, "The Evaluation of Bessel Functions
33;;; via Exp-Arc Integrals", http://web.cs.dal.ca/~jborwein/bessel.pdf
34;;;
35
36(defvar *debug-exparc* nil)
37
38;; B[k](p) = 1/2^(k+3/2)*integrate(exp(-p*u)*u^(k-1/2),u,0,1)
39;;         = 1/2^(k+3/2)/p^(k+1/2)*integrate(t^(k-1/2)*exp(-t),t,0,p)
40;;         = 1/2^(k+3/2)/p^(k+1/2) * G(k+1/2, p)
41;;
42;; where G(a,z) is the lower incomplete gamma function.
43;;
44;; There is the continued fraction expansion for G(a,z) (see
45;; cf-incomplete-gamma in qd-gamma.lisp):
46;;
47;;  G(a,z) = z^a*exp(-z)/ CF
48;;
49;; So
50;;
51;;  B[k](p) = 1/2^(k+3/2)/p^(k+1/2)*p^(k+1/2)*exp(-p)/CF
52;;          = exp(-p)/2^(k+3/2)/CF
53;;
54;;
55;; Note also that [2] gives a recurrence relationship for B[k](p) in
56;; eq (2.6), but there is an error there.  The correct relationship is
57;;
58;;  B[k](p) = -exp(-p)/(p*sqrt(2)*2^(k+1)) + (k-1/2)*B[k-1](p)/(2*p)
59;;
60;; The paper is missing the division by p in the term containing
61;; B[k-1](p).  This is easily derived from the recurrence relationship
62;; for the (lower) incomplete gamma function.
63;;
64;; Note too that as k increases, the recurrence appears to be unstable
65;; and B[k](p) begins to increase even though it is strictly bounded.
66;; (This is also easy to see from the integral.)  Hence, we do not use
67;; the recursion.  However, it might be stable for use with
68;; double-float precision; this has not been tested.
69;;
70(defun bk (k p)
71  (/ (exp (- p))
72     (* (sqrt (float 2 (realpart p))) (ash 1 (+ k 1)))
73     (let ((a (float (+ k 1/2) (realpart p))))
74       (lentz #'(lambda (n)
75                  (+ n a))
76              #'(lambda (n)
77                  (if (evenp n)
78                      (* (ash n -1) p)
79                      (- (* (+ a (ash n -1)) p))))))))
80
81;; exp-arc I function, as given in the Laguerre paper
82;;
83;; I(p, q) = 4*exp(p) * sum(g[k](-2*%i*q)/(2*k)!*B[k](p), k, 0, inf)
84;;
85;; where g[k](p) = product(p^2+(2*j-1)^2, j, 1, k) and B[k](p) as above.
86;;
87;; For computation, note that g[k](p) = g[k-1](p) * (p^2 + (2*k-1)^2)
88;; and (2*k)! = (2*k-2)! * (2*k-1) * (2*k).  Then, let
89;;
90;;  R[k](p) = g[k](p)/(2*k)!
91;;
92;; Then
93;;
94;;  R[k](p) = g[k](p)/(2*k)!
95;;          = g[k-1](p)/(2*k-2)! * (p^2 + (2*k-1)^2)/((2*k-1)*(2*k)
96;;          = R[k-1](p) * (p^2 + (2*k-1)^2)/((2*k-1)*(2*k)
97;;
98;; In the exp-arc paper, the function is defined (equivalently) as
99;;
100;; I(p, q) = 2*%i*exp(p)/q * sum(r[2*k+1](-2*%i*q)/(2*k)!*B[k](p), k, 0, inf)
101;;
102;; where r[2*k+1](p) = p*product(p^2 + (2*j-1)^2, j, 1, k)
103;;
104;; Let's note some properties of I(p, q).
105;;
106;; I(-%i*z, v) = 2*%i*exp(-%i*z)/q * sum(r[2*k+1](-2*%i*v)/(2*k)!*B[k](-%i*z))
107;;
108;; Note thate B[k](-%i*z) = 1/2^(k+3/2)*integrate(exp(%i*z*u)*u^(k-1/2),u,0,1)
109;;                        = conj(B[k](%i*z).
110;;
111;; Hence I(-%i*z, v) = conj(I(%i*z, v)) when both z and v are real.
112;;
113;; Also note that when v is an integer of the form (2*m+1)/2, then
114;;   r[2*k+1](-2*%i*v) = r[2*k+1](-%i*(2*m+1))
115;;                     = -%i*(2*m+1)*product(-(2*m+1)^2+(2*j-1)^2, j, 1, k)
116;; so the product is zero when k >= m and the series I(p, q) is
117;; finite.
118(defun exp-arc-i (p q)
119  (let* ((sqrt2 (sqrt (float 2 (realpart p))))
120         (exp/p/sqrt2 (/ (exp (- p)) p sqrt2))
121         (v (* #c(0 -2) q))
122         (v2 (expt v 2))
123         (eps (epsilon (realpart p))))
124    (when *debug-exparc*
125      (format t "sqrt2 = ~S~%" sqrt2)
126      (format t "exp/p/sqrt2 = ~S~%" exp/p/sqrt2))
127    (do* ((k 0 (1+ k))
128          (bk (/ (incomplete-gamma 1/2 p)
129                 2 sqrt2 (sqrt p))
130              (- (/ (* bk (- k 1/2)) 2 p)
131                 (/ exp/p/sqrt2 (ash 1 (+ k 1)))))
132          ;; ratio[k] = r[2*k+1](v)/(2*k)!.
133          ;; r[1] = v and r[2*k+1](v) = r[2*k-1](v)*(v^2 + (2*k-1)^2)
134          ;; ratio[0] = v
135          ;; and ratio[k] = r[2*k-1](v)*(v^2+(2*k-1)^2) / ((2*k-2)! * (2*k-1) * 2*k)
136          ;;              = ratio[k]*(v^2+(2*k-1)^2)/((2*k-1) * 2 * k)
137          (ratio v
138                 (* ratio (/ (+ v2 (expt (1- (* 2 k)) 2))
139                             (* 2 k (1- (* 2 k))))))
140          (term (* ratio bk)
141                (* ratio bk))
142          (sum term (+ sum term)))
143         ((< (abs term) (* (abs sum) eps))
144          (* sum #c(0 2) (/ (exp p) q)))
145      (when *debug-exparc*
146        (format t "k      = ~D~%" k)
147        (format t " bk    = ~S~%" bk)
148        (format t " ratio = ~S~%" ratio)
149        (format t " term  = ~S~%" term)
150        (format t " sum   - ~S~%" sum)))))
151
152(defun exp-arc-i-2 (p q)
153  (let* ((v (* #c(0 -2) q))
154         (v2 (expt v 2))
155         (eps (epsilon (realpart p))))
156    (do* ((k 0 (1+ k))
157          (bk (bk 0 p)
158              (bk k p))
159          (ratio v
160                 (* ratio (/ (+ v2 (expt (1- (* 2 k)) 2))
161                             (* 2 k (1- (* 2 k))))))
162          (term (* ratio bk)
163                (* ratio bk))
164          (sum term (+ sum term)))
165         ((< (abs term) (* (abs sum) eps))
166          (when *debug-exparc*
167            (format t "Final k= ~D~%" k)
168            (format t " bk    = ~S~%" bk)
169            (format t " ratio = ~S~%" ratio)
170            (format t " term  = ~S~%" term)
171            (format t " sum   - ~S~%" sum))
172          (* sum #c(0 2) (/ (exp p) q)))
173      (when *debug-exparc*
174        (format t "k      = ~D~%" k)
175        (format t " bk    = ~S~%" bk)
176        (format t " ratio = ~S~%" ratio)
177        (format t " term  = ~S~%" term)
178        (format t " sum   - ~S~%" sum)))))
179
180
181;;
182(defun integer-bessel-j-exp-arc (v z)
183  (let* ((iz (* #c(0 1) z))
184         (i+ (exp-arc-i-2 iz v)))
185    (cond ((= v (ftruncate v))
186           ;; We can simplify the result
187           (let ((c (cis (* v (float-pi i+) -1/2))))
188             (/ (+ (* c i+)
189                   (* (conjugate c) (conjugate i+)))
190                (float-pi i+)
191                2)))
192          (t
193           (let ((i- (exp-arc-i-2 (- iz ) v)))
194             (/ (+ (* (cis (* v (float-pi i+) -1/2))
195                      i+)
196                   (* (cis (* v (float-pi i+) 1/2))
197                      i-))
198                (float-pi i+)
199                2))))))
200
201;; alpha[n](z) = integrate(exp(-z*s)*s^n, s, 0, 1/2)
202;; beta[n](z)  = integrate(exp(-z*s)*s^n, s, -1/2, 1/2)
203;;
204;; The recurrence in [2] is
205;;
206;; alpha[n](z) = - exp(-z/2)/2^n/z + n/z*alpha[n-1](z)
207;; beta[n]z)   = ((-1)^n*exp(z/2)-exp(-z/2))/2^n/z + n/z*beta[n-1](z)
208;;
209;; We also note that
210;;
211;; alpha[n](z) = G(n+1,z/2)/z^(n+1)
212;; beta[n](z)  = G(n+1,z/2)/z^(n+1) - G(n+1,-z/2)/z^(n+1)
213
214(defun alpha (n z)
215  (let ((n (float n (realpart z))))
216    (/ (cf-incomplete-gamma (1+ n) (/ z 2))
217       (expt z (1+ n)))))
218
219(defun beta (n z)
220  (let ((n (float n (realpart z))))
221    (/ (- (cf-incomplete-gamma (1+ n) (/ z 2))
222          (cf-incomplete-gamma (1+ n) (/ z -2)))
223       (expt z (1+ n)))))
224
225;; a[0](k,v) := (k+sqrt(k^2+1))^(-v);
226;; a[1](k,v) := -v*a[0](k,v)/sqrt(k^2+1);
227;; 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));
228
229;; Convert this to iteration instead of using this quick-and-dirty
230;; memoization?
231(let ((hash (make-hash-table :test 'equal)))
232  (defun an-clrhash ()
233    (clrhash hash))
234  (defun an-dump-hash ()
235    (maphash #'(lambda (k v)
236                 (format t "~S -> ~S~%" k v))
237             hash))
238  (defun an (n k v)
239    (or (gethash (list n k v) hash)
240        (let ((result
241                (cond ((= n 0)
242                       (expt (+ k (sqrt (float (1+ (* k k)) (realpart v)))) (- v)))
243                      ((= n 1)
244                       (- (/ (* v (an 0 k v))
245                             (sqrt (float (1+ (* k k)) (realpart v))))))
246                      (t
247                       (/ (- (* (- (* v v) (expt (- n 2) 2)) (an (- n 2) k v))
248                             (* k (- n 1) (+ n n -3) (an (- n 1) k v)))
249                          (+ 1 (* k k))
250                          (- n 1)
251                          n)))))
252          (setf (gethash (list n k v) hash) result)
253          result))))
254
255;; SUM-AN computes the series
256;;
257;; sum(exp(-k*z)*a[n](k,v), k, 1, N)
258;;
259(defun sum-an (big-n n v z)
260  (let ((sum 0))
261    (loop for k from 1 upto big-n
262          do
263             (incf sum (* (exp (- (* k z)))
264                          (an n k v))))
265    sum))
266
267;; SUM-AB computes the series
268;;
269;; sum(alpha[n](z)*a[n](0,v) + beta[n](z)*sum_an(N, n, v, z), n, 0, inf)
270(defun sum-ab (big-n v z)
271  (let ((eps (epsilon (realpart z))))
272    (an-clrhash)
273    (do* ((n 0 (+ 1 n))
274          (term (+ (* (alpha n z) (an n 0 v))
275                   (* (beta n z) (sum-an big-n n v z)))
276                (+ (* (alpha n z) (an n 0 v))
277                   (* (beta n z) (sum-an big-n n v z))))
278          (sum term (+ sum term)))
279         ((<= (abs term) (* eps (abs sum)))
280          sum)
281      (when nil
282        (format t "n = ~D~%" n)
283        (format t " term = ~S~%" term)
284        (format t " sum  = ~S~%" sum)))))
285
286;; Convert to iteration instead of this quick-and-dirty memoization?
287(let ((hash (make-hash-table :test 'equal)))
288  (defun %big-a-clrhash ()
289    (clrhash hash))
290  (defun %big-a-dump-hash ()
291    (maphash #'(lambda (k v)
292                 (format t "~S -> ~S~%" k v))
293             hash))
294  (defun %big-a (n v)
295    (or (gethash (list n v) hash)
296        (let ((result
297                (cond ((zerop n)
298                       (expt 2 (- v)))
299                      (t
300                       (* (%big-a (- n 1) v)
301                          (/ (* (+ v n n -2) (+ v n n -1))
302                             (* 4 n (+ n v))))))))
303          (setf (gethash (list n v) hash) result)
304          result))))
305
306;; Computes A[n](v) =
307;; (-1)^n*v*2^(-v)*pochhammer(v+n+1,n-1)/(2^(2*n)*n!)  If v is a
308;; negative integer -m, use A[n](-m) = (-1)^(m+1)*A[n-m](m) for n >=
309;; m.
310(defun big-a (n v)
311  (let ((m (ftruncate v)))
312    (cond ((and (= m v) (minusp m))
313           (if (< n m)
314               (%big-a n v)
315               (let ((result (%big-a (+ n m) (- v))))
316                 (if (oddp (truncate m))
317                     result
318                     (- result)))))
319          (t
320           (%big-a n v)))))
321
322;; I[n](t, z, v) = exp(-t*z)/t^(2*n+v-1) *
323;;                  integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf)
324;;
325;; Use the substitution u=1+s to get a new integral
326;;
327;;   integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf)
328;;     = exp(t*z) * integrate(u^(-v-2*n)*exp(-t*u*z), u, 1, inf)
329;;     = exp(t*z)*t^(v+2*n-1)*z^(v+2*n-1)*incomplete_gamma_tail(1-v-2*n,t*z)
330;;
331;; Thus,
332;;
333;;   I[n](t, z, v) = z^(v+2*n-1)*incomplete_gamma_tail(1-v-2*n,t*z)
334;;
335(defun big-i (n theta z v)
336  (let* ((a (- 1 v n n)))
337    (* (expt z (- a))
338       (incomplete-gamma-tail a (* theta z)))))
339
340(defun sum-big-ia (big-n v z)
341  (let ((big-n-1/2 (+ big-n 1/2))
342        (eps (epsilon z)))
343    (do* ((n 0 (1+ n))
344          (term (* (big-a 0 v)
345                   (big-i 0 big-n-1/2 z v))
346                (* (big-a n v)
347                   (big-i n big-n-1/2 z v)))
348          (sum term (+ sum term)))
349         ((<= (abs term) (* eps (abs sum)))
350          sum)
351      #+nil
352      (progn
353        (format t "n = ~D~%" n)
354        (format t " term = ~S~%" term)
355        (format t " sum  = ~S~%" sum)))))
356
357;; Series for bessel J:
358;;
359;; (z/2)^v*sum((-1)^k/Gamma(k+v+1)/k!*(z^2//4)^k, k, 0, inf)
360(defun s-bessel-j (v z)
361  (with-floating-point-contagion (v z)
362    (let ((z2/4 (* z z 1/4))
363          (eps (epsilon z)))
364      (do* ((k 0 (+ 1 k))
365            (f (gamma (+ v 1))
366               (* k (+ v k)))
367            (term (/ f)
368                  (/ (* (- term) z2/4) f))
369            (sum term (+ sum term)))
370           ((<= (abs term) (* eps (abs sum)))
371            (* sum (expt (* z 1/2) v)))
372        #+nil
373        (progn
374          (format t "k = ~D~%" k)
375          (format t " f    = ~S~%" f)
376          (format t " term = ~S~%" term)
377          (format t " sum  = ~S~%" sum))))))
378 
379(defun bessel-j (v z)
380  (let ((vv (ftruncate v)))
381    (cond ((= vv v)
382           ;; v is an integer
383           (integer-bessel-j-exp-arc v z))
384          (t
385           (let ((big-n 100)
386                 (vpi (* v (float-pi (realpart z)))))
387             (+ (integer-bessel-j-exp-arc v z)
388                (* z
389                   (/ (sin vpi) vpi)
390                   (+ (/ -1 z)
391                      (sum-ab big-n v z)
392                      (sum-big-ia big-n v z)))))))))
393
394(defun paris-series (v z n)
395  (labels ((pochhammer (a k)
396             (/ (gamma (+ a k))
397                (gamma a)))
398           (a (v k)
399             (* (/ (pochhammer (+ 1/2 v) k)
400                   (gamma (float (1+ k) z)))
401                (pochhammer (- 1/2 v) k))))
402    (* (loop for k from 0 below n
403             sum (* (/ (a v k)
404                       (expt (* 2 z) k))
405                    (/ (cf-incomplete-gamma (+ k v 1/2) (* 2 z))
406                       (gamma (+ k v 1/2)))))
407       (/ (exp z)
408          (sqrt (* 2 (float-pi z) z))))))
409
Note: See TracBrowser for help on using the repository browser.