Changeset a5a4c7acd95d1946e7cf426f9a927a918c9e2afc

Show
Ignore:
Timestamp:
04/07/12 09:28:16 (2 years ago)
Author:
Raymond Toy <toy.raymond@…>
Children:
015558a3df8fae2686c7b6e2f049da3f4b1a2cd8
Parents:
1d404aca1e6652ad2456ba14e8bbdf5d329c06a0
git-committer:
Raymond Toy <toy.raymond@…> (04/07/12 09:28:16)
Message:

Add more parts of the exp-arc algorithm. Needs lots of work, but it
seems that bessel_j(n,z) mostly works.

Files:
1 modified

Legend:

Unmodified
Added
Removed
  • qd-bessel.lisp

    r1d404a ra5a4c7  
    4040;;         = 1/2^(k+3/2)/p^(k+1/2) * g(k+1/2, p) 
    4141;; 
    42 ;; where g(a,z) is the lower incomplete gamma function. 
    43 ;; 
    44 ;; There is the continued fraction expansion for g(a,z) (see 
     42;; where G(a,z) is the lower incomplete gamma function. 
     43;; 
     44;; There is the continued fraction expansion for G(a,z) (see 
    4545;; cf-incomplete-gamma in qd-gamma.lisp): 
    4646;; 
    47 ;;  g(a,z) = z^a*exp(-z)/ CF 
     47;;  G(a,z) = z^a*exp(-z)/ CF 
    4848;; 
    4949;; So 
     
    184184       (float-pi i+) 
    185185       2))) 
     186 
     187;; alpha[n](z) = integrate(exp(-z*s)*s^n, s, 0, 1/2) 
     188;; beta[n](z)  = integrate(exp(-z*s)*s^n, s, -1/2, 1/2) 
     189;; 
     190;; The recurrence in [2] is 
     191;; 
     192;; alpha[n](z) = - exp(-z/2)/2^n/z + n/z*alpha[n-1](z) 
     193;; beta[n]z)   = ((-1)^n*exp(z/2)-exp(-z/2))/2^n/z + n/z*beta[n-1](z) 
     194;; 
     195;; We also note that 
     196;; 
     197;; alpha[n](z) = G(n+1,z/2)/z^(n+1) 
     198;; beta[n](z)  = G(n+1,z/2)/z^(n+1) - G(n+1,-z/2)/z^(n+1) 
     199 
     200(defun alpha (n z) 
     201  (let ((n (float n (realpart z)))) 
     202    (/ (cf-incomplete-gamma (1+ n) (/ z 2)) 
     203       (expt z (1+ n))))) 
     204 
     205(defun beta (n z) 
     206  (let ((n (float n (realpart z)))) 
     207    (/ (- (cf-incomplete-gamma (1+ n) (/ z 2)) 
     208          (cf-incomplete-gamma (1+ n) (/ z -2))) 
     209       (expt z (1+ n))))) 
     210 
     211;; a[0](k,v) := (k+sqrt(k^2+1))^(-v); 
     212;; a[1](k,v) := -v*a[0](k,v)/sqrt(k^2+1); 
     213;; a[n](k,v) := 1/(k^2+1)/(n-1)/n*((v^2-(n-2)^2)*a[n-2](k,v)-k*(n-1)*(2*n-3)*a[n-1](k,v)); 
     214 
     215;; Convert this to iteration instead of using this quick-and-dirty 
     216;; memoization? 
     217(let ((hash (make-hash-table :test 'equal))) 
     218  (defun an-clrhash () 
     219    (clrhash hash)) 
     220  (defun an-dump-hash () 
     221    (maphash #'(lambda (k v) 
     222                 (format t "~S -> ~S~%" k v)) 
     223             hash)) 
     224  (defun an (n k v) 
     225    (or (gethash (list n k v) hash) 
     226        (let ((result 
     227                (cond ((= n 0) 
     228                       (expt (+ k (sqrt (float (1+ (* k k)) (realpart v)))) (- v))) 
     229                      ((= n 1) 
     230                       (- (/ (* v (an 0 k v)) 
     231                             (sqrt (float (1+ (* k k)) (realpart v)))))) 
     232                      (t 
     233                       (/ (- (* (- (* v v) (expt (- n 2) 2)) (an (- n 2) k v)) 
     234                             (* k (- n 1) (+ n n -3) (an (- n 1) k v))) 
     235                          (+ 1 (* k k)) 
     236                          (- n 1) 
     237                          n))))) 
     238          (setf (gethash (list n k v) hash) result) 
     239          result)))) 
     240 
     241;; SUM-AN computes the series 
     242;; 
     243;; sum(exp(-k*z)*a[n](k,v), k, 1, N) 
     244;; 
     245(defun sum-an (big-n n v z) 
     246  (let ((sum 0)) 
     247    (loop for k from 1 upto big-n 
     248          do 
     249             (incf sum (* (exp (- (* k z))) 
     250                          (an n k v)))) 
     251    sum)) 
     252 
     253;; SUM-AB computes the series 
     254;; 
     255;; sum(alpha[n](z)*a[n](0,v) + beta[n](z)*sum_an(N, n, v, z), n, 0, inf) 
     256(defun sum-ab (big-n v z) 
     257  (let ((eps (epsilon (realpart z)))) 
     258    (an-clrhash) 
     259    (do* ((n 0 (+ 1 n)) 
     260          (term (+ (* (alpha n z) (an n 0 v)) 
     261                   (* (beta n z) (sum-an big-n n v z))) 
     262                (+ (* (alpha n z) (an n 0 v)) 
     263                   (* (beta n z) (sum-an big-n n v z)))) 
     264          (sum term (+ sum term))) 
     265         ((<= (abs term) (* eps (abs sum))) 
     266          sum) 
     267      (when nil 
     268        (format t "n = ~D~%" n) 
     269        (format t " term = ~S~%" term) 
     270        (format t " sum  = ~S~%" sum))))) 
     271 
     272;; Convert to iteration instead of this quick-and-dirty memoization? 
     273(let ((hash (make-hash-table :test 'equal))) 
     274  (defun %big-a-clrhash () 
     275    (clrhash hash)) 
     276  (defun %big-a-dump-hash () 
     277    (maphash #'(lambda (k v) 
     278                 (format t "~S -> ~S~%" k v)) 
     279             hash)) 
     280  (defun %big-a (n v) 
     281    (or (gethash (list n v) hash) 
     282        (let ((result 
     283                (cond ((zerop n) 
     284                       (expt 2 (- v))) 
     285                      (t 
     286                       (* (%big-a (- n 1) v) 
     287                          (/ (* (+ v n n -2) (+ v n n -1)) 
     288                             (* 4 n (+ n v)))))))) 
     289          (setf (gethash (list n v) hash) result) 
     290          result)))) 
     291 
     292;; Computes A[n](v) = 
     293;; (-1)^n*v*2^(-v)*pochhammer(v+n+1,n-1)/(2^(2*n)*n!)  If v is a 
     294;; negative integer -m, use A[n](-m) = (-1)^(m+1)*A[n-m](m) for n >= 
     295;; m. 
     296(defun big-a (n v) 
     297  (let ((m (ftruncate v))) 
     298    (cond ((and (= m v) (minusp m)) 
     299           (if (< n m) 
     300               (%big-a n v) 
     301               (let ((result (%big-a (+ n m) v))) 
     302                 (if (oddp (truncate m)) 
     303                     result 
     304                     (- result))))) 
     305          (t 
     306           (%big-a n v))))) 
     307 
     308;; I[n](t, z, v) = exp(-t*z)/t^(2*n+v-1) * 
     309;;                  integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf) 
     310;; 
     311;; Use the substitution u=1+s to get a new integral 
     312;; 
     313;; integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf) 
     314;;   = exp(t*z) * integrate(u^(-v-2*n)*exp(-t*u*z), u, 1, inf) 
     315;;   = exp(t*z)*t^(v+2*n-1)*z^(v+2*n-1)*incomplete_gamma_tail(1-v-2*n,t*z) 
     316;; 
     317;; The continued fraction for incomplete_gamma_tail(a,z) is 
     318;; 
     319;;   z^a*exp(-z)/CF 
     320;; 
     321;; So incomplete_gamma_tail(1-v-2*n, t*z) is 
     322;; 
     323;;   (t*z)^(1-v-2*n)*exp(-t*z)/CF 
     324;; 
     325;; which finally gives 
     326;; 
     327;; integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf) 
     328;;   = CF 
     329;; 
     330;; and I[n](t, z, v) = exp(-t*z)/t^(2*n+v-1)/CF 
     331(defun big-i (n t z v) 
     332  (/ (exp (- (* t z))) 
     333     (expt t (+ n n v -1)) 
     334     (let* ((a (- 1 v n n)) 
     335            (z-a (- z a))) 
     336       (lentz #'(lambda (n) 
     337                  (+ n n 1 z-a)) 
     338              #'(lambda (n) 
     339                  (* n (- a n))))))) 
     340 
     341(defun sum-big-ia (big-n v z) 
     342  ) 
     343 
     344(defun bessel-j (v z) 
     345  (let ((vv (ftruncate v))) 
     346    (cond ((= vv v) 
     347           ;; v is an integer 
     348           (integer-bessel-j-exp-arc v z)) 
     349          (t 
     350           (let ((big-n 100) 
     351                 (vpi (* v (float-pi (realpart z))))) 
     352             (+ (integer-bessel-j-exp-arc v z) 
     353                (* z 
     354                   (/ (sin vpi) vpi) 
     355                   (+ (/ -1 z) 
     356                      (sum-ab big-n v z))))))))) 
    186357 
    187358(defun paris-series (v z n)