Commit 64c424b7f41413424bbab93d4f5dd4f8461b784d

Authored by Raymond Toy
1 parent 8c5195a8

Add iterative versions for some functions.

Showing 1 changed file with 87 additions and 0 deletions
qd-bessel.lisp
... ... @@ -78,6 +78,21 @@
78 78 (* (ash n -1) p)
79 79 (- (* (+ a (ash n -1)) p))))))))
80 80  
  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 +
81 96 ;; exp-arc I function, as given in the Laguerre paper
82 97 ;;
83 98 ;; I(p, q) = 4*exp(p) * sum(g[k](-2*%i*q)/(2*k)!*B[k](p), k, 0, inf)
... ... @@ -178,6 +193,35 @@
178 193 (format t " term = ~S~%" term)
179 194 (format t " sum - ~S~%" sum)))))
180 195  
  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 +
181 225  
182 226 ;; Not really just for Bessel J for integer orders, but in that case,
183 227 ;; this is all that's needed to compute Bessel J. For other values,
... ... @@ -212,6 +256,7 @@
212 256 ;;
213 257 ;; alpha[n](z) = - exp(-z/2)/2^n/z + n/z*alpha[n-1](z)
214 258 ;; beta[n]z) = ((-1)^n*exp(z/2)-exp(-z/2))/2^n/z + n/z*beta[n-1](z)
  259 +;; = (-1)^n/(2^n)*2*sinh(z/2)/z + n/z*beta[n-1](z)
215 260 ;;
216 261 ;; We also note that
217 262 ;;
... ... @@ -223,12 +268,34 @@
223 268 (/ (incomplete-gamma (1+ n) (/ z 2))
224 269 (expt z (1+ n)))))
225 270  
  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 +
226 281 (defun beta (n z)
227 282 (let ((n (float n (realpart z))))
228 283 (/ (- (incomplete-gamma (1+ n) (/ z 2))
229 284 (incomplete-gamma (1+ n) (/ z -2)))
230 285 (expt z (1+ n)))))
231 286  
  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 +
232 299 ;; a[0](k,v) := (k+sqrt(k^2+1))^(-v);
233 300 ;; a[1](k,v) := -v*a[0](k,v)/sqrt(k^2+1);
234 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));
... ... @@ -305,6 +372,26 @@
305 372 (format t " term = ~S~%" term)
306 373 (format t " sum = ~S~%" sum)))))
307 374  
  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 +
308 395 ;; Convert to iteration instead of this quick-and-dirty memoization?
309 396 (let ((hash (make-hash-table :test 'equal)))
310 397 (defun %big-a-clrhash ()
... ...