Changeset 64c424b7f41413424bbab93d4f5dd4f8461b784d

Show
Ignore:
Timestamp:
04/16/12 20:45:12 (2 years ago)
Author:
Raymond Toy <toy.raymond@…>
Children:
e6a3f286122ea99d65a6448c3c41d0bdc411b2d7
Parents:
8c5195a87137fd61952a175a5dae676c70b480ef
git-committer:
Raymond Toy <toy.raymond@…> (04/16/12 20:45:12)
Message:

Add iterative versions for some functions.

Files:
1 modified

Legend:

Unmodified
Added
Removed
  • qd-bessel.lisp

    r8c5195 r64c424  
    7878                      (* (ash n -1) p) 
    7979                      (- (* (+ a (ash n -1)) p)))))))) 
     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)))))))) 
    8095 
    8196;; exp-arc I function, as given in the Laguerre paper 
     
    179194        (format t " sum   - ~S~%" sum))))) 
    180195 
     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 
    181225 
    182226;; Not really just for Bessel J for integer orders, but in that case, 
     
    213257;; alpha[n](z) = - exp(-z/2)/2^n/z + n/z*alpha[n-1](z) 
    214258;; 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) 
    215260;; 
    216261;; We also note that 
     
    223268    (/ (incomplete-gamma (1+ n) (/ z 2)) 
    224269       (expt z (1+ n))))) 
     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))))) 
    225280 
    226281(defun beta (n z) 
     
    229284          (incomplete-gamma (1+ n) (/ z -2))) 
    230285       (expt z (1+ n))))) 
     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 
    231298 
    232299;; a[0](k,v) := (k+sqrt(k^2+1))^(-v); 
     
    298365                (+ (* (alpha n z) (an n 0 v)) 
    299366                   (* (beta n z) (sum-an big-n 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 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)))) 
    300387          (sum term (+ sum term))) 
    301388         ((<= (abs term) (* eps (abs sum)))