Changeset 64c424


Ignore:
Timestamp:
04/17/12 03:45:12 (3 years ago)
Author:
Raymond Toy <toy.raymond@…>
Branches:
master
Children:
e6a3f2
Parents:
8c5195
Message:

Add iterative versions for some functions.

File:
1 edited

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)))
Note: See TracChangeset for help on using the changeset viewer.