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 | |
---|