Commit f203389e5e78d1f001f68447ac2a9dd86dcfbbf6

Authored by Raymond Toy
1 parent f8943af6

* Implement sum-big-ia.

* Add series for Bessel J.  (Not working yet.)
Showing 1 changed file with 36 additions and 2 deletions
qd-bessel.lisp
... ... @@ -324,8 +324,41 @@
324 324 (incomplete-gamma-tail a (* theta z)))))
325 325  
326 326 (defun sum-big-ia (big-n v z)
327   - )
  327 + (let ((big-n-1/2 (+ big-n 1/2))
  328 + (eps (epsilon z)))
  329 + (do* ((n 0 (1+ n))
  330 + (term (* (big-a 0 v)
  331 + (big-i 0 big-n-1/2 z v))
  332 + (* (big-a n v)
  333 + (big-i n big-n-1/2 z v)))
  334 + (sum term (+ sum term)))
  335 + ((<= (abs term) (* eps (abs sum)))
  336 + sum)
  337 + #+nil
  338 + (progn
  339 + (format t "n = ~D~%" n)
  340 + (format t " term = ~S~%" term)
  341 + (format t " sum = ~S~%" sum)))))
328 342  
  343 +;; Series for bessel J:
  344 +;;
  345 +;; (z/2)^v*sum((-1)^k/Gamma(k+v+1)/k!*(z^2//4)^k, k, 0, inf)
  346 +(defun s-bessel-j (v z)
  347 + (with-floating-point-contagion (v z)
  348 + (let ((z2/4 (* z z 1/4))
  349 + (eps (epsilon z)))
  350 + (do* ((k 0 (+ 1 k))
  351 + (f (gamma (+ v 1))
  352 + (* f (* k (+ v k))))
  353 + (term (/ f)
  354 + (/ (* (- term) z2/4) f))
  355 + (sum term (+ sum term)))
  356 + ((<= (abs term) (* eps (abs sum)))
  357 + (* sum (expt (* z 1/2) v)))
  358 + (format t "k = ~D~%" k)
  359 + (format t " term = ~S~%" term)
  360 + (format t " sum = ~S~%" sum)))))
  361 +
329 362 (defun bessel-j (v z)
330 363 (let ((vv (ftruncate v)))
331 364 (cond ((= vv v)
... ... @@ -338,7 +371,8 @@
338 371 (* z
339 372 (/ (sin vpi) vpi)
340 373 (+ (/ -1 z)
341   - (sum-ab big-n v z)))))))))
  374 + (sum-ab big-n v z)
  375 + (sum-big-ia big-n v z)))))))))
342 376  
343 377 (defun paris-series (v z n)
344 378 (labels ((pochhammer (a k)
... ...