Changeset 64c424b7f41413424bbab93d4f5dd4f8461b784d
 Timestamp:
 04/16/12 20:45:12 (2 years ago)
 Author:
 Raymond Toy <toy.raymond@…>
 Children:
 e6a3f286122ea99d65a6448c3c41d0bdc411b2d7
 Parents:
 8c5195a87137fd61952a175a5dae676c70b480ef
 gitcommitter:
 Raymond Toy <toy.raymond@…> (04/16/12 20:45:12)
 Message:

Add iterative versions for some functions.
 Files:

Legend:
 Unmodified
 Added
 Removed

r8c5195

r64c424


78  78  (* (ash n 1) p) 
79  79  ( (* (+ a (ash n 1)) p)))))))) 
 80  
 81  ;; Use the recursion 
 82  (defun bkiter (k p oldbk) 
 83  (withfloatingpointcontagion (p oldbk) 
 84  (if (zerop k) 
 85  (* (sqrt (/ (floatpi p) 8)) 
 86  (let ((rp (sqrt p))) 
 87  (/ (erf rp) 
 88  rp))) 
 89  ( (* ( k 1/2) 
 90  (/ oldbk (* 2 p))) 
 91  (/ (exp ( p)) 
 92  p 
 93  (ash 1 (+ k 1)) 
 94  (sqrt (float 2 (realpart p)))))))) 
80  95  
81  96  ;; exparc I function, as given in the Laguerre paper 
… 
… 

179  194  (format t " sum  ~S~%" sum))))) 
180  195  
 196  (defun exparci3 (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  (bkiter 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 *debugexparc* 
 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 *debugexparc* 
 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, 
… 
… 

213  257  ;; alpha[n](z) =  exp(z/2)/2^n/z + n/z*alpha[n1](z) 
214  258  ;; beta[n]z) = ((1)^n*exp(z/2)exp(z/2))/2^n/z + n/z*beta[n1](z) 
 259  ;; = (1)^n/(2^n)*2*sinh(z/2)/z + n/z*beta[n1](z) 
215  260  ;; 
216  261  ;; We also note that 
… 
… 

223  268  (/ (incompletegamma (1+ n) (/ z 2)) 
224  269  (expt z (1+ n))))) 
 270  
 271  (defun alphaiter (n z alphaold) 
 272  (if (zerop n) 
 273  ;; (1 exp(z/2))/z. 
 274  (/ ( 1 (exp (* z 1/2))) 
 275  z) 
 276  ( (* (/ n z) alphaold) 
 277  (/ (exp ( (* z 1/2))) 
 278  z 
 279  (ash 1 n))))) 
225  280  
226  281  (defun beta (n z) 
… 
… 

229  284  (incompletegamma (1+ n) (/ z 2))) 
230  285  (expt z (1+ n))))) 
 286  
 287  (defun betaiter (n z oldbeta) 
 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 (/ oldbeta z)) 
 295  (* (/ (sinh (* 1/2 z)) (* 1/2 z)) 
 296  (scalefloat (float (if (evenp n) 1 1) (realpart z)) ( n)))))) 
 297  
231  298  
232  299  ;; a[0](k,v) := (k+sqrt(k^2+1))^(v); 
… 
… 

298  365  (+ (* (alpha n z) (an n 0 v)) 
299  366  (* (beta n z) (suman bign 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 sumab2 (bign v z) 
 376  (let ((eps (epsilon (realpart z)))) 
 377  (anclrhash) 
 378  (do* ((n 0 (+ 1 n)) 
 379  (alphan (alphaiter 0 z 0) 
 380  (alphaiter n z alphan)) 
 381  (betan (betaiter 0 z 0) 
 382  (betaiter n z betan)) 
 383  (term (+ (* alphan (an n 0 v)) 
 384  (* betan (suman bign n v z))) 
 385  (+ (* alphan (an n 0 v)) 
 386  (* betan (suman bign n v z)))) 
300  387  (sum term (+ sum term))) 
301  388  ((<= (abs term) (* eps (abs sum))) 