source: qd-bessel.lisp @ e6a3f2

Last change on this file since e6a3f2 was e6a3f2, checked in by Raymond Toy <toy.raymond@…>, 3 years ago

sum-an was summing one too many terms.

  • Property mode set to 100644
File size: 16.6 KB
RevLine 
[4c1ed0]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)
[7c5a31]40;;         = 1/2^(k+3/2)/p^(k+1/2) * G(k+1/2, p)
[4c1ed0]41;;
[a5a4c7]42;; where G(a,z) is the lower incomplete gamma function.
[4c1ed0]43;;
[a5a4c7]44;; There is the continued fraction expansion for G(a,z) (see
[4c1ed0]45;; cf-incomplete-gamma in qd-gamma.lisp):
46;;
[a5a4c7]47;;  G(a,z) = z^a*exp(-z)/ CF
[4c1ed0]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;;
[1d404a]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;;
[4c1ed0]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
[64c424]81;; Use the recursion
82(defun bk-iter (k p old-bk)
83  (with-floating-point-contagion (p old-bk)
84    (if (zerop k)
85        (* (sqrt (/ (float-pi p) 8))
86           (let ((rp (sqrt p)))
87             (/ (erf rp)
88                rp)))
89        (- (* (- k 1/2)
90              (/ old-bk (* 2 p)))
91           (/ (exp (- p))
92              p
93              (ash 1 (+ k 1))
94              (sqrt  (float 2 (realpart p))))))))
95
[4c1ed0]96;; exp-arc I function, as given in the Laguerre paper
97;;
98;; I(p, q) = 4*exp(p) * sum(g[k](-2*%i*q)/(2*k)!*B[k](p), k, 0, inf)
99;;
100;; where g[k](p) = product(p^2+(2*j-1)^2, j, 1, k) and B[k](p) as above.
101;;
102;; For computation, note that g[k](p) = g[k-1](p) * (p^2 + (2*k-1)^2)
103;; and (2*k)! = (2*k-2)! * (2*k-1) * (2*k).  Then, let
104;;
105;;  R[k](p) = g[k](p)/(2*k)!
106;;
107;; Then
108;;
109;;  R[k](p) = g[k](p)/(2*k)!
110;;          = g[k-1](p)/(2*k-2)! * (p^2 + (2*k-1)^2)/((2*k-1)*(2*k)
111;;          = R[k-1](p) * (p^2 + (2*k-1)^2)/((2*k-1)*(2*k)
112;;
113;; In the exp-arc paper, the function is defined (equivalently) as
114;;
115;; I(p, q) = 2*%i*exp(p)/q * sum(r[2*k+1](-2*%i*q)/(2*k)!*B[k](p), k, 0, inf)
116;;
117;; where r[2*k+1](p) = p*product(p^2 + (2*j-1)^2, j, 1, k)
118;;
119;; Let's note some properties of I(p, q).
120;;
121;; I(-%i*z, v) = 2*%i*exp(-%i*z)/q * sum(r[2*k+1](-2*%i*v)/(2*k)!*B[k](-%i*z))
122;;
123;; Note thate B[k](-%i*z) = 1/2^(k+3/2)*integrate(exp(%i*z*u)*u^(k-1/2),u,0,1)
124;;                        = conj(B[k](%i*z).
125;;
126;; Hence I(-%i*z, v) = conj(I(%i*z, v)) when both z and v are real.
[7c5a31]127;;
128;; Also note that when v is an integer of the form (2*m+1)/2, then
129;;   r[2*k+1](-2*%i*v) = r[2*k+1](-%i*(2*m+1))
130;;                     = -%i*(2*m+1)*product(-(2*m+1)^2+(2*j-1)^2, j, 1, k)
131;; so the product is zero when k >= m and the series I(p, q) is
132;; finite.
[4c1ed0]133(defun exp-arc-i (p q)
134  (let* ((sqrt2 (sqrt (float 2 (realpart p))))
135         (exp/p/sqrt2 (/ (exp (- p)) p sqrt2))
136         (v (* #c(0 -2) q))
137         (v2 (expt v 2))
138         (eps (epsilon (realpart p))))
139    (when *debug-exparc*
140      (format t "sqrt2 = ~S~%" sqrt2)
141      (format t "exp/p/sqrt2 = ~S~%" exp/p/sqrt2))
142    (do* ((k 0 (1+ k))
143          (bk (/ (incomplete-gamma 1/2 p)
144                 2 sqrt2 (sqrt p))
145              (- (/ (* bk (- k 1/2)) 2 p)
146                 (/ exp/p/sqrt2 (ash 1 (+ k 1)))))
147          ;; ratio[k] = r[2*k+1](v)/(2*k)!.
148          ;; r[1] = v and r[2*k+1](v) = r[2*k-1](v)*(v^2 + (2*k-1)^2)
149          ;; ratio[0] = v
150          ;; and ratio[k] = r[2*k-1](v)*(v^2+(2*k-1)^2) / ((2*k-2)! * (2*k-1) * 2*k)
151          ;;              = ratio[k]*(v^2+(2*k-1)^2)/((2*k-1) * 2 * k)
152          (ratio v
153                 (* ratio (/ (+ v2 (expt (1- (* 2 k)) 2))
154                             (* 2 k (1- (* 2 k))))))
155          (term (* ratio bk)
156                (* ratio bk))
157          (sum term (+ sum term)))
158         ((< (abs term) (* (abs sum) eps))
159          (* sum #c(0 2) (/ (exp p) q)))
160      (when *debug-exparc*
161        (format t "k      = ~D~%" k)
162        (format t " bk    = ~S~%" bk)
163        (format t " ratio = ~S~%" ratio)
164        (format t " term  = ~S~%" term)
165        (format t " sum   - ~S~%" sum)))))
166
167(defun exp-arc-i-2 (p q)
[7c5a31]168  (let* ((v (* #c(0 -2) q))
[4c1ed0]169         (v2 (expt v 2))
170         (eps (epsilon (realpart p))))
171    (do* ((k 0 (1+ k))
172          (bk (bk 0 p)
173              (bk k p))
[9fd2eb]174          ;; Compute g[k](p)/(2*k)!, not r[2*k+1](p)/(2*k)!
175          (ratio 1
[4c1ed0]176                 (* ratio (/ (+ v2 (expt (1- (* 2 k)) 2))
177                             (* 2 k (1- (* 2 k))))))
178          (term (* ratio bk)
179                (* ratio bk))
180          (sum term (+ sum term)))
181         ((< (abs term) (* (abs sum) eps))
[7c5a31]182          (when *debug-exparc*
183            (format t "Final k= ~D~%" k)
184            (format t " bk    = ~S~%" bk)
185            (format t " ratio = ~S~%" ratio)
186            (format t " term  = ~S~%" term)
187            (format t " sum   - ~S~%" sum))
[9fd2eb]188          (* sum 4 (exp p)))
[4c1ed0]189      (when *debug-exparc*
190        (format t "k      = ~D~%" k)
191        (format t " bk    = ~S~%" bk)
192        (format t " ratio = ~S~%" ratio)
193        (format t " term  = ~S~%" term)
194        (format t " sum   - ~S~%" sum)))))
195
[64c424]196(defun exp-arc-i-3 (p q)
197  (let* ((v (* #c(0 -2) q))
198         (v2 (expt v 2))
199         (eps (epsilon (realpart p))))
200    (do* ((k 0 (1+ k))
201          (bk (bk 0 p)
202              (bk-iter k p bk))
203          ;; Compute g[k](p)/(2*k)!, not r[2*k+1](p)/(2*k)!
204          (ratio 1
205                 (* ratio (/ (+ v2 (expt (1- (* 2 k)) 2))
206                             (* 2 k (1- (* 2 k))))))
207          (term (* ratio bk)
208                (* ratio bk))
209          (sum term (+ sum term)))
210         ((< (abs term) (* (abs sum) eps))
211          (when *debug-exparc*
212            (format t "Final k= ~D~%" k)
213            (format t " bk    = ~S~%" bk)
214            (format t " ratio = ~S~%" ratio)
215            (format t " term  = ~S~%" term)
216            (format t " sum   - ~S~%" sum))
217          (* sum 4 (exp p)))
218      (when *debug-exparc*
219        (format t "k      = ~D~%" k)
220        (format t " bk    = ~S~%" bk)
221        (format t " ratio = ~S~%" ratio)
222        (format t " term  = ~S~%" term)
223        (format t " sum   - ~S~%" sum)))))
224
[4c1ed0]225
[78801d]226;; Not really just for Bessel J for integer orders, but in that case,
227;; this is all that's needed to compute Bessel J.  For other values,
228;; this is just part of the computation needed.
[4c1ed0]229;;
[78801d]230;; Compute
231;;
232;;  1/(2*%pi) * (exp(-%i*v*%pi/2) * I(%i*z, v) + exp(%i*v*%pi/2) * I(-%i*z, v))
[1d404a]233(defun integer-bessel-j-exp-arc (v z)
[4c1ed0]234  (let* ((iz (* #c(0 1) z))
[e9cd1a]235         (i+ (exp-arc-i-2 iz v)))
[235ac2]236    (cond ((and (= v (ftruncate v)) (realp z))
[e9cd1a]237           ;; We can simplify the result
[235ac2]238           (let ((c (exp (* v (float-pi i+) #c(0 -1/2)))))
[e9cd1a]239             (/ (+ (* c i+)
240                   (* (conjugate c) (conjugate i+)))
241                (float-pi i+)
242                2)))
243          (t
244           (let ((i- (exp-arc-i-2 (- iz ) v)))
[235ac2]245             (/ (+ (* (exp (* v (float-pi i+) #c(0 -1/2)))
[e9cd1a]246                      i+)
[235ac2]247                   (* (exp (* v (float-pi i+) #c(0 1/2)))
[e9cd1a]248                      i-))
249                (float-pi i+)
250                2))))))
[a5a4c7]251
252;; alpha[n](z) = integrate(exp(-z*s)*s^n, s, 0, 1/2)
253;; beta[n](z)  = integrate(exp(-z*s)*s^n, s, -1/2, 1/2)
254;;
255;; The recurrence in [2] is
256;;
257;; alpha[n](z) = - exp(-z/2)/2^n/z + n/z*alpha[n-1](z)
258;; beta[n]z)   = ((-1)^n*exp(z/2)-exp(-z/2))/2^n/z + n/z*beta[n-1](z)
[64c424]259;;             = (-1)^n/(2^n)*2*sinh(z/2)/z + n/z*beta[n-1](z)
[a5a4c7]260;;
261;; We also note that
262;;
263;; alpha[n](z) = G(n+1,z/2)/z^(n+1)
264;; beta[n](z)  = G(n+1,z/2)/z^(n+1) - G(n+1,-z/2)/z^(n+1)
265
266(defun alpha (n z)
267  (let ((n (float n (realpart z))))
[c0e12d]268    (/ (incomplete-gamma (1+ n) (/ z 2))
[a5a4c7]269       (expt z (1+ n)))))
270
[64c424]271(defun alpha-iter (n z alpha-old)
272  (if (zerop n)
273      ;; (1- exp(-z/2))/z.
274      (/ (- 1 (exp (* z -1/2)))
275         z)
276      (- (* (/ n z) alpha-old)
277         (/ (exp (- (* z 1/2)))
278            z
279            (ash 1 n)))))
280
[a5a4c7]281(defun beta (n z)
282  (let ((n (float n (realpart z))))
[c0e12d]283    (/ (- (incomplete-gamma (1+ n) (/ z 2))
284          (incomplete-gamma (1+ n) (/ z -2)))
[a5a4c7]285       (expt z (1+ n)))))
286
[64c424]287(defun beta-iter (n z old-beta)
288  (if (zerop n)
289      ;; integrate(exp(-z*s),s,-1/2,1/2)
290      ;;   = (exp(z/2)-exp(-z/2)/z
291      ;;   = 2*sinh(z/2)/z
292      ;;   = sinh(z/2)/(z/2)
293      (* 2 (/ (sinh (* 1/2 z)) z))
294      (+ (* n (/ old-beta z))
295         (* (/ (sinh (* 1/2 z)) (* 1/2 z))
296            (scale-float (float (if (evenp n) 1 -1) (realpart z)) (- n))))))
297
298
[a5a4c7]299;; a[0](k,v) := (k+sqrt(k^2+1))^(-v);
300;; a[1](k,v) := -v*a[0](k,v)/sqrt(k^2+1);
301;; 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));
302
303;; Convert this to iteration instead of using this quick-and-dirty
304;; memoization?
305(let ((hash (make-hash-table :test 'equal)))
306  (defun an-clrhash ()
307    (clrhash hash))
308  (defun an-dump-hash ()
309    (maphash #'(lambda (k v)
310                 (format t "~S -> ~S~%" k v))
311             hash))
312  (defun an (n k v)
313    (or (gethash (list n k v) hash)
314        (let ((result
315                (cond ((= n 0)
316                       (expt (+ k (sqrt (float (1+ (* k k)) (realpart v)))) (- v)))
317                      ((= n 1)
318                       (- (/ (* v (an 0 k v))
319                             (sqrt (float (1+ (* k k)) (realpart v))))))
320                      (t
321                       (/ (- (* (- (* v v) (expt (- n 2) 2)) (an (- n 2) k v))
322                             (* k (- n 1) (+ n n -3) (an (- n 1) k v)))
323                          (+ 1 (* k k))
324                          (- n 1)
325                          n)))))
326          (setf (gethash (list n k v) hash) result)
327          result))))
328
329;; SUM-AN computes the series
330;;
331;; sum(exp(-k*z)*a[n](k,v), k, 1, N)
332;;
[78801d]333#+nil
[a5a4c7]334(defun sum-an (big-n n v z)
335  (let ((sum 0))
336    (loop for k from 1 upto big-n
337          do
338             (incf sum (* (exp (- (* k z)))
339                          (an n k v))))
340    sum))
341
[78801d]342;; Like above, but we just stop when the terms no longer contribute to
343;; the sum.
344(defun sum-an (big-n n v z)
345  (let ((eps (epsilon (realpart z))))
346    (do* ((k 1 (+ 1 k))
347          (term (* (exp (- (* k z)))
348                   (an n k v))
349                (* (exp (- (* k z)))
350                   (an n k v)))
351          (sum term (+ sum term)))
352         ((or (<= (abs term) (* eps (abs sum)))
[e6a3f2]353              (>= k big-n))
[78801d]354          sum))))
355
[a5a4c7]356;; SUM-AB computes the series
357;;
358;; sum(alpha[n](z)*a[n](0,v) + beta[n](z)*sum_an(N, n, v, z), n, 0, inf)
359(defun sum-ab (big-n v z)
360  (let ((eps (epsilon (realpart z))))
361    (an-clrhash)
362    (do* ((n 0 (+ 1 n))
363          (term (+ (* (alpha n z) (an n 0 v))
364                   (* (beta n z) (sum-an big-n n v z)))
365                (+ (* (alpha n z) (an n 0 v))
366                   (* (beta n z) (sum-an big-n n v z))))
367          (sum term (+ sum term)))
368         ((<= (abs term) (* eps (abs sum)))
369          sum)
370      (when nil
371        (format t "n = ~D~%" n)
372        (format t " term = ~S~%" term)
373        (format t " sum  = ~S~%" sum)))))
374
[64c424]375(defun sum-ab-2 (big-n v z)
376  (let ((eps (epsilon (realpart z))))
377    (an-clrhash)
378    (do* ((n 0 (+ 1 n))
379          (alphan (alpha-iter 0 z 0)
380                  (alpha-iter n z alphan))
381          (betan (beta-iter 0 z 0)
382                 (beta-iter n z betan))
383          (term (+ (* alphan (an n 0 v))
384                   (* betan (sum-an big-n n v z)))
385                (+ (* alphan (an n 0 v))
386                   (* betan (sum-an big-n n v z))))
387          (sum term (+ sum term)))
388         ((<= (abs term) (* eps (abs sum)))
389          sum)
390      (when nil
391        (format t "n = ~D~%" n)
392        (format t " term = ~S~%" term)
393        (format t " sum  = ~S~%" sum)))))
394
[a5a4c7]395;; Convert to iteration instead of this quick-and-dirty memoization?
396(let ((hash (make-hash-table :test 'equal)))
397  (defun %big-a-clrhash ()
398    (clrhash hash))
399  (defun %big-a-dump-hash ()
400    (maphash #'(lambda (k v)
401                 (format t "~S -> ~S~%" k v))
402             hash))
403  (defun %big-a (n v)
404    (or (gethash (list n v) hash)
405        (let ((result
406                (cond ((zerop n)
407                       (expt 2 (- v)))
408                      (t
409                       (* (%big-a (- n 1) v)
410                          (/ (* (+ v n n -2) (+ v n n -1))
411                             (* 4 n (+ n v))))))))
412          (setf (gethash (list n v) hash) result)
413          result))))
414
415;; Computes A[n](v) =
416;; (-1)^n*v*2^(-v)*pochhammer(v+n+1,n-1)/(2^(2*n)*n!)  If v is a
417;; negative integer -m, use A[n](-m) = (-1)^(m+1)*A[n-m](m) for n >=
418;; m.
419(defun big-a (n v)
420  (let ((m (ftruncate v)))
421    (cond ((and (= m v) (minusp m))
422           (if (< n m)
423               (%big-a n v)
[c7fc98]424               (let ((result (%big-a (+ n m) (- v))))
[a5a4c7]425                 (if (oddp (truncate m))
426                     result
427                     (- result)))))
428          (t
429           (%big-a n v)))))
430
431;; I[n](t, z, v) = exp(-t*z)/t^(2*n+v-1) *
432;;                  integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf)
433;;
434;; Use the substitution u=1+s to get a new integral
435;;
[c7fc98]436;;   integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf)
437;;     = exp(t*z) * integrate(u^(-v-2*n)*exp(-t*u*z), u, 1, inf)
438;;     = exp(t*z)*t^(v+2*n-1)*z^(v+2*n-1)*incomplete_gamma_tail(1-v-2*n,t*z)
[a5a4c7]439;;
[c7fc98]440;; Thus,
[a5a4c7]441;;
[c7fc98]442;;   I[n](t, z, v) = z^(v+2*n-1)*incomplete_gamma_tail(1-v-2*n,t*z)
[a5a4c7]443;;
[95aa58]444(defun big-i (n theta z v)
[c7fc98]445  (let* ((a (- 1 v n n)))
446    (* (expt z (- a))
447       (incomplete-gamma-tail a (* theta z)))))
[a5a4c7]448
449(defun sum-big-ia (big-n v z)
[f20338]450  (let ((big-n-1/2 (+ big-n 1/2))
451        (eps (epsilon z)))
452    (do* ((n 0 (1+ n))
453          (term (* (big-a 0 v)
454                   (big-i 0 big-n-1/2 z v))
455                (* (big-a n v)
456                   (big-i n big-n-1/2 z v)))
457          (sum term (+ sum term)))
458         ((<= (abs term) (* eps (abs sum)))
459          sum)
460      #+nil
461      (progn
462        (format t "n = ~D~%" n)
463        (format t " term = ~S~%" term)
464        (format t " sum  = ~S~%" sum)))))
[a5a4c7]465
[f20338]466;; Series for bessel J:
467;;
468;; (z/2)^v*sum((-1)^k/Gamma(k+v+1)/k!*(z^2//4)^k, k, 0, inf)
469(defun s-bessel-j (v z)
470  (with-floating-point-contagion (v z)
471    (let ((z2/4 (* z z 1/4))
472          (eps (epsilon z)))
473      (do* ((k 0 (+ 1 k))
474            (f (gamma (+ v 1))
[e9cd1a]475               (* k (+ v k)))
[f20338]476            (term (/ f)
477                  (/ (* (- term) z2/4) f))
478            (sum term (+ sum term)))
479           ((<= (abs term) (* eps (abs sum)))
480            (* sum (expt (* z 1/2) v)))
[e9cd1a]481        #+nil
482        (progn
483          (format t "k = ~D~%" k)
484          (format t " f    = ~S~%" f)
485          (format t " term = ~S~%" term)
486          (format t " sum  = ~S~%" sum))))))
[d795ba]487
[6ab522]488;;
[d795ba]489;; TODO:
490;;  o For |z| <= 1 use the series.
491;;  o Currently accuracy is not good for large z and half-integer
492;;    order.
493;;  o For real v and z, return a real number instead of complex.
494;;  o Handle the case of Re(z) < 0. (The formulas are for Re(z) > 0:
495;;    bessel_j(v,z*exp(m*%pi*%i)) = exp(m*v*%pi*%i)*bessel_j(v, z)
496;;  o The paper suggests using
497;;      bessel_i(v,z) = exp(-v*%pi*%i/2)*bessel_j(v, %i*z)
498;;    when Im(z) >> Re(z)
499;;
[8c5195]500(defvar *big-n* 100)
[a5a4c7]501(defun bessel-j (v z)
502  (let ((vv (ftruncate v)))
[9fd2eb]503    ;; Clear the caches for now.
504    (an-clrhash)
505    (%big-a-clrhash)
[235ac2]506    (cond ((and (= vv v) (realp z))
507           ;; v is an integer and z is real
[a5a4c7]508           (integer-bessel-j-exp-arc v z))
509          (t
[d795ba]510           ;; Need to fine-tune the value of big-n.
[8c5195]511           (let ((big-n *big-n*)
[a5a4c7]512                 (vpi (* v (float-pi (realpart z)))))
513             (+ (integer-bessel-j-exp-arc v z)
[235ac2]514                (if (= vv v)
515                    0
516                    (* z
517                       (/ (sin vpi) vpi)
518                       (+ (/ -1 z)
519                          (sum-ab big-n v z)
520                          (sum-big-ia big-n v z))))))))))
[6ab522]521
522;; Bessel Y
523;;
524;; bessel_y(v, z) = 1/(2*%pi*%i)*(exp(-%i*v*%pi/2)*I(%i*v,z) - exp(%i*v*%pi/2)*I(-%i*z, v))
525;;                   + z/v/%pi*((1-cos(v*%pi)/z) + S(N,z,v)*cos(v*%pi)-S(N,z,-v))
526;;
527;; where
528;;
529;;   S(N,z,v) = sum(alpha[n](z)*a[n](0,v) + beta[n](z)*sum(exp(-k*z)*a[n](k,v),k,1,N),n,0,inf)
530;;               + sum(A[n](v)*I[n](N+1/2,z,v),n,0,inf)
531;;
532(defun bessel-y (v z)
533  (flet ((ipart (v z)
534           (let* ((iz (* #c(0 1) z))
535                  (c+ (exp (* v (float-pi z) 1/2)))
536                  (c- (exp (* v (float-pi z) -1/2)))
537                  (i+ (exp-arc-i-2 iz v))
538                  (i- (exp-arc-i-2 (- iz) v)))
539             (/ (- (* c- i+) (* c+ i-))
540                (* #c(0 2) (float-pi z)))))
541         (s (big-n z v)
542           (+ (sum-ab big-n v z)
543              (sum-big-ia big-n v z))))
544    (let* ((big-n 100)
545           (vpi (* v (float-pi z)))
546           (c (cos vpi)))
547      (+ (ipart v z)
548         (* (/ z vpi)
549            (+ (/ (- 1 c)
550                  z)
551               (* c
552                  (s big-n z v))
553               (- (s big-n z (- v)))))))))
554           
555 
[4c1ed0]556
557(defun paris-series (v z n)
558  (labels ((pochhammer (a k)
559             (/ (gamma (+ a k))
560                (gamma a)))
561           (a (v k)
562             (* (/ (pochhammer (+ 1/2 v) k)
563                   (gamma (float (1+ k) z)))
564                (pochhammer (- 1/2 v) k))))
565    (* (loop for k from 0 below n
566             sum (* (/ (a v k)
567                       (expt (* 2 z) k))
568                    (/ (cf-incomplete-gamma (+ k v 1/2) (* 2 z))
569                       (gamma (+ k v 1/2)))))
570       (/ (exp z)
571          (sqrt (* 2 (float-pi z) z))))))
572
Note: See TracBrowser for help on using the repository browser.