Commit a5a4c7acd95d1946e7cf426f9a927a918c9e2afc

Authored by Raymond Toy
1 parent 1d404aca

Add more parts of the exp-arc algorithm. Needs lots of work, but it

seems that bessel_j(n,z) mostly works.
Showing 1 changed file with 174 additions and 3 deletions
qd-bessel.lisp
... ... @@ -39,12 +39,12 @@
39 39 ;; = 1/2^(k+3/2)/p^(k+1/2)*integrate(t^(k-1/2)*exp(-t),t,0,p)
40 40 ;; = 1/2^(k+3/2)/p^(k+1/2) * g(k+1/2, p)
41 41 ;;
42   -;; where g(a,z) is the lower incomplete gamma function.
  42 +;; where G(a,z) is the lower incomplete gamma function.
43 43 ;;
44   -;; There is the continued fraction expansion for g(a,z) (see
  44 +;; There is the continued fraction expansion for G(a,z) (see
45 45 ;; cf-incomplete-gamma in qd-gamma.lisp):
46 46 ;;
47   -;; g(a,z) = z^a*exp(-z)/ CF
  47 +;; G(a,z) = z^a*exp(-z)/ CF
48 48 ;;
49 49 ;; So
50 50 ;;
... ... @@ -183,6 +183,177 @@
183 183 i-))
184 184 (float-pi i+)
185 185 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)))))))))
186 357  
187 358 (defun paris-series (v z n)
188 359 (labels ((pochhammer (a k)
... ...