| 186 | |
| 187 | ;; alpha[n](z) = integrate(exp(-z*s)*s^n, s, 0, 1/2) |
| 188 | ;; beta[n](z) = integrate(exp(-z*s)*s^n, s, -1/2, 1/2) |
| 189 | ;; |
| 190 | ;; The recurrence in [2] is |
| 191 | ;; |
| 192 | ;; alpha[n](z) = - exp(-z/2)/2^n/z + n/z*alpha[n-1](z) |
| 193 | ;; beta[n]z) = ((-1)^n*exp(z/2)-exp(-z/2))/2^n/z + n/z*beta[n-1](z) |
| 194 | ;; |
| 195 | ;; We also note that |
| 196 | ;; |
| 197 | ;; alpha[n](z) = G(n+1,z/2)/z^(n+1) |
| 198 | ;; beta[n](z) = G(n+1,z/2)/z^(n+1) - G(n+1,-z/2)/z^(n+1) |
| 199 | |
| 200 | (defun alpha (n z) |
| 201 | (let ((n (float n (realpart z)))) |
| 202 | (/ (cf-incomplete-gamma (1+ n) (/ z 2)) |
| 203 | (expt z (1+ n))))) |
| 204 | |
| 205 | (defun beta (n z) |
| 206 | (let ((n (float n (realpart z)))) |
| 207 | (/ (- (cf-incomplete-gamma (1+ n) (/ z 2)) |
| 208 | (cf-incomplete-gamma (1+ n) (/ z -2))) |
| 209 | (expt z (1+ n))))) |
| 210 | |
| 211 | ;; a[0](k,v) := (k+sqrt(k^2+1))^(-v); |
| 212 | ;; a[1](k,v) := -v*a[0](k,v)/sqrt(k^2+1); |
| 213 | ;; 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)); |
| 214 | |
| 215 | ;; Convert this to iteration instead of using this quick-and-dirty |
| 216 | ;; memoization? |
| 217 | (let ((hash (make-hash-table :test 'equal))) |
| 218 | (defun an-clrhash () |
| 219 | (clrhash hash)) |
| 220 | (defun an-dump-hash () |
| 221 | (maphash #'(lambda (k v) |
| 222 | (format t "~S -> ~S~%" k v)) |
| 223 | hash)) |
| 224 | (defun an (n k v) |
| 225 | (or (gethash (list n k v) hash) |
| 226 | (let ((result |
| 227 | (cond ((= n 0) |
| 228 | (expt (+ k (sqrt (float (1+ (* k k)) (realpart v)))) (- v))) |
| 229 | ((= n 1) |
| 230 | (- (/ (* v (an 0 k v)) |
| 231 | (sqrt (float (1+ (* k k)) (realpart v)))))) |
| 232 | (t |
| 233 | (/ (- (* (- (* v v) (expt (- n 2) 2)) (an (- n 2) k v)) |
| 234 | (* k (- n 1) (+ n n -3) (an (- n 1) k v))) |
| 235 | (+ 1 (* k k)) |
| 236 | (- n 1) |
| 237 | n))))) |
| 238 | (setf (gethash (list n k v) hash) result) |
| 239 | result)))) |
| 240 | |
| 241 | ;; SUM-AN computes the series |
| 242 | ;; |
| 243 | ;; sum(exp(-k*z)*a[n](k,v), k, 1, N) |
| 244 | ;; |
| 245 | (defun sum-an (big-n n v z) |
| 246 | (let ((sum 0)) |
| 247 | (loop for k from 1 upto big-n |
| 248 | do |
| 249 | (incf sum (* (exp (- (* k z))) |
| 250 | (an n k v)))) |
| 251 | sum)) |
| 252 | |
| 253 | ;; SUM-AB computes the series |
| 254 | ;; |
| 255 | ;; sum(alpha[n](z)*a[n](0,v) + beta[n](z)*sum_an(N, n, v, z), n, 0, inf) |
| 256 | (defun sum-ab (big-n v z) |
| 257 | (let ((eps (epsilon (realpart z)))) |
| 258 | (an-clrhash) |
| 259 | (do* ((n 0 (+ 1 n)) |
| 260 | (term (+ (* (alpha n z) (an n 0 v)) |
| 261 | (* (beta n z) (sum-an big-n n v z))) |
| 262 | (+ (* (alpha n z) (an n 0 v)) |
| 263 | (* (beta n z) (sum-an big-n n v z)))) |
| 264 | (sum term (+ sum term))) |
| 265 | ((<= (abs term) (* eps (abs sum))) |
| 266 | sum) |
| 267 | (when nil |
| 268 | (format t "n = ~D~%" n) |
| 269 | (format t " term = ~S~%" term) |
| 270 | (format t " sum = ~S~%" sum))))) |
| 271 | |
| 272 | ;; Convert to iteration instead of this quick-and-dirty memoization? |
| 273 | (let ((hash (make-hash-table :test 'equal))) |
| 274 | (defun %big-a-clrhash () |
| 275 | (clrhash hash)) |
| 276 | (defun %big-a-dump-hash () |
| 277 | (maphash #'(lambda (k v) |
| 278 | (format t "~S -> ~S~%" k v)) |
| 279 | hash)) |
| 280 | (defun %big-a (n v) |
| 281 | (or (gethash (list n v) hash) |
| 282 | (let ((result |
| 283 | (cond ((zerop n) |
| 284 | (expt 2 (- v))) |
| 285 | (t |
| 286 | (* (%big-a (- n 1) v) |
| 287 | (/ (* (+ v n n -2) (+ v n n -1)) |
| 288 | (* 4 n (+ n v)))))))) |
| 289 | (setf (gethash (list n v) hash) result) |
| 290 | result)))) |
| 291 | |
| 292 | ;; Computes A[n](v) = |
| 293 | ;; (-1)^n*v*2^(-v)*pochhammer(v+n+1,n-1)/(2^(2*n)*n!) If v is a |
| 294 | ;; negative integer -m, use A[n](-m) = (-1)^(m+1)*A[n-m](m) for n >= |
| 295 | ;; m. |
| 296 | (defun big-a (n v) |
| 297 | (let ((m (ftruncate v))) |
| 298 | (cond ((and (= m v) (minusp m)) |
| 299 | (if (< n m) |
| 300 | (%big-a n v) |
| 301 | (let ((result (%big-a (+ n m) v))) |
| 302 | (if (oddp (truncate m)) |
| 303 | result |
| 304 | (- result))))) |
| 305 | (t |
| 306 | (%big-a n v))))) |
| 307 | |
| 308 | ;; I[n](t, z, v) = exp(-t*z)/t^(2*n+v-1) * |
| 309 | ;; integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf) |
| 310 | ;; |
| 311 | ;; Use the substitution u=1+s to get a new integral |
| 312 | ;; |
| 313 | ;; integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf) |
| 314 | ;; = exp(t*z) * integrate(u^(-v-2*n)*exp(-t*u*z), u, 1, inf) |
| 315 | ;; = exp(t*z)*t^(v+2*n-1)*z^(v+2*n-1)*incomplete_gamma_tail(1-v-2*n,t*z) |
| 316 | ;; |
| 317 | ;; The continued fraction for incomplete_gamma_tail(a,z) is |
| 318 | ;; |
| 319 | ;; z^a*exp(-z)/CF |
| 320 | ;; |
| 321 | ;; So incomplete_gamma_tail(1-v-2*n, t*z) is |
| 322 | ;; |
| 323 | ;; (t*z)^(1-v-2*n)*exp(-t*z)/CF |
| 324 | ;; |
| 325 | ;; which finally gives |
| 326 | ;; |
| 327 | ;; integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf) |
| 328 | ;; = CF |
| 329 | ;; |
| 330 | ;; and I[n](t, z, v) = exp(-t*z)/t^(2*n+v-1)/CF |
| 331 | (defun big-i (n t z v) |
| 332 | (/ (exp (- (* t z))) |
| 333 | (expt t (+ n n v -1)) |
| 334 | (let* ((a (- 1 v n n)) |
| 335 | (z-a (- z a))) |
| 336 | (lentz #'(lambda (n) |
| 337 | (+ n n 1 z-a)) |
| 338 | #'(lambda (n) |
| 339 | (* n (- a n))))))) |
| 340 | |
| 341 | (defun sum-big-ia (big-n v z) |
| 342 | ) |
| 343 | |
| 344 | (defun bessel-j (v z) |
| 345 | (let ((vv (ftruncate v))) |
| 346 | (cond ((= vv v) |
| 347 | ;; v is an integer |
| 348 | (integer-bessel-j-exp-arc v z)) |
| 349 | (t |
| 350 | (let ((big-n 100) |
| 351 | (vpi (* v (float-pi (realpart z))))) |
| 352 | (+ (integer-bessel-j-exp-arc v z) |
| 353 | (* z |
| 354 | (/ (sin vpi) vpi) |
| 355 | (+ (/ -1 z) |
| 356 | (sum-ab big-n v z))))))))) |