Changeset 64c424b7f41413424bbab93d4f5dd4f8461b784d for qdbessel.lisp
 Timestamp:
 04/16/12 20:45:12 (2 years ago)
 Children:
 e6a3f286122ea99d65a6448c3c41d0bdc411b2d7
 Parents:
 8c5195a87137fd61952a175a5dae676c70b480ef
 gitcommitter:
 Raymond Toy <toy.raymond@…> (04/16/12 20:45:12)
 Files:

 1 modified
Legend:
 Unmodified
 Added
 Removed

qdbessel.lisp
r8c5195 r64c424 78 78 (* (ash n 1) p) 79 79 ( (* (+ a (ash n 1)) p)))))))) 80 81 ;; Use the recursion 82 (defun bkiter (k p oldbk) 83 (withfloatingpointcontagion (p oldbk) 84 (if (zerop k) 85 (* (sqrt (/ (floatpi p) 8)) 86 (let ((rp (sqrt p))) 87 (/ (erf rp) 88 rp))) 89 ( (* ( k 1/2) 90 (/ oldbk (* 2 p))) 91 (/ (exp ( p)) 92 p 93 (ash 1 (+ k 1)) 94 (sqrt (float 2 (realpart p)))))))) 80 95 81 96 ;; exparc I function, as given in the Laguerre paper … … 179 194 (format t " sum  ~S~%" sum))))) 180 195 196 (defun exparci3 (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 (bkiter 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 *debugexparc* 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 *debugexparc* 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 181 225 182 226 ;; Not really just for Bessel J for integer orders, but in that case, … … 213 257 ;; alpha[n](z) =  exp(z/2)/2^n/z + n/z*alpha[n1](z) 214 258 ;; beta[n]z) = ((1)^n*exp(z/2)exp(z/2))/2^n/z + n/z*beta[n1](z) 259 ;; = (1)^n/(2^n)*2*sinh(z/2)/z + n/z*beta[n1](z) 215 260 ;; 216 261 ;; We also note that … … 223 268 (/ (incompletegamma (1+ n) (/ z 2)) 224 269 (expt z (1+ n))))) 270 271 (defun alphaiter (n z alphaold) 272 (if (zerop n) 273 ;; (1 exp(z/2))/z. 274 (/ ( 1 (exp (* z 1/2))) 275 z) 276 ( (* (/ n z) alphaold) 277 (/ (exp ( (* z 1/2))) 278 z 279 (ash 1 n))))) 225 280 226 281 (defun beta (n z) … … 229 284 (incompletegamma (1+ n) (/ z 2))) 230 285 (expt z (1+ n))))) 286 287 (defun betaiter (n z oldbeta) 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 (/ oldbeta z)) 295 (* (/ (sinh (* 1/2 z)) (* 1/2 z)) 296 (scalefloat (float (if (evenp n) 1 1) (realpart z)) ( n)))))) 297 231 298 232 299 ;; a[0](k,v) := (k+sqrt(k^2+1))^(v); … … 298 365 (+ (* (alpha n z) (an n 0 v)) 299 366 (* (beta n z) (suman bign 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 375 (defun sumab2 (bign v z) 376 (let ((eps (epsilon (realpart z)))) 377 (anclrhash) 378 (do* ((n 0 (+ 1 n)) 379 (alphan (alphaiter 0 z 0) 380 (alphaiter n z alphan)) 381 (betan (betaiter 0 z 0) 382 (betaiter n z betan)) 383 (term (+ (* alphan (an n 0 v)) 384 (* betan (suman bign n v z))) 385 (+ (* alphan (an n 0 v)) 386 (* betan (suman bign n v z)))) 300 387 (sum term (+ sum term))) 301 388 ((<= (abs term) (* eps (abs sum)))