Commit 9fd2ebcbeed3ecae899f732b15fd279f6fb0f14f

Authored by Raymond Toy
Committed by Raymond Toy
1 parent 6ab5226c

* Use the g[k] formula instead of r[2*k+1] because we fail to

converge when v = 0.
 * Clear hash tables in bessel-j.
Showing 1 changed file with 6 additions and 2 deletions
qd-bessel.lisp
... ... @@ -156,7 +156,8 @@
156 156 (do* ((k 0 (1+ k))
157 157 (bk (bk 0 p)
158 158 (bk k p))
159   - (ratio v
  159 + ;; Compute g[k](p)/(2*k)!, not r[2*k+1](p)/(2*k)!
  160 + (ratio 1
160 161 (* ratio (/ (+ v2 (expt (1- (* 2 k)) 2))
161 162 (* 2 k (1- (* 2 k))))))
162 163 (term (* ratio bk)
... ... @@ -169,7 +170,7 @@
169 170 (format t " ratio = ~S~%" ratio)
170 171 (format t " term = ~S~%" term)
171 172 (format t " sum - ~S~%" sum))
172   - (* sum #c(0 2) (/ (exp p) q)))
  173 + (* sum 4 (exp p)))
173 174 (when *debug-exparc*
174 175 (format t "k = ~D~%" k)
175 176 (format t " bk = ~S~%" bk)
... ... @@ -390,6 +391,9 @@
390 391 ;;
391 392 (defun bessel-j (v z)
392 393 (let ((vv (ftruncate v)))
  394 + ;; Clear the caches for now.
  395 + (an-clrhash)
  396 + (%big-a-clrhash)
393 397 (cond ((= vv v)
394 398 ;; v is an integer
395 399 (integer-bessel-j-exp-arc v z))
... ...