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