This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "OCT: A portable Lisp implementation for quad-double precision floats". The branch, master has been updated via 6cfb0ac4b6bcc1a25bc119e87fd2b57bfa1f4355 (commit) via 8ec0d2004ed4063fd568250b00d940b208573a92 (commit) via bd6d814332fdf6a831873bb33c9e4753f29d414f (commit) via c7fc989dad090ada2c046b28bf43406810474a77 (commit) via 95aa580ce0dd62113535c30c4c0c82f8656c2d86 (commit) via dbc1e376add1d492dc35c37b3098c2ffb796d6f1 (commit) via 6cd96249b94dc29962cecda996a5799d8ac27271 (commit) via c86217a5f7191f1a29566d6ee24cb57344dd546d (commit) via 015558a3df8fae2686c7b6e2f049da3f4b1a2cd8 (commit) via a5a4c7acd95d1946e7cf426f9a927a918c9e2afc (commit) via 1d404aca1e6652ad2456ba14e8bbdf5d329c06a0 (commit) via 4c1ed0f45e79bbe467f7d4694f417ff92ae46adb (commit) from dd5c2a4a6628648c5fe9a7eec98c82a9d284daa3 (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit 6cfb0ac4b6bcc1a25bc119e87fd2b57bfa1f4355 Merge: 8ec0d20 dd5c2a4 Author: Raymond Toy <toy.raymond@gmail.com> Date: Mon Apr 9 08:37:51 2012 -0700 Merge branch 'master' of ssh://common-lisp.net/var/git/projects/oct/oct diff --cc qd-gamma.lisp index 410c315,f0e2418..fe79813 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@@ -378,10 -372,13 +378,11 @@@ "Tail of the incomplete gamma function defined by: integrate(t^(a-1)*exp(-t), t, z, inf)" - (let* ((prec (float-contagion a z)) - (a (apply-contagion a prec)) - (z (apply-contagion z prec))) + (with-floating-point-contagion (a z) - (if (zerop a) - ;; incomplete_gamma_tail(0, z) = exp_integral_e(1,z) - (exp-integral-e 1 z) + (if (and (realp a) (<= a 0)) + ;; incomplete_gamma_tail(v, z) = z^v*exp_integral_e(1-a,z) + (* (expt z a) + (exp-integral-e (- 1 a) z)) (if (and (zerop (imagpart a)) (zerop (imagpart z))) ;; For real values, we split the result to compute either the @@@ -531,20 -528,30 +532,36 @@@ ;; for |arg(z)| < pi. ;; ;; - (cond ((and (realp v) (minusp v)) - ;; E(-v, z) = z^(-v-1)*incomplete_gamma_tail(v+1,z) - (let ((-v (- v))) + (let* ((prec (float-contagion v z)) + (v (apply-contagion v prec)) + (z (apply-contagion z prec))) + (cond ((and (realp v) (minusp v)) + ;; E(-v, z) = z^(-v-1)*incomplete_gamma_tail(v+1,z) + (let ((-v (- v))) + (* (expt z (- v 1)) + (incomplete-gamma-tail (+ -v 1) z)))) + ((< (abs z) 1) + ;; Use series for small z + (s-exp-integral-e v z)) + ((>= (abs (phase z)) 3.1) + ;; The continued fraction doesn't converge on the negative + ;; real axis, and converges very slowly near the negative + ;; real axis, so use the incomplete-gamma-tail function in + ;; this region. "Closeness" to the negative real axis is + ;; teken to mean that z is in a sector near the axis. + ;; + ;; E(v,z) = z^(v-1)*incomplete_gamma_tail(1-v,z) (* (expt z (- v 1)) - (incomplete-gamma-tail (- 1 v) z))) - (t - ;; Use continued fraction for everything else. - (cf-exp-integral-e v z))))) + (incomplete-gamma-tail (+ -v 1) z)))) + ((or (< (abs z) 1) (>= (abs (phase z)) 3.1)) + ;; Use series for small z or if z is near the negative real + ;; axis because the continued fraction does not converge on + ;; the negative axis and converges slowly near the negative + ;; axis. + (s-exp-integral-e v z)) + (t + ;; Use continued fraction for everything else. + (cf-exp-integral-e v z)))) ;; Series for Fresnel S ;; diff --cc qd-methods.lisp index f73c02b,4fac522..2a38e58 --- a/qd-methods.lisp +++ b/qd-methods.lisp @@@ -1139,4 -1132,4 +1146,4 @@@ the same precision as the argument. Th (pprint-logical-block (stream nil :per-line-prefix " ") (apply #'format stream (simple-condition-format-control condition) -- (simple-condition-format-arguments condition)))))) ++ (simple-condition-format-arguments condition)))))) commit 8ec0d2004ed4063fd568250b00d940b208573a92 Author: Raymond Toy <toy.raymond@gmail.com> Date: Sun Apr 8 10:14:37 2012 -0700 Define FLOATP, fix bugs in FLOAT. qd-methods.lisp: * Define FLOATP * Fix bugs in FLOAT: * (FLOAT float nil) is an error * (FLOAT float) returns the float * (FLOAT rational) returns a single-float. qd-package.lisp: o Export FLOATP, shadowing CL:FLOAT. rt-tests.lisp: o Add a few tests for FLOAT. diff --git a/qd-methods.lisp b/qd-methods.lisp index b54fb4b..f73c02b 100644 --- a/qd-methods.lisp +++ b/qd-methods.lisp @@ -255,6 +255,9 @@ (make-instance 'qd-real :value (rational-to-qd bignum))) +(defun floatp (x) + (typep x '(or short-float single-float double-float long-float qd-real))) + (defmethod qfloat ((x real) (num-type cl:float)) (cl:float x num-type)) @@ -299,8 +302,12 @@ x) (declaim (inline float)) -(defun float (x &optional num-type) - (qfloat x (or num-type 1.0))) +(defun float (x &optional (other nil otherp)) + (if otherp + (qfloat x other) + (if (floatp x) + x + (qfloat x 1.0)))) (defmethod qrealpart ((x number)) (cl:realpart x)) diff --git a/qd-package.lisp b/qd-package.lisp index 7db5cd3..ed71347 100644 --- a/qd-package.lisp +++ b/qd-package.lisp @@ -180,6 +180,7 @@ #:decode-float #:scale-float #:float + #:floatp #:floor #:ffloor #:ceiling @@ -252,6 +253,7 @@ #:decode-float #:scale-float #:float + #:floatp #:floor #:ffloor #:ceiling diff --git a/rt-tests.lisp b/rt-tests.lisp index 173c528..eda6b0e 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -55,6 +55,22 @@ ;;; Some simple tests from the Yozo Hida's qd package. +(rt:deftest float.1 + (float 3/2) + 1.5) + +(rt:deftest float.2 + (float 3/2 1d0) + 1.5d0) + +(rt:deftest float.3 + (float 1.5d0) + 1.5d0) + +(rt:deftest float.4 + (= (float #q1.5) #q1.5) + t) + (rt:deftest ceiling-d.1 (multiple-value-list (ceiling -50d0)) (-50 0d0)) commit bd6d814332fdf6a831873bb33c9e4753f29d414f Author: Raymond Toy <toy.raymond@gmail.com> Date: Sun Apr 8 09:59:39 2012 -0700 Oops. The second arg to FLOAT is optional! diff --git a/qd-methods.lisp b/qd-methods.lisp index c461ba3..b54fb4b 100644 --- a/qd-methods.lisp +++ b/qd-methods.lisp @@ -299,8 +299,8 @@ x) (declaim (inline float)) -(defun float (x num-type) - (qfloat x num-type)) +(defun float (x &optional num-type) + (qfloat x (or num-type 1.0))) (defmethod qrealpart ((x number)) (cl:realpart x)) commit c7fc989dad090ada2c046b28bf43406810474a77 Author: Raymond Toy <toy.raymond@gmail.com> Date: Sun Apr 8 09:57:12 2012 -0700 Fix typo in %big-a. Just use incomplete-gamma-tail in big-i. diff --git a/qd-bessel.lisp b/qd-bessel.lisp index ee5b9ef..fbe9760 100644 --- a/qd-bessel.lisp +++ b/qd-bessel.lisp @@ -298,7 +298,7 @@ (cond ((and (= m v) (minusp m)) (if (< n m) (%big-a n v) - (let ((result (%big-a (+ n m) v))) + (let ((result (%big-a (+ n m) (- v)))) (if (oddp (truncate m)) result (- result))))) @@ -310,32 +310,18 @@ ;; ;; Use the substitution u=1+s to get a new integral ;; -;; integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf) -;; = exp(t*z) * integrate(u^(-v-2*n)*exp(-t*u*z), u, 1, inf) -;; = exp(t*z)*t^(v+2*n-1)*z^(v+2*n-1)*incomplete_gamma_tail(1-v-2*n,t*z) -;; -;; The continued fraction for incomplete_gamma_tail(a,z) is -;; -;; z^a*exp(-z)/CF -;; -;; So incomplete_gamma_tail(1-v-2*n, t*z) is -;; -;; (t*z)^(1-v-2*n)/CF +;; integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf) +;; = exp(t*z) * integrate(u^(-v-2*n)*exp(-t*u*z), u, 1, inf) +;; = exp(t*z)*t^(v+2*n-1)*z^(v+2*n-1)*incomplete_gamma_tail(1-v-2*n,t*z) ;; -;; which finally gives +;; Thus, ;; -;; integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf) -;; = (t*z)^(1-v-2*n)/CF +;; I[n](t, z, v) = z^(v+2*n-1)*incomplete_gamma_tail(1-v-2*n,t*z) ;; -;; and I[n](t, z, v) = exp(-t*z)/CF (defun big-i (n theta z v) - (/ (exp (- (* theta z))) - (let* ((a (- 1 v n n)) - (z-a (- (* theta z) a))) - (lentz #'(lambda (n) - (+ n n 1 z-a)) - #'(lambda (n) - (* n (- a n))))))) + (let* ((a (- 1 v n n))) + (* (expt z (- a)) + (incomplete-gamma-tail a (* theta z))))) (defun sum-big-ia (big-n v z) ) commit 95aa580ce0dd62113535c30c4c0c82f8656c2d86 Author: Raymond Toy <toy.raymond@gmail.com> Date: Sun Apr 8 09:37:22 2012 -0700 T is not a variable. Use a different name. diff --git a/qd-bessel.lisp b/qd-bessel.lisp index 6d153a7..ee5b9ef 100644 --- a/qd-bessel.lisp +++ b/qd-bessel.lisp @@ -320,19 +320,18 @@ ;; ;; So incomplete_gamma_tail(1-v-2*n, t*z) is ;; -;; (t*z)^(1-v-2*n)*exp(-t*z)/CF +;; (t*z)^(1-v-2*n)/CF ;; ;; which finally gives ;; -;; integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf) -;; = CF +;; integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf) +;; = (t*z)^(1-v-2*n)/CF ;; -;; and I[n](t, z, v) = exp(-t*z)/t^(2*n+v-1)/CF -(defun big-i (n t z v) - (/ (exp (- (* t z))) - (expt t (+ n n v -1)) +;; and I[n](t, z, v) = exp(-t*z)/CF +(defun big-i (n theta z v) + (/ (exp (- (* theta z))) (let* ((a (- 1 v n n)) - (z-a (- z a))) + (z-a (- (* theta z) a))) (lentz #'(lambda (n) (+ n n 1 z-a)) #'(lambda (n) commit dbc1e376add1d492dc35c37b3098c2ffb796d6f1 Author: Raymond Toy <toy.raymond@gmail.com> Date: Sun Apr 8 09:36:46 2012 -0700 Build qd-bessel now. (qd-bessel is a work in progress.) diff --git a/oct.asd b/oct.asd index 3f8d70c..0ca30e5 100644 --- a/oct.asd +++ b/oct.asd @@ -67,7 +67,8 @@ :depends-on ("qd-methods" "qd-reader")) (:file "qd-gamma" :depends-on ("qd-methods" "qd-reader")) - )) + (:file "qd-bessel" + :depends-on ("qd-methods")))) (defmethod perform ((op test-op) (c (eql (asdf:find-system :oct)))) (oos 'test-op 'oct-tests)) commit 6cd96249b94dc29962cecda996a5799d8ac27271 Author: Raymond Toy <toy.raymond@gmail.com> Date: Sun Apr 8 09:35:35 2012 -0700 Define macro WITH-FLOATING-POINT-CONTAGION. qd-methods.lisp: o Define the macro qd-gamma.lisp: o Use it. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index eb8c55b..410c315 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -108,10 +108,11 @@ ;; log(gamma(-z)) = log(pi)-log(-z)-log(sin(pi*z))-log(gamma(z)) ;; Or ;; log(gamma(z)) = log(pi)-log(-z)-log(sin(pi*z))-log(gamma(-z)) - (- (apply-contagion (log pi) precision) - (log (- z)) - (apply-contagion (log (sin (* pi z))) precision) - (log-gamma (- z)))) + (let ((p (float-pi z))) + (- (log p) + (log (- z)) + (log (sin (* p z))) + (log-gamma (- z))))) (t (let ((absz (abs z))) (cond ((>= absz limit) @@ -377,9 +378,7 @@ "Tail of the incomplete gamma function defined by: integrate(t^(a-1)*exp(-t), t, z, inf)" - (let* ((prec (float-contagion a z)) - (a (apply-contagion a prec)) - (z (apply-contagion z prec))) + (with-floating-point-contagion (a z) (if (zerop a) ;; incomplete_gamma_tail(0, z) = exp_integral_e(1,z) (exp-integral-e 1 z) @@ -409,9 +408,7 @@ "Incomplete gamma function defined by: integrate(t^(a-1)*exp(-t), t, 0, z)" - (let* ((prec (float-contagion a z)) - (a (apply-contagion a prec)) - (z (apply-contagion z prec))) + (with-floating-point-contagion (a z) (if (and (< (abs a) 1) (< (abs z) 1)) (s-incomplete-gamma a z) (if (and (realp a) (realp z)) @@ -539,19 +536,12 @@ (let ((-v (- v))) (* (expt z (- v 1)) (incomplete-gamma-tail (+ -v 1) z)))) - ((< (abs z) 1) - ;; Use series for small z + ((or (< (abs z) 1) (>= (abs (phase z)) 3.1)) + ;; Use series for small z or if z is near the negative real + ;; axis because the continued fraction does not converge on + ;; the negative axis and converges slowly near the negative + ;; axis. (s-exp-integral-e v z)) - ((>= (abs (phase z)) 3.1) - ;; The continued fraction doesn't converge on the negative - ;; real axis, and converges very slowly near the negative - ;; real axis, so use the incomplete-gamma-tail function in - ;; this region. "Closeness" to the negative real axis is - ;; teken to mean that z is in a sector near the axis. - ;; - ;; E(v,z) = z^(v-1)*incomplete_gamma_tail(1-v,z) - (* (expt z (- v 1)) - (incomplete-gamma-tail (- 1 v) z))) (t ;; Use continued fraction for everything else. (cf-exp-integral-e v z)))) diff --git a/qd-methods.lisp b/qd-methods.lisp index 0eef82d..c461ba3 100644 --- a/qd-methods.lisp +++ b/qd-methods.lisp @@ -79,6 +79,17 @@ (complex (coerce (realpart number) precision) (coerce (imagpart number) precision))))) +;; WITH-FLOATING-POINT-CONTAGION - macro +;; +;; Determines the highest precision of the variables in VARLIST and +;; converts each of the values to that precision. +(defmacro with-floating-point-contagion (varlist &body body) + (let ((precision (gensym "PRECISION-"))) + `(let ((,precision (float-contagion ,@varlist))) + (let (,@(mapcar #'(lambda (v) + `(,v (apply-contagion ,v ,precision))) + varlist)) + ,@body)))) (defmethod add1 ((a number)) (cl::1+ a)) commit c86217a5f7191f1a29566d6ee24cb57344dd546d Author: Raymond Toy <toy.raymond@gmail.com> Date: Sun Apr 8 09:04:18 2012 -0700 Fix some comments. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index be832bc..eb8c55b 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -249,13 +249,18 @@ :format-control "~<Continued fraction failed to converge after ~D iterations.~% Delta = ~S~>" :format-arguments (list *max-cf-iterations* (/ c d)))))))) -;; Continued fraction for erf(b): +;; Continued fraction for erf(z): ;; -;; z[n] = 1+2*n-2*z^2 -;; a[n] = 4*n*z^2 +;; erf(z) = 2*z/sqrt(pi)*exp(-z^2)/K +;; +;; where K is the continued fraction with +;; +;; b[n] = 1+2*n-2*z^2 +;; a[n] = 4*n*z^2 ;; ;; This works ok, but has problems for z > 3 where sometimes the -;; result is greater than 1. +;; result is greater than 1 and for larger values, negative numbers +;; are returned. #+nil (defun cf-erf (z) (let* ((z2 (* z z)) commit 015558a3df8fae2686c7b6e2f049da3f4b1a2cd8 Author: Raymond Toy <toy.raymond@gmail.com> Date: Sat Apr 7 23:28:08 2012 -0700 Fix a divide by zero error in s-exp-integral-e for v = 1. We need to skip the first term in the series. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index e353ff5..be832bc 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -257,7 +257,7 @@ ;; This works ok, but has problems for z > 3 where sometimes the ;; result is greater than 1. #+nil -(defun erf (z) +(defun cf-erf (z) (let* ((z2 (* z z)) (twoz2 (* 2 z2))) (* (/ (* 2 z) @@ -504,10 +504,13 @@ (- (psi v) (log z))) (loop for k from 0 for term = 1 then (* term (/ -z k)) - for sum = (/ (- 1 v)) then (+ sum (let ((denom (- k n-1))) - (if (zerop denom) - 0 - (/ term denom)))) + for sum = (if (zerop n-1) + 0 + (/ (- 1 v))) + then (+ sum (let ((denom (- k n-1))) + (if (zerop denom) + 0 + (/ term denom)))) when (< (abs term) (* (abs sum) eps)) return sum))) (loop for k from 0 commit a5a4c7acd95d1946e7cf426f9a927a918c9e2afc Author: Raymond Toy <toy.raymond@gmail.com> Date: Sat Apr 7 09:28:16 2012 -0700 Add more parts of the exp-arc algorithm. Needs lots of work, but it seems that bessel_j(n,z) mostly works. diff --git a/qd-bessel.lisp b/qd-bessel.lisp index e85d2d1..6d153a7 100644 --- a/qd-bessel.lisp +++ b/qd-bessel.lisp @@ -39,12 +39,12 @@ ;; = 1/2^(k+3/2)/p^(k+1/2)*integrate(t^(k-1/2)*exp(-t),t,0,p) ;; = 1/2^(k+3/2)/p^(k+1/2) * g(k+1/2, p) ;; -;; where g(a,z) is the lower incomplete gamma function. +;; where G(a,z) is the lower incomplete gamma function. ;; -;; There is the continued fraction expansion for g(a,z) (see +;; There is the continued fraction expansion for G(a,z) (see ;; cf-incomplete-gamma in qd-gamma.lisp): ;; -;; g(a,z) = z^a*exp(-z)/ CF +;; G(a,z) = z^a*exp(-z)/ CF ;; ;; So ;; @@ -183,6 +183,177 @@ i-)) (float-pi i+) 2))) + +;; alpha[n](z) = integrate(exp(-z*s)*s^n, s, 0, 1/2) +;; beta[n](z) = integrate(exp(-z*s)*s^n, s, -1/2, 1/2) +;; +;; The recurrence in [2] is +;; +;; alpha[n](z) = - exp(-z/2)/2^n/z + n/z*alpha[n-1](z) +;; beta[n]z) = ((-1)^n*exp(z/2)-exp(-z/2))/2^n/z + n/z*beta[n-1](z) +;; +;; We also note that +;; +;; alpha[n](z) = G(n+1,z/2)/z^(n+1) +;; beta[n](z) = G(n+1,z/2)/z^(n+1) - G(n+1,-z/2)/z^(n+1) + +(defun alpha (n z) + (let ((n (float n (realpart z)))) + (/ (cf-incomplete-gamma (1+ n) (/ z 2)) + (expt z (1+ n))))) + +(defun beta (n z) + (let ((n (float n (realpart z)))) + (/ (- (cf-incomplete-gamma (1+ n) (/ z 2)) + (cf-incomplete-gamma (1+ n) (/ z -2))) + (expt z (1+ n))))) + +;; a[0](k,v) := (k+sqrt(k^2+1))^(-v); +;; a[1](k,v) := -v*a[0](k,v)/sqrt(k^2+1); +;; 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)); + +;; Convert this to iteration instead of using this quick-and-dirty +;; memoization? +(let ((hash (make-hash-table :test 'equal))) + (defun an-clrhash () + (clrhash hash)) + (defun an-dump-hash () + (maphash #'(lambda (k v) + (format t "~S -> ~S~%" k v)) + hash)) + (defun an (n k v) + (or (gethash (list n k v) hash) + (let ((result + (cond ((= n 0) + (expt (+ k (sqrt (float (1+ (* k k)) (realpart v)))) (- v))) + ((= n 1) + (- (/ (* v (an 0 k v)) + (sqrt (float (1+ (* k k)) (realpart v)))))) + (t + (/ (- (* (- (* v v) (expt (- n 2) 2)) (an (- n 2) k v)) + (* k (- n 1) (+ n n -3) (an (- n 1) k v))) + (+ 1 (* k k)) + (- n 1) + n))))) + (setf (gethash (list n k v) hash) result) + result)))) + +;; SUM-AN computes the series +;; +;; sum(exp(-k*z)*a[n](k,v), k, 1, N) +;; +(defun sum-an (big-n n v z) + (let ((sum 0)) + (loop for k from 1 upto big-n + do + (incf sum (* (exp (- (* k z))) + (an n k v)))) + sum)) + +;; SUM-AB computes the series +;; +;; sum(alpha[n](z)*a[n](0,v) + beta[n](z)*sum_an(N, n, v, z), n, 0, inf) +(defun sum-ab (big-n v z) + (let ((eps (epsilon (realpart z)))) + (an-clrhash) + (do* ((n 0 (+ 1 n)) + (term (+ (* (alpha n z) (an n 0 v)) + (* (beta n z) (sum-an big-n n v z))) + (+ (* (alpha n z) (an n 0 v)) + (* (beta n z) (sum-an big-n n v z)))) + (sum term (+ sum term))) + ((<= (abs term) (* eps (abs sum))) + sum) + (when nil + (format t "n = ~D~%" n) + (format t " term = ~S~%" term) + (format t " sum = ~S~%" sum))))) + +;; Convert to iteration instead of this quick-and-dirty memoization? +(let ((hash (make-hash-table :test 'equal))) + (defun %big-a-clrhash () + (clrhash hash)) + (defun %big-a-dump-hash () + (maphash #'(lambda (k v) + (format t "~S -> ~S~%" k v)) + hash)) + (defun %big-a (n v) + (or (gethash (list n v) hash) + (let ((result + (cond ((zerop n) + (expt 2 (- v))) + (t + (* (%big-a (- n 1) v) + (/ (* (+ v n n -2) (+ v n n -1)) + (* 4 n (+ n v)))))))) + (setf (gethash (list n v) hash) result) + result)))) + +;; Computes A[n](v) = +;; (-1)^n*v*2^(-v)*pochhammer(v+n+1,n-1)/(2^(2*n)*n!) If v is a +;; negative integer -m, use A[n](-m) = (-1)^(m+1)*A[n-m](m) for n >= +;; m. +(defun big-a (n v) + (let ((m (ftruncate v))) + (cond ((and (= m v) (minusp m)) + (if (< n m) + (%big-a n v) + (let ((result (%big-a (+ n m) v))) + (if (oddp (truncate m)) + result + (- result))))) + (t + (%big-a n v))))) + +;; I[n](t, z, v) = exp(-t*z)/t^(2*n+v-1) * +;; integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf) +;; +;; Use the substitution u=1+s to get a new integral +;; +;; integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf) +;; = exp(t*z) * integrate(u^(-v-2*n)*exp(-t*u*z), u, 1, inf) +;; = exp(t*z)*t^(v+2*n-1)*z^(v+2*n-1)*incomplete_gamma_tail(1-v-2*n,t*z) +;; +;; The continued fraction for incomplete_gamma_tail(a,z) is +;; +;; z^a*exp(-z)/CF +;; +;; So incomplete_gamma_tail(1-v-2*n, t*z) is +;; +;; (t*z)^(1-v-2*n)*exp(-t*z)/CF +;; +;; which finally gives +;; +;; integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf) +;; = CF +;; +;; and I[n](t, z, v) = exp(-t*z)/t^(2*n+v-1)/CF +(defun big-i (n t z v) + (/ (exp (- (* t z))) + (expt t (+ n n v -1)) + (let* ((a (- 1 v n n)) + (z-a (- z a))) + (lentz #'(lambda (n) + (+ n n 1 z-a)) + #'(lambda (n) + (* n (- a n))))))) + +(defun sum-big-ia (big-n v z) + ) + +(defun bessel-j (v z) + (let ((vv (ftruncate v))) + (cond ((= vv v) + ;; v is an integer + (integer-bessel-j-exp-arc v z)) + (t + (let ((big-n 100) + (vpi (* v (float-pi (realpart z))))) + (+ (integer-bessel-j-exp-arc v z) + (* z + (/ (sin vpi) vpi) + (+ (/ -1 z) + (sum-ab big-n v z)))))))))  (defun paris-series (v z n) (labels ((pochhammer (a k) commit 1d404aca1e6652ad2456ba14e8bbdf5d329c06a0 Author: Raymond Toy <toy.raymond@gmail.com> Date: Fri Apr 6 19:28:23 2012 -0700 Add some comments, rename bessel-j-exp-arc to integer-bessel-j-exp-arc. diff --git a/qd-bessel.lisp b/qd-bessel.lisp index 95a61c0..e85d2d1 100644 --- a/qd-bessel.lisp +++ b/qd-bessel.lisp @@ -51,6 +51,22 @@ ;; B[k](p) = 1/2^(k+3/2)/p^(k+1/2)*p^(k+1/2)*exp(-p)/CF ;; = exp(-p)/2^(k+3/2)/CF ;; +;; +;; Note also that [2] gives a recurrence relationship for B[k](p) in +;; eq (2.6), but there is an error there. The correct relationship is +;; +;; B[k](p) = -exp(-p)/(p*sqrt(2)*2^(k+1)) + (k-1/2)*B[k-1](p)/(2*p) +;; +;; The paper is missing the division by p in the term containing +;; B[k-1](p). This is easily derived from the recurrence relationship +;; for the (lower) incomplete gamma function. +;; +;; Note too that as k increases, the recurrence appears to be unstable +;; and B[k](p) begins to increase even though it is strictly bounded. +;; (This is also easy to see from the integral.) Hence, we do not use +;; the recursion. However, it might be stable for use with +;; double-float precision; this has not been tested. +;; (defun bk (k p) (/ (exp (- p)) (* (sqrt (float 2 (realpart p))) (ash 1 (+ k 1))) @@ -157,7 +173,7 @@ ;; This currently only works for v an integer. ;; -(defun bessel-j-exp-arc (v z) +(defun integer-bessel-j-exp-arc (v z) (let* ((iz (* #c(0 1) z)) (i+ (exp-arc-i-2 iz v)) (i- (exp-arc-i-2 (- iz ) v))) commit 4c1ed0f45e79bbe467f7d4694f417ff92ae46adb Author: Raymond Toy <toy.raymond@gmail.com> Date: Fri Apr 6 19:21:46 2012 -0700 First cut at Bessel functions. Needs lots of work. diff --git a/qd-bessel.lisp b/qd-bessel.lisp new file mode 100644 index 0000000..95a61c0 --- /dev/null +++ b/qd-bessel.lisp @@ -0,0 +1,186 @@ +;;;; -*- Mode: lisp -*- +;;;; +;;;; Copyright (c) 2011 Raymond Toy +;;;; Permission is hereby granted, free of charge, to any person +;;;; obtaining a copy of this software and associated documentation +;;;; files (the "Software"), to deal in the Software without +;;;; restriction, including without limitation the rights to use, +;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell +;;;; copies of the Software, and to permit persons to whom the +;;;; Software is furnished to do so, subject to the following +;;;; conditions: +;;;; +;;;; The above copyright notice and this permission notice shall be +;;;; included in all copies or substantial portions of the Software. +;;;; +;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +;;;; OTHER DEALINGS IN THE SOFTWARE. + +(in-package #:oct) + +;;; References: +;;; +;;; [1] Borwein, Borwein, Crandall, "Effective Laguerre Asymptotics", +;;; http://people.reed.edu/~crandall/papers/Laguerre-f.pdf +;;; +;;; [2] Borwein, Borwein, Chan, "The Evaluation of Bessel Functions +;;; via Exp-Arc Integrals", http://web.cs.dal.ca/~jborwein/bessel.pdf +;;; + +(defvar *debug-exparc* nil) + +;; B[k](p) = 1/2^(k+3/2)*integrate(exp(-p*u)*u^(k-1/2),u,0,1) +;; = 1/2^(k+3/2)/p^(k+1/2)*integrate(t^(k-1/2)*exp(-t),t,0,p) +;; = 1/2^(k+3/2)/p^(k+1/2) * g(k+1/2, p) +;; +;; where g(a,z) is the lower incomplete gamma function. +;; +;; There is the continued fraction expansion for g(a,z) (see +;; cf-incomplete-gamma in qd-gamma.lisp): +;; +;; g(a,z) = z^a*exp(-z)/ CF +;; +;; So +;; +;; B[k](p) = 1/2^(k+3/2)/p^(k+1/2)*p^(k+1/2)*exp(-p)/CF +;; = exp(-p)/2^(k+3/2)/CF +;; +(defun bk (k p) + (/ (exp (- p)) + (* (sqrt (float 2 (realpart p))) (ash 1 (+ k 1))) + (let ((a (float (+ k 1/2) (realpart p)))) + (lentz #'(lambda (n) + (+ n a)) + #'(lambda (n) + (if (evenp n) + (* (ash n -1) p) + (- (* (+ a (ash n -1)) p)))))))) + +;; exp-arc I function, as given in the Laguerre paper +;; +;; I(p, q) = 4*exp(p) * sum(g[k](-2*%i*q)/(2*k)!*B[k](p), k, 0, inf) +;; +;; where g[k](p) = product(p^2+(2*j-1)^2, j, 1, k) and B[k](p) as above. +;; +;; For computation, note that g[k](p) = g[k-1](p) * (p^2 + (2*k-1)^2) +;; and (2*k)! = (2*k-2)! * (2*k-1) * (2*k). Then, let +;; +;; R[k](p) = g[k](p)/(2*k)! +;; +;; Then +;; +;; R[k](p) = g[k](p)/(2*k)! +;; = g[k-1](p)/(2*k-2)! * (p^2 + (2*k-1)^2)/((2*k-1)*(2*k) +;; = R[k-1](p) * (p^2 + (2*k-1)^2)/((2*k-1)*(2*k) +;; +;; In the exp-arc paper, the function is defined (equivalently) as +;; +;; I(p, q) = 2*%i*exp(p)/q * sum(r[2*k+1](-2*%i*q)/(2*k)!*B[k](p), k, 0, inf) +;; +;; where r[2*k+1](p) = p*product(p^2 + (2*j-1)^2, j, 1, k) +;; +;; Let's note some properties of I(p, q). +;; +;; I(-%i*z, v) = 2*%i*exp(-%i*z)/q * sum(r[2*k+1](-2*%i*v)/(2*k)!*B[k](-%i*z)) +;; +;; Note thate B[k](-%i*z) = 1/2^(k+3/2)*integrate(exp(%i*z*u)*u^(k-1/2),u,0,1) +;; = conj(B[k](%i*z). +;; +;; Hence I(-%i*z, v) = conj(I(%i*z, v)) when both z and v are real. +(defun exp-arc-i (p q) + (let* ((sqrt2 (sqrt (float 2 (realpart p)))) + (exp/p/sqrt2 (/ (exp (- p)) p sqrt2)) + (v (* #c(0 -2) q)) + (v2 (expt v 2)) + (eps (epsilon (realpart p)))) + (when *debug-exparc* + (format t "sqrt2 = ~S~%" sqrt2) + (format t "exp/p/sqrt2 = ~S~%" exp/p/sqrt2)) + (do* ((k 0 (1+ k)) + (bk (/ (incomplete-gamma 1/2 p) + 2 sqrt2 (sqrt p)) + (- (/ (* bk (- k 1/2)) 2 p) + (/ exp/p/sqrt2 (ash 1 (+ k 1))))) + ;; ratio[k] = r[2*k+1](v)/(2*k)!. + ;; r[1] = v and r[2*k+1](v) = r[2*k-1](v)*(v^2 + (2*k-1)^2) + ;; ratio[0] = v + ;; and ratio[k] = r[2*k-1](v)*(v^2+(2*k-1)^2) / ((2*k-2)! * (2*k-1) * 2*k) + ;; = ratio[k]*(v^2+(2*k-1)^2)/((2*k-1) * 2 * k) + (ratio v + (* ratio (/ (+ v2 (expt (1- (* 2 k)) 2)) + (* 2 k (1- (* 2 k)))))) + (term (* ratio bk) + (* ratio bk)) + (sum term (+ sum term))) + ((< (abs term) (* (abs sum) eps)) + (* sum #c(0 2) (/ (exp p) q))) + (when *debug-exparc* + (format t "k = ~D~%" k) + (format t " bk = ~S~%" bk) + (format t " ratio = ~S~%" ratio) + (format t " term = ~S~%" term) + (format t " sum - ~S~%" sum))))) + +(defun exp-arc-i-2 (p q) + (let* ((sqrt2 (sqrt (float 2 (realpart p)))) + (exp/p/sqrt2 (/ (exp (- p)) p sqrt2)) + (v (* #c(0 -2) q)) + (v2 (expt v 2)) + (eps (epsilon (realpart p)))) + (when *debug-exparc* + (format t "sqrt2 = ~S~%" sqrt2) + (format t "exp/p/sqrt2 = ~S~%" exp/p/sqrt2)) + (do* ((k 0 (1+ k)) + (bk (bk 0 p) + (bk k p)) + (ratio v + (* ratio (/ (+ v2 (expt (1- (* 2 k)) 2)) + (* 2 k (1- (* 2 k)))))) + (term (* ratio bk) + (* ratio bk)) + (sum term (+ sum term))) + ((< (abs term) (* (abs sum) eps)) + (* sum #c(0 2) (/ (exp p) q))) + (when *debug-exparc* + (format t "k = ~D~%" k) + (format t " bk = ~S~%" bk) + (format t " ratio = ~S~%" ratio) + (format t " term = ~S~%" term) + (format t " sum - ~S~%" sum))))) + + +;; This currently only works for v an integer. +;; +(defun bessel-j-exp-arc (v z) + (let* ((iz (* #c(0 1) z)) + (i+ (exp-arc-i-2 iz v)) + (i- (exp-arc-i-2 (- iz ) v))) + (/ (+ (* (cis (* v (float-pi i+) -1/2)) + i+) + (* (cis (* v (float-pi i+) 1/2)) + i-)) + (float-pi i+) + 2))) + +(defun paris-series (v z n) + (labels ((pochhammer (a k) + (/ (gamma (+ a k)) + (gamma a))) + (a (v k) + (* (/ (pochhammer (+ 1/2 v) k) + (gamma (float (1+ k) z))) + (pochhammer (- 1/2 v) k)))) + (* (loop for k from 0 below n + sum (* (/ (a v k) + (expt (* 2 z) k)) + (/ (cf-incomplete-gamma (+ k v 1/2) (* 2 z)) + (gamma (+ k v 1/2))))) + (/ (exp z) + (sqrt (* 2 (float-pi z) z)))))) + ----------------------------------------------------------------------- Summary of changes: oct.asd | 3 +- qd-bessel.lisp | 358 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ qd-gamma.lisp | 50 +++++--- qd-methods.lisp | 26 +++- qd-package.lisp | 2 + rt-tests.lisp | 16 +++ 6 files changed, 428 insertions(+), 27 deletions(-) create mode 100644 qd-bessel.lisp hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats