April 2012 - oct-cvs - mailman3.common-lisp.net

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. de65a54b545af76b6df16f3de36e54763ce97235
by Raymond Toy 17 Apr '12

17 Apr '12

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 de65a54b545af76b6df16f3de36e54763ce97235 (commit) via fb7212315ca5f592ef768c2636e7b3cd745473e4 (commit) via 581a0a08f04a985424c09b6a7b3661d2eb58c3e9 (commit) via 76ea9b133d74861db511d7c3e3a29653ce4ebda7 (commit) via 25129063aaa6168763b447db6f033359a19b60e0 (commit) via 22db89aa72af64360ed2d0f97709fc26ca01368b (commit) via fd281d5fcf9fcddb459b4fb12e6922ff2a2121a8 (commit) from e6a3f286122ea99d65a6448c3c41d0bdc411b2d7 (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 de65a54b545af76b6df16f3de36e54763ce97235 Merge: fb72123 e6a3f28 Author: Raymond Toy <rtoy(a)google.com> Date: Tue Apr 17 16:16:07 2012 -0700 Merge branch 'master' of ssh://common-lisp.net/var/git/projects/oct/oct commit fb7212315ca5f592ef768c2636e7b3cd745473e4 Author: Raymond Toy <rtoy(a)google.com> Date: Tue Apr 17 16:15:37 2012 -0700 Add function to compute J(1/2,z) and use it in the tests. diff --git a/rt-tests.lisp b/rt-tests.lisp index 5b9f5ad..af1d00c 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1617,6 +1617,11 @@ (check-accuracy 51.89 b true)) nil) +(defun bessel-j-1/2 (z) + ;; bessel_j(1/2,z) = sin(x)/sqrt(x)*sqrt(2/pi) + (* (/ (sin z) (sqrt z)) + (sqrt (/ 2 (float-pi z))))) + ;; Bessel J for half integer order and real args (rt:deftest bessel-j-1/2.d.1 (loop for k from 0 below 100 @@ -1625,7 +1630,7 @@ ;; = 0, where we will lose accuracy. for x = (+ 1 (random (/ pi 2))) for b = (bessel-j 0.5d0 x) - for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 pi))) + for true = (bessel-j-1/2 x) for result = (check-accuracy 48.42 b true) when result append (list (list (list k x) result))) @@ -1634,7 +1639,7 @@ (rt:deftest bessel-j-1/2.d.1.a (let* ((x 2.3831631289164497d0) (b (bessel-j 0.5d0 x)) - (true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 pi))))) + (true (bessel-j-1/2 x))) (check-accuracy 48.42 b true)) nil) @@ -1645,7 +1650,7 @@ ;; = 0, where we will lose accuracy. for x = (+ 1 (random (/ (float-pi #q1) 2))) for b = (bessel-j #q0.5 x) - for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1)))) + for true = (bessel-j-1/2 x) for result = (check-accuracy 169.45 b true) when result append (list (list (list k x) result))) @@ -1654,21 +1659,21 @@ (rt:deftest bessel-j-1/2.q.1.a (let* ((x #q1.1288834862545916200627583005758663687705443417892789067029865493882q0) (b (bessel-j #q0.5 x)) - (true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1)))))) + (true (bessel-j-1/2 x))) (check-accuracy 182.92 b true)) nil) (rt:deftest bessel-j-1/2.q.1.b (let* ((x #q1.1288834862545916200627583005758663687705443417892789067029865493882q0) (b (bessel-j #q0.5 x)) - (true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1)))))) + (true (bessel-j-1/2 x))) (check-accuracy 173.28 b true)) nil) (rt:deftest bessel-j-1/2.q.1.c (let* ((x #q1.0360263937639582798798376485114581552570020473846457752365459851056q0) (b (bessel-j #q0.5 x)) - (true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1)))))) + (true (bessel-j-1/2 x))) (check-accuracy 169.45 b true)) nil) @@ -1708,7 +1713,7 @@ for x = (complex (random (/ pi 2)) (random (/ pi 2))) for b = (bessel-j 0.5d0 x) - for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 pi))) + for true = (bessel-j-1/2 x) for result = (check-accuracy 49.8 b true) when result append (list (list (list k x) result))) @@ -1719,7 +1724,7 @@ for x = (complex (random (/ (float-pi #q1) 2)) (random (/ (float-pi #q1) 2))) for b = (bessel-j #q0.5 x) - for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1)))) + for true = (bessel-j-1/2 x) for result = (check-accuracy 212 b true) when result append (list (list (list k x) result))) commit 581a0a08f04a985424c09b6a7b3661d2eb58c3e9 Merge: 76ea9b1 78801d6 Author: Raymond Toy <rtoy(a)google.com> Date: Mon Apr 16 10:05:21 2012 -0700 Merge branch 'master' of ssh://common-lisp.net/var/git/projects/oct/oct commit 76ea9b133d74861db511d7c3e3a29653ce4ebda7 Merge: 2512906 5daefcc Author: Raymond Toy <rtoy(a)google.com> Date: Thu Apr 12 15:32:00 2012 -0700 Merge branch 'master' of ssh://common-lisp.net/var/git/projects/oct/oct commit 25129063aaa6168763b447db6f033359a19b60e0 Merge: 22db89a 7c5a318 Author: Raymond Toy <rtoy(a)google.com> Date: Wed Apr 11 13:53:07 2012 -0700 Merge branch 'master' of ssh://common-lisp.net/var/git/projects/oct/oct commit 22db89aa72af64360ed2d0f97709fc26ca01368b Merge: fd281d5 f8943af Author: Raymond Toy <rtoy(a)google.com> Date: Mon Apr 9 11:14:16 2012 -0700 Merge branch 'master' of ssh://common-lisp.net/var/git/projects/oct/oct diff --cc oct.asd index ccabbb1,0ca30e5..10cf0f9 --- a/oct.asd +++ b/oct.asd @@@ -66,9 -66,10 +66,11 @@@ (:file "qd-theta" :depends-on ("qd-methods" "qd-reader")) (:file "qd-gamma" - :depends-on ("qd-methods" "qd-reader")) + :depends-on ("qd-complex" "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 fd281d5fcf9fcddb459b4fb12e6922ff2a2121a8 Author: Raymond Toy <rtoy(a)google.com> Date: Fri Mar 23 20:48:40 2012 -0700 qd-gamma depends on qd-complex. diff --git a/oct.asd b/oct.asd index 3f8d70c..ccabbb1 100644 --- a/oct.asd +++ b/oct.asd @@ -66,7 +66,7 @@ (:file "qd-theta" :depends-on ("qd-methods" "qd-reader")) (:file "qd-gamma" - :depends-on ("qd-methods" "qd-reader")) + :depends-on ("qd-complex" "qd-methods" "qd-reader")) )) (defmethod perform ((op test-op) (c (eql (asdf:find-system :oct)))) ----------------------------------------------------------------------- Summary of changes: oct.asd | 3 ++- rt-tests.lisp | 21 +++++++++++++-------- 2 files changed, 15 insertions(+), 9 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. e6a3f286122ea99d65a6448c3c41d0bdc411b2d7
by Raymond Toy 17 Apr '12

17 Apr '12

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 e6a3f286122ea99d65a6448c3c41d0bdc411b2d7 (commit) from 64c424b7f41413424bbab93d4f5dd4f8461b784d (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 e6a3f286122ea99d65a6448c3c41d0bdc411b2d7 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Tue Apr 17 09:32:27 2012 -0700 sum-an was summing one too many terms. diff --git a/qd-bessel.lisp b/qd-bessel.lisp index 37a4c8e..bf3ca04 100644 --- a/qd-bessel.lisp +++ b/qd-bessel.lisp @@ -350,7 +350,7 @@ (an n k v))) (sum term (+ sum term))) ((or (<= (abs term) (* eps (abs sum))) - (> k big-n)) + (>= k big-n)) sum)))) ;; SUM-AB computes the series ----------------------------------------------------------------------- Summary of changes: qd-bessel.lisp | 2 +- 1 files changed, 1 insertions(+), 1 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 64c424b7f41413424bbab93d4f5dd4f8461b784d
by Raymond Toy 17 Apr '12

17 Apr '12

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 64c424b7f41413424bbab93d4f5dd4f8461b784d (commit) via 8c5195a87137fd61952a175a5dae676c70b480ef (commit) from 78801d6381aaaf4f21967582680f26889582db60 (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 64c424b7f41413424bbab93d4f5dd4f8461b784d Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Mon Apr 16 20:45:12 2012 -0700 Add iterative versions for some functions. diff --git a/qd-bessel.lisp b/qd-bessel.lisp index 14743af..37a4c8e 100644 --- a/qd-bessel.lisp +++ b/qd-bessel.lisp @@ -78,6 +78,21 @@ (* (ash n -1) p) (- (* (+ a (ash n -1)) p)))))))) +;; Use the recursion +(defun bk-iter (k p old-bk) + (with-floating-point-contagion (p old-bk) + (if (zerop k) + (* (sqrt (/ (float-pi p) 8)) + (let ((rp (sqrt p))) + (/ (erf rp) + rp))) + (- (* (- k 1/2) + (/ old-bk (* 2 p))) + (/ (exp (- p)) + p + (ash 1 (+ k 1)) + (sqrt (float 2 (realpart 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) @@ -178,6 +193,35 @@ (format t " term = ~S~%" term) (format t " sum - ~S~%" sum))))) +(defun exp-arc-i-3 (p q) + (let* ((v (* #c(0 -2) q)) + (v2 (expt v 2)) + (eps (epsilon (realpart p)))) + (do* ((k 0 (1+ k)) + (bk (bk 0 p) + (bk-iter k p bk)) + ;; Compute g[k](p)/(2*k)!, not r[2*k+1](p)/(2*k)! + (ratio 1 + (* 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)) + (when *debug-exparc* + (format t "Final k= ~D~%" k) + (format t " bk = ~S~%" bk) + (format t " ratio = ~S~%" ratio) + (format t " term = ~S~%" term) + (format t " sum - ~S~%" sum)) + (* sum 4 (exp p))) + (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))))) + ;; Not really just for Bessel J for integer orders, but in that case, ;; this is all that's needed to compute Bessel J. For other values, @@ -212,6 +256,7 @@ ;; ;; 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) +;; = (-1)^n/(2^n)*2*sinh(z/2)/z + n/z*beta[n-1](z) ;; ;; We also note that ;; @@ -223,12 +268,34 @@ (/ (incomplete-gamma (1+ n) (/ z 2)) (expt z (1+ n))))) +(defun alpha-iter (n z alpha-old) + (if (zerop n) + ;; (1- exp(-z/2))/z. + (/ (- 1 (exp (* z -1/2))) + z) + (- (* (/ n z) alpha-old) + (/ (exp (- (* z 1/2))) + z + (ash 1 n))))) + (defun beta (n z) (let ((n (float n (realpart z)))) (/ (- (incomplete-gamma (1+ n) (/ z 2)) (incomplete-gamma (1+ n) (/ z -2))) (expt z (1+ n))))) +(defun beta-iter (n z old-beta) + (if (zerop n) + ;; integrate(exp(-z*s),s,-1/2,1/2) + ;; = (exp(z/2)-exp(-z/2)/z + ;; = 2*sinh(z/2)/z + ;; = sinh(z/2)/(z/2) + (* 2 (/ (sinh (* 1/2 z)) z)) + (+ (* n (/ old-beta z)) + (* (/ (sinh (* 1/2 z)) (* 1/2 z)) + (scale-float (float (if (evenp n) 1 -1) (realpart z)) (- 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)); @@ -305,6 +372,26 @@ (format t " term = ~S~%" term) (format t " sum = ~S~%" sum))))) +(defun sum-ab-2 (big-n v z) + (let ((eps (epsilon (realpart z)))) + (an-clrhash) + (do* ((n 0 (+ 1 n)) + (alphan (alpha-iter 0 z 0) + (alpha-iter n z alphan)) + (betan (beta-iter 0 z 0) + (beta-iter n z betan)) + (term (+ (* alphan (an n 0 v)) + (* betan (sum-an big-n n v z))) + (+ (* alphan (an n 0 v)) + (* betan (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 () commit 8c5195a87137fd61952a175a5dae676c70b480ef Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Mon Apr 16 19:30:12 2012 -0700 Make big-n a defvar so we can change it easily. diff --git a/qd-bessel.lisp b/qd-bessel.lisp index 48a0dac..14743af 100644 --- a/qd-bessel.lisp +++ b/qd-bessel.lisp @@ -410,6 +410,7 @@ ;; bessel_i(v,z) = exp(-v*%pi*%i/2)*bessel_j(v, %i*z) ;; when Im(z) >> Re(z) ;; +(defvar *big-n* 100) (defun bessel-j (v z) (let ((vv (ftruncate v))) ;; Clear the caches for now. @@ -420,7 +421,7 @@ (integer-bessel-j-exp-arc v z)) (t ;; Need to fine-tune the value of big-n. - (let ((big-n 10) + (let ((big-n *big-n*) (vpi (* v (float-pi (realpart z))))) (+ (integer-bessel-j-exp-arc v z) (if (= vv v) ----------------------------------------------------------------------- Summary of changes: qd-bessel.lisp | 90 +++++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 files changed, 89 insertions(+), 1 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 78801d6381aaaf4f21967582680f26889582db60
by Raymond Toy 15 Apr '12

15 Apr '12

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 78801d6381aaaf4f21967582680f26889582db60 (commit) from 5566bc397937c2e8979ee6847e8f59f279f1f643 (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 78801d6381aaaf4f21967582680f26889582db60 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Sun Apr 15 10:50:16 2012 -0700 * Add comments for {{{integer-bessel-j-exp-arc}}}. * Simplify {{{sum-an}}} so we stop the sum when the terms no longer contribute to the sum. * Change {{{big-n}}}. This still needs work. diff --git a/qd-bessel.lisp b/qd-bessel.lisp index 0c06f8a..48a0dac 100644 --- a/qd-bessel.lisp +++ b/qd-bessel.lisp @@ -179,7 +179,13 @@ (format t " sum - ~S~%" sum))))) +;; Not really just for Bessel J for integer orders, but in that case, +;; this is all that's needed to compute Bessel J. For other values, +;; this is just part of the computation needed. ;; +;; Compute +;; +;; 1/(2*%pi) * (exp(-%i*v*%pi/2) * I(%i*z, v) + exp(%i*v*%pi/2) * I(-%i*z, v)) (defun integer-bessel-j-exp-arc (v z) (let* ((iz (* #c(0 1) z)) (i+ (exp-arc-i-2 iz v))) @@ -257,6 +263,7 @@ ;; ;; sum(exp(-k*z)*a[n](k,v), k, 1, N) ;; +#+nil (defun sum-an (big-n n v z) (let ((sum 0)) (loop for k from 1 upto big-n @@ -265,6 +272,20 @@ (an n k v)))) sum)) +;; Like above, but we just stop when the terms no longer contribute to +;; the sum. +(defun sum-an (big-n n v z) + (let ((eps (epsilon (realpart z)))) + (do* ((k 1 (+ 1 k)) + (term (* (exp (- (* k z))) + (an n k v)) + (* (exp (- (* k z))) + (an n k v))) + (sum term (+ sum term))) + ((or (<= (abs term) (* eps (abs sum))) + (> k big-n)) + 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) @@ -399,7 +420,7 @@ (integer-bessel-j-exp-arc v z)) (t ;; Need to fine-tune the value of big-n. - (let ((big-n 100) + (let ((big-n 10) (vpi (* v (float-pi (realpart z))))) (+ (integer-bessel-j-exp-arc v z) (if (= vv v) ----------------------------------------------------------------------- Summary of changes: qd-bessel.lisp | 23 ++++++++++++++++++++++- 1 files changed, 22 insertions(+), 1 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 5566bc397937c2e8979ee6847e8f59f279f1f643
by Raymond Toy 13 Apr '12

13 Apr '12

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 5566bc397937c2e8979ee6847e8f59f279f1f643 (commit) via 2f48eb872f71a8f41187decbc57878cb67fa86de (commit) from 5daefcc43f8b88741634f7a982bdb8b976739a1c (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 5566bc397937c2e8979ee6847e8f59f279f1f643 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Thu Apr 12 23:24:03 2012 -0700 Add more tests for Bessel J. Not all of them pass. diff --git a/rt-tests.lisp b/rt-tests.lisp index a4dac02..5b9f5ad 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1673,28 +1673,54 @@ nil) ;; Bessel J for complex args -(rt:deftest bessel-j-complex.pos-order.d.1 +(rt:deftest bessel-j-complex-arg.d.1 (let ((b (bessel-j 0d0 #c(1d0 1))) (true #c(0.9376084768060293d0 -0.4965299476091221d0))) (check-accuracy 50.73 b true)) nil) -(rt:deftest bessel-j-complex.pos-order.d.2 +(rt:deftest bessel-j-complex-arg.d.2 (let ((b (bessel-j 1d0 #c(1d0 1))) (true #c(0.6141603349229036d0 0.3650280288270878d0))) (check-accuracy 52.51 b true)) nil) -(rt:deftest bessel-j-complex.pos-order.d.3 +(rt:deftest bessel-j-complex-arg.d.3 (let ((b (bessel-j 2d0 #c(1d0 1))) (true #c(0.0415798869439621d0 0.2473976415133063d0))) (check-accuracy 50.41 b true)) nil) -#+nil -(rt:deftest bessel-j-complex.pos-order.d.4 +(rt:deftest bessel-j-complex-arg.d.4 (let ((b (bessel-j 2.3d0 #c(1d0 1))) + (true #c(-0.0141615213034667d0 0.1677798241687935d0))) + (check-accuracy 48.56 b true)) + nil) + +(rt:deftest bessel-j-complex-arg.d.5 + (let ((b (bessel-j -2.3d0 #c(1d0 1))) (true #c(0.1920598664138632d0 -0.5158676904105332d0))) - (check-accuracy 53 b true)) + (check-accuracy 50.97 b true)) nil) +(rt:deftest bessel-j-1/2-complex.d.1 + (loop for k from 0 below 10 + for x = (complex (random (/ pi 2)) + (random (/ pi 2))) + for b = (bessel-j 0.5d0 x) + for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 pi))) + for result = (check-accuracy 49.8 b true) + when result + append (list (list (list k x) result))) + nil) + +(rt:deftest bessel-j-1/2-complex.q.1 + (loop for k from 0 below 10 + for x = (complex (random (/ (float-pi #q1) 2)) + (random (/ (float-pi #q1) 2))) + for b = (bessel-j #q0.5 x) + for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1)))) + for result = (check-accuracy 212 b true) + when result + append (list (list (list k x) result))) + nil) commit 2f48eb872f71a8f41187decbc57878cb67fa86de Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Thu Apr 12 18:27:08 2012 -0700 Add more tests. diff --git a/rt-tests.lisp b/rt-tests.lisp index 712b96c..a4dac02 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1679,3 +1679,22 @@ (check-accuracy 50.73 b true)) nil) +(rt:deftest bessel-j-complex.pos-order.d.2 + (let ((b (bessel-j 1d0 #c(1d0 1))) + (true #c(0.6141603349229036d0 0.3650280288270878d0))) + (check-accuracy 52.51 b true)) + nil) + +(rt:deftest bessel-j-complex.pos-order.d.3 + (let ((b (bessel-j 2d0 #c(1d0 1))) + (true #c(0.0415798869439621d0 0.2473976415133063d0))) + (check-accuracy 50.41 b true)) + nil) + +#+nil +(rt:deftest bessel-j-complex.pos-order.d.4 + (let ((b (bessel-j 2.3d0 #c(1d0 1))) + (true #c(0.1920598664138632d0 -0.5158676904105332d0))) + (check-accuracy 53 b true)) + nil) + ----------------------------------------------------------------------- Summary of changes: rt-tests.lisp | 47 ++++++++++++++++++++++++++++++++++++++++++++++- 1 files changed, 46 insertions(+), 1 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 5daefcc43f8b88741634f7a982bdb8b976739a1c
by Raymond Toy 12 Apr '12

12 Apr '12

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 5daefcc43f8b88741634f7a982bdb8b976739a1c (commit) via 235ac2e4756b7c3f03815b267d292a26e8e1df37 (commit) from 0823acf49a91cbdb9a3e7f1faae2efe128580038 (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 5daefcc43f8b88741634f7a982bdb8b976739a1c Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Thu Apr 12 14:37:38 2012 -0700 Enable test of complex args. diff --git a/rt-tests.lisp b/rt-tests.lisp index 29a6501..712b96c 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1673,10 +1673,9 @@ nil) ;; Bessel J for complex args -#+nil (rt:deftest bessel-j-complex.pos-order.d.1 (let ((b (bessel-j 0d0 #c(1d0 1))) (true #c(0.9376084768060293d0 -0.4965299476091221d0))) - (check-accuracy 53 b true)) + (check-accuracy 50.73 b true)) nil) commit 235ac2e4756b7c3f03815b267d292a26e8e1df37 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Thu Apr 12 14:37:27 2012 -0700 Fix for complex args. diff --git a/qd-bessel.lisp b/qd-bessel.lisp index fb02f75..0c06f8a 100644 --- a/qd-bessel.lisp +++ b/qd-bessel.lisp @@ -183,18 +183,18 @@ (defun integer-bessel-j-exp-arc (v z) (let* ((iz (* #c(0 1) z)) (i+ (exp-arc-i-2 iz v))) - (cond ((= v (ftruncate v)) + (cond ((and (= v (ftruncate v)) (realp z)) ;; We can simplify the result - (let ((c (cis (* v (float-pi i+) -1/2)))) + (let ((c (exp (* v (float-pi i+) #c(0 -1/2))))) (/ (+ (* c i+) (* (conjugate c) (conjugate i+))) (float-pi i+) 2))) (t (let ((i- (exp-arc-i-2 (- iz ) v))) - (/ (+ (* (cis (* v (float-pi i+) -1/2)) + (/ (+ (* (exp (* v (float-pi i+) #c(0 -1/2))) i+) - (* (cis (* v (float-pi i+) 1/2)) + (* (exp (* v (float-pi i+) #c(0 1/2))) i-)) (float-pi i+) 2)))))) @@ -394,19 +394,21 @@ ;; Clear the caches for now. (an-clrhash) (%big-a-clrhash) - (cond ((= vv v) - ;; v is an integer + (cond ((and (= vv v) (realp z)) + ;; v is an integer and z is real (integer-bessel-j-exp-arc v z)) (t ;; Need to fine-tune the value of big-n. (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) - (sum-big-ia big-n v z))))))))) + (if (= vv v) + 0 + (* z + (/ (sin vpi) vpi) + (+ (/ -1 z) + (sum-ab big-n v z) + (sum-big-ia big-n v z)))))))))) ;; Bessel Y ;; ----------------------------------------------------------------------- Summary of changes: qd-bessel.lisp | 24 +++++++++++++----------- rt-tests.lisp | 3 +-- 2 files changed, 14 insertions(+), 13 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 0823acf49a91cbdb9a3e7f1faae2efe128580038
by Raymond Toy 12 Apr '12

12 Apr '12

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 0823acf49a91cbdb9a3e7f1faae2efe128580038 (commit) from 139e1571a14f8cbdaa8da16aab06bf37ad98d3a3 (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 0823acf49a91cbdb9a3e7f1faae2efe128580038 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Thu Apr 12 09:46:59 2012 -0700 Update accuracy for bessel-j-1/2.q test. diff --git a/rt-tests.lisp b/rt-tests.lisp index 86ae66c..29a6501 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1646,7 +1646,7 @@ for x = (+ 1 (random (/ (float-pi #q1) 2))) for b = (bessel-j #q0.5 x) for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1)))) - for result = (check-accuracy 173.28 b true) + for result = (check-accuracy 169.45 b true) when result append (list (list (list k x) result))) nil) @@ -1665,6 +1665,13 @@ (check-accuracy 173.28 b true)) nil) +(rt:deftest bessel-j-1/2.q.1.c + (let* ((x #q1.0360263937639582798798376485114581552570020473846457752365459851056q0) + (b (bessel-j #q0.5 x)) + (true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1)))))) + (check-accuracy 169.45 b true)) + nil) + ;; Bessel J for complex args #+nil (rt:deftest bessel-j-complex.pos-order.d.1 ----------------------------------------------------------------------- Summary of changes: rt-tests.lisp | 9 ++++++++- 1 files changed, 8 insertions(+), 1 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 139e1571a14f8cbdaa8da16aab06bf37ad98d3a3
by Raymond Toy 12 Apr '12

12 Apr '12

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 139e1571a14f8cbdaa8da16aab06bf37ad98d3a3 (commit) via 9fd2ebcbeed3ecae899f732b15fd279f6fb0f14f (commit) from 6ab5226ca3e32b443d87934ec138ff0efc8aaecc (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 139e1571a14f8cbdaa8da16aab06bf37ad98d3a3 Author: Raymond Toy <rtoy(a)google.com> Date: Thu Apr 12 09:30:28 2012 -0700 Add tests for bessel-j. diff --git a/rt-tests.lisp b/rt-tests.lisp index 325d202..86ae66c 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1542,3 +1542,134 @@ (check-accuracy 219.1 e true)) nil) + +;; Bessel J tests for negative order +(rt:deftest bessel-j.neg-order.d.1 + (let ((b (bessel-j -1d0 2d0)) + (true -0.5767248077568734d0)) + (check-accuracy 50.2 b true)) + nil) + +(rt:deftest bessel-j.neg-order.d.2 + (let ((b (bessel-j -1d0 1.5d0)) + (true -0.5579365079100996d0)) + (check-accuracy 50.5 b true)) + nil) + +(rt:deftest bessel-j.neg-order.d.3 + (let ((b (bessel-j -1.5d0 2d0)) + (true -0.3956232813587035d0)) + (check-accuracy 50.59 b true)) + nil) + +(rt:deftest bessel-j.neg-order.d.4 + (let ((b (bessel-j -1.8d0 1.5d0)) + (true -0.251327217627129314d0)) + (check-accuracy 49.98 b true)) + nil) + +(rt:deftest bessel-j.neg-order.d.5 + (let ((b (bessel-j -2d0 1.5d0)) + (true 0.2320876721442147d0)) + (check-accuracy 51.89 b true)) + nil) + +(rt:deftest bessel-j.neg-order.d.6 + (let ((b (bessel-j -2.5d0 1.5d0)) + (true 1.315037204805194d0)) + (check-accuracy 52.37 b true)) + nil) + +(rt:deftest bessel-j.neg-order.d.7 + (let ((b (bessel-j -2.3d0 1.5d0)) + (true 1.012178926325313d0)) + (check-accuracy 50.01 b true)) + nil) + +;; Bessel-J tests for positive order +(rt:deftest bessel-j.pos-order.d.1 + (let ((b (bessel-j 1.5d0 1d0)) + (true 0.2402978391234270d0)) + (check-accuracy 51.83 b true)) + nil) + +(rt:deftest bessel-j.pos-order.d.2 + (let ((b (bessel-j 1.8d0 1d0)) + (true 0.1564953153109239d0)) + (check-accuracy 51.97 b true)) + nil) + +(rt:deftest bessel-j.pos-order.d.3 + (let ((b (bessel-j 2d0 1d0)) + (true 0.1149034849319005d0)) + (check-accuracy 51.87 b true)) + nil) + +(rt:deftest bessel-j.pos-order.d.4 + (let ((b (bessel-j 2.5d0 1d0)) + (true 0.04949681022847794d0)) + (check-accuracy 47.17 b true)) + nil) + +(rt:deftest bessel-j.pos-order.d.5 + (let ((b (bessel-j -2d0 1.5d0)) + (true 0.2320876721442147d0)) + (check-accuracy 51.89 b true)) + nil) + +;; Bessel J for half integer order and real args +(rt:deftest bessel-j-1/2.d.1 + (loop for k from 0 below 100 + ;; x in [1,1+pi/2] because we don't want to test the Bessel + ;; series and we don't want to test near pi because sin(pi) + ;; = 0, where we will lose accuracy. + for x = (+ 1 (random (/ pi 2))) + for b = (bessel-j 0.5d0 x) + for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 pi))) + for result = (check-accuracy 48.42 b true) + when result + append (list (list (list k x) result))) + nil) + +(rt:deftest bessel-j-1/2.d.1.a + (let* ((x 2.3831631289164497d0) + (b (bessel-j 0.5d0 x)) + (true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 pi))))) + (check-accuracy 48.42 b true)) + nil) + +(rt:deftest bessel-j-1/2.q.1 + (loop for k from 0 below 10 + ;; x in [1,1+pi/2] because we don't want to test the Bessel + ;; series and we don't want to test near pi because sin(pi) + ;; = 0, where we will lose accuracy. + for x = (+ 1 (random (/ (float-pi #q1) 2))) + for b = (bessel-j #q0.5 x) + for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1)))) + for result = (check-accuracy 173.28 b true) + when result + append (list (list (list k x) result))) + nil) + +(rt:deftest bessel-j-1/2.q.1.a + (let* ((x #q1.1288834862545916200627583005758663687705443417892789067029865493882q0) + (b (bessel-j #q0.5 x)) + (true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1)))))) + (check-accuracy 182.92 b true)) + nil) + +(rt:deftest bessel-j-1/2.q.1.b + (let* ((x #q1.1288834862545916200627583005758663687705443417892789067029865493882q0) + (b (bessel-j #q0.5 x)) + (true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1)))))) + (check-accuracy 173.28 b true)) + nil) + +;; Bessel J for complex args +#+nil +(rt:deftest bessel-j-complex.pos-order.d.1 + (let ((b (bessel-j 0d0 #c(1d0 1))) + (true #c(0.9376084768060293d0 -0.4965299476091221d0))) + (check-accuracy 53 b true)) + nil) + commit 9fd2ebcbeed3ecae899f732b15fd279f6fb0f14f Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Thu Apr 12 09:30:10 2012 -0700 * Use the g[k] formula instead of r[2*k+1] because we fail to converge when v = 0. * Clear hash tables in bessel-j. diff --git a/qd-bessel.lisp b/qd-bessel.lisp index 588b5af..fb02f75 100644 --- a/qd-bessel.lisp +++ b/qd-bessel.lisp @@ -156,7 +156,8 @@ (do* ((k 0 (1+ k)) (bk (bk 0 p) (bk k p)) - (ratio v + ;; Compute g[k](p)/(2*k)!, not r[2*k+1](p)/(2*k)! + (ratio 1 (* ratio (/ (+ v2 (expt (1- (* 2 k)) 2)) (* 2 k (1- (* 2 k)))))) (term (* ratio bk) @@ -169,7 +170,7 @@ (format t " ratio = ~S~%" ratio) (format t " term = ~S~%" term) (format t " sum - ~S~%" sum)) - (* sum #c(0 2) (/ (exp p) q))) + (* sum 4 (exp p))) (when *debug-exparc* (format t "k = ~D~%" k) (format t " bk = ~S~%" bk) @@ -390,6 +391,9 @@ ;; (defun bessel-j (v z) (let ((vv (ftruncate v))) + ;; Clear the caches for now. + (an-clrhash) + (%big-a-clrhash) (cond ((= vv v) ;; v is an integer (integer-bessel-j-exp-arc v z)) ----------------------------------------------------------------------- Summary of changes: qd-bessel.lisp | 8 +++- rt-tests.lisp | 131 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 137 insertions(+), 2 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 6ab5226ca3e32b443d87934ec138ff0efc8aaecc
by Raymond Toy 12 Apr '12

12 Apr '12

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 6ab5226ca3e32b443d87934ec138ff0efc8aaecc (commit) via d795ba718dc53f591c82994811f50250aceec1d7 (commit) via c0e12ddf6b61f571555d46c6f168e6bebabc80b1 (commit) from 7c5a3186070096ee93e16f2ddf51b2c84e7c5895 (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 6ab5226ca3e32b443d87934ec138ff0efc8aaecc Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Wed Apr 11 20:32:23 2012 -0700 First cut at Bessel Y. Not working yet. diff --git a/qd-bessel.lisp b/qd-bessel.lisp index e7ca034..588b5af 100644 --- a/qd-bessel.lisp +++ b/qd-bessel.lisp @@ -376,6 +376,7 @@ (format t " term = ~S~%" term) (format t " sum = ~S~%" sum)))))) +;; ;; TODO: ;; o For |z| <= 1 use the series. ;; o Currently accuracy is not good for large z and half-integer @@ -402,6 +403,41 @@ (+ (/ -1 z) (sum-ab big-n v z) (sum-big-ia big-n v z))))))))) + +;; Bessel Y +;; +;; bessel_y(v, z) = 1/(2*%pi*%i)*(exp(-%i*v*%pi/2)*I(%i*v,z) - exp(%i*v*%pi/2)*I(-%i*z, v)) +;; + z/v/%pi*((1-cos(v*%pi)/z) + S(N,z,v)*cos(v*%pi)-S(N,z,-v)) +;; +;; where +;; +;; S(N,z,v) = sum(alpha[n](z)*a[n](0,v) + beta[n](z)*sum(exp(-k*z)*a[n](k,v),k,1,N),n,0,inf) +;; + sum(A[n](v)*I[n](N+1/2,z,v),n,0,inf) +;; +(defun bessel-y (v z) + (flet ((ipart (v z) + (let* ((iz (* #c(0 1) z)) + (c+ (exp (* v (float-pi z) 1/2))) + (c- (exp (* v (float-pi z) -1/2))) + (i+ (exp-arc-i-2 iz v)) + (i- (exp-arc-i-2 (- iz) v))) + (/ (- (* c- i+) (* c+ i-)) + (* #c(0 2) (float-pi z))))) + (s (big-n z v) + (+ (sum-ab big-n v z) + (sum-big-ia big-n v z)))) + (let* ((big-n 100) + (vpi (* v (float-pi z))) + (c (cos vpi))) + (+ (ipart v z) + (* (/ z vpi) + (+ (/ (- 1 c) + z) + (* c + (s big-n z v)) + (- (s big-n z (- v))))))))) + + (defun paris-series (v z n) (labels ((pochhammer (a k) commit d795ba718dc53f591c82994811f50250aceec1d7 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Wed Apr 11 20:15:47 2012 -0700 Add TODO list for bessel-j. diff --git a/qd-bessel.lisp b/qd-bessel.lisp index 3958ff2..e7ca034 100644 --- a/qd-bessel.lisp +++ b/qd-bessel.lisp @@ -375,13 +375,25 @@ (format t " f = ~S~%" f) (format t " term = ~S~%" term) (format t " sum = ~S~%" sum)))))) - + +;; TODO: +;; o For |z| <= 1 use the series. +;; o Currently accuracy is not good for large z and half-integer +;; order. +;; o For real v and z, return a real number instead of complex. +;; o Handle the case of Re(z) < 0. (The formulas are for Re(z) > 0: +;; bessel_j(v,z*exp(m*%pi*%i)) = exp(m*v*%pi*%i)*bessel_j(v, z) +;; o The paper suggests using +;; bessel_i(v,z) = exp(-v*%pi*%i/2)*bessel_j(v, %i*z) +;; when Im(z) >> Re(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 + ;; Need to fine-tune the value of big-n. (let ((big-n 100) (vpi (* v (float-pi (realpart z))))) (+ (integer-bessel-j-exp-arc v z) commit c0e12ddf6b61f571555d46c6f168e6bebabc80b1 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Wed Apr 11 20:06:38 2012 -0700 Correct the computation of alpha and beta. diff --git a/qd-bessel.lisp b/qd-bessel.lisp index 0197e7f..3958ff2 100644 --- a/qd-bessel.lisp +++ b/qd-bessel.lisp @@ -213,13 +213,13 @@ (defun alpha (n z) (let ((n (float n (realpart z)))) - (/ (cf-incomplete-gamma (1+ n) (/ z 2)) + (/ (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))) + (/ (- (incomplete-gamma (1+ n) (/ z 2)) + (incomplete-gamma (1+ n) (/ z -2))) (expt z (1+ n))))) ;; a[0](k,v) := (k+sqrt(k^2+1))^(-v); ----------------------------------------------------------------------- Summary of changes: qd-bessel.lisp | 56 ++++++++++++++++++++++++++++++++++++++++++++++++++++---- 1 files changed, 52 insertions(+), 4 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 7c5a3186070096ee93e16f2ddf51b2c84e7c5895
by Raymond Toy 11 Apr '12

11 Apr '12

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 7c5a3186070096ee93e16f2ddf51b2c84e7c5895 (commit) from bba9f8940c9f904bf14adc405d795a38ac333c24 (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 7c5a3186070096ee93e16f2ddf51b2c84e7c5895 Author: Raymond Toy <rtoy(a)google.com> Date: Wed Apr 11 09:18:31 2012 -0700 Correct some comments, remove unused code. diff --git a/qd-bessel.lisp b/qd-bessel.lisp index 367b847..0197e7f 100644 --- a/qd-bessel.lisp +++ b/qd-bessel.lisp @@ -37,7 +37,7 @@ ;; 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) +;; = 1/2^(k+3/2)/p^(k+1/2) * G(k+1/2, p) ;; ;; where G(a,z) is the lower incomplete gamma function. ;; @@ -109,6 +109,12 @@ ;; = conj(B[k](%i*z). ;; ;; Hence I(-%i*z, v) = conj(I(%i*z, v)) when both z and v are real. +;; +;; Also note that when v is an integer of the form (2*m+1)/2, then +;; r[2*k+1](-2*%i*v) = r[2*k+1](-%i*(2*m+1)) +;; = -%i*(2*m+1)*product(-(2*m+1)^2+(2*j-1)^2, j, 1, k) +;; so the product is zero when k >= m and the series I(p, q) is +;; finite. (defun exp-arc-i (p q) (let* ((sqrt2 (sqrt (float 2 (realpart p)))) (exp/p/sqrt2 (/ (exp (- p)) p sqrt2)) @@ -144,14 +150,9 @@ (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)) + (let* ((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)) @@ -162,6 +163,12 @@ (* ratio bk)) (sum term (+ sum term))) ((< (abs term) (* (abs sum) eps)) + (when *debug-exparc* + (format t "Final k= ~D~%" k) + (format t " bk = ~S~%" bk) + (format t " ratio = ~S~%" ratio) + (format t " term = ~S~%" term) + (format t " sum - ~S~%" sum)) (* sum #c(0 2) (/ (exp p) q))) (when *debug-exparc* (format t "k = ~D~%" k) ----------------------------------------------------------------------- Summary of changes: qd-bessel.lisp | 21 ++++++++++++++------- 1 files changed, 14 insertions(+), 7 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. bba9f8940c9f904bf14adc405d795a38ac333c24
by Raymond Toy 11 Apr '12

11 Apr '12

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 bba9f8940c9f904bf14adc405d795a38ac333c24 (commit) from 104efdeae6ef52c5d370e0b3a048a5087cdb1ea2 (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 bba9f8940c9f904bf14adc405d795a38ac333c24 Author: Raymond Toy <rtoy(a)google.com> Date: Wed Apr 11 08:47:28 2012 -0700 Update accuracy for oct.elliptic-pi.n2.d and add new test that caused oct.elliptic-pi.n2.d to fail. diff --git a/rt-tests.lisp b/rt-tests.lisp index 192b118..325d202 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1114,13 +1114,28 @@ for epi = (elliptic-pi n phi 0) for true = (/ (atanh (* (tan phi) (sqrt (- n 1)))) (sqrt (- n 1))) - for result = (check-accuracy 47 epi true) + for result = (check-accuracy 45.85 epi true) ;; Not sure if this formula holds when atanh gives a complex ;; result. Wolfram doesn't say when (and (not (complexp true)) result) append (list (list (list k n phi) result))) nil) +;; Failed test case: +;; ((89 66.68551748022054d0 0.12266024127708153d0) +;; (45.868614757480834d0 47 0.47787458521306514d0 +;; 0.4778745852130726d0)) +;; New threshold is 45.85 bits. +(rt:deftest oct.elliptic-pi.n2.d-1 + (let* ((n 66.68551748022054d0) + (phi 0.12266024127708153d0) + (epi (elliptic-pi n phi 0)) + (true (/ (atanh (* (tan phi) (sqrt (- n 1)))) + (sqrt (- n 1))))) + (check-accuracy 45.8686d0 epi true)) + nil) + + (rt:deftest oct.elliptic-pi.n0.q ;; Tests for random values for phi in [0, pi/2] and n in [0, 1] (loop for k from 0 below 100 @@ -1241,7 +1256,7 @@ ;; zero because value-or-tiny was too tiny. (let ((g (cf-incomplete-gamma 3d0 5d0)) (true (- 2 (* 37 (exp -5d0))))) - (check-accuracy 53 g true)) + (check-accuracy 51.2 g true)) nil) (rt:deftest gamma.1.d ----------------------------------------------------------------------- Summary of changes: rt-tests.lisp | 19 +++++++++++++++++-- 1 files changed, 17 insertions(+), 2 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 104efdeae6ef52c5d370e0b3a048a5087cdb1ea2
by Raymond Toy 11 Apr '12

11 Apr '12

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 104efdeae6ef52c5d370e0b3a048a5087cdb1ea2 (commit) from e9cd1a46fcf2a5c0101b3473648cb242b55987e8 (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 104efdeae6ef52c5d370e0b3a048a5087cdb1ea2 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Tue Apr 10 20:36:11 2012 -0700 VALUE-OR-TINY returned a value that was too tiny. * qd-gamma.lisp:: * Return sqrt(least-positive-normalized-double) instead of least-positive-normalized-double. * rt-tests.lisp:: * Add test for this case. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index 25fa38e..40a45eb 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -216,9 +216,9 @@ (incf tiny-value-count) (etypecase v ((or double-float cl:complex) - least-positive-normalized-double-float) + (sqrt least-positive-normalized-double-float)) ((or qd-real qd-complex) - (make-qd least-positive-normalized-double-float)))) + (make-qd (sqrt least-positive-normalized-double-float))))) v))) (let* ((f (value-or-tiny (funcall bf 0))) (c f) @@ -241,7 +241,7 @@ (setf d (/ d)) (setf f (* f delta)) (when *debug-cf-eval* - (format t " dl= ~S~%" delta) + (format t " dl= ~S (|dl - 1| = ~S)~%" delta (abs (1- delta))) (format t " f = ~S~%" f)) (when (<= (abs (- delta 1)) eps) (return-from lentz (values f j tiny-value-count))))) diff --git a/rt-tests.lisp b/rt-tests.lisp index 0544ea8..192b118 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -1234,6 +1234,16 @@ append (list (list (list k m) result))) nil) +(rt:deftest lentz + ;; This isn't really a test of cf-incomplete-gamma. It's a test + ;; that Lentz's algorithm works in this case. For these args, + ;; cf-incomplete-gamma used to generate an overflow or division by + ;; zero because value-or-tiny was too tiny. + (let ((g (cf-incomplete-gamma 3d0 5d0)) + (true (- 2 (* 37 (exp -5d0))))) + (check-accuracy 53 g true)) + nil) + (rt:deftest gamma.1.d (let ((g (gamma 0.5d0)) (true (sqrt pi))) ----------------------------------------------------------------------- Summary of changes: qd-gamma.lisp | 6 +++--- rt-tests.lisp | 10 ++++++++++ 2 files changed, 13 insertions(+), 3 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. e9cd1a46fcf2a5c0101b3473648cb242b55987e8
by Raymond Toy 10 Apr '12

10 Apr '12

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 e9cd1a46fcf2a5c0101b3473648cb242b55987e8 (commit) from f203389e5e78d1f001f68447ac2a9dd86dcfbbf6 (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 e9cd1a46fcf2a5c0101b3473648cb242b55987e8 Author: Raymond Toy <rtoy(a)google.com> Date: Tue Apr 10 09:51:55 2012 -0700 Fix bug in s-bessel-j. Microoptimize integer-bessel-j-exp-arc for the case where v is an integer. diff --git a/qd-bessel.lisp b/qd-bessel.lisp index 9759044..367b847 100644 --- a/qd-bessel.lisp +++ b/qd-bessel.lisp @@ -171,18 +171,25 @@ (format t " sum - ~S~%" sum))))) -;; This currently only works for v an integer. ;; (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))) - (/ (+ (* (cis (* v (float-pi i+) -1/2)) - i+) - (* (cis (* v (float-pi i+) 1/2)) - i-)) - (float-pi i+) - 2))) + (i+ (exp-arc-i-2 iz v))) + (cond ((= v (ftruncate v)) + ;; We can simplify the result + (let ((c (cis (* v (float-pi i+) -1/2)))) + (/ (+ (* c i+) + (* (conjugate c) (conjugate i+))) + (float-pi i+) + 2))) + (t + (let ((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)))))) ;; 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) @@ -349,15 +356,18 @@ (eps (epsilon z))) (do* ((k 0 (+ 1 k)) (f (gamma (+ v 1)) - (* f (* k (+ v k)))) + (* k (+ v k))) (term (/ f) (/ (* (- term) z2/4) f)) (sum term (+ sum term))) ((<= (abs term) (* eps (abs sum))) (* sum (expt (* z 1/2) v))) - (format t "k = ~D~%" k) - (format t " term = ~S~%" term) - (format t " sum = ~S~%" sum))))) + #+nil + (progn + (format t "k = ~D~%" k) + (format t " f = ~S~%" f) + (format t " term = ~S~%" term) + (format t " sum = ~S~%" sum)))))) (defun bessel-j (v z) (let ((vv (ftruncate v))) ----------------------------------------------------------------------- Summary of changes: qd-bessel.lisp | 36 +++++++++++++++++++++++------------- 1 files changed, 23 insertions(+), 13 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. f203389e5e78d1f001f68447ac2a9dd86dcfbbf6
by Raymond Toy 10 Apr '12

10 Apr '12

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 f203389e5e78d1f001f68447ac2a9dd86dcfbbf6 (commit) from f8943af6bff60e23d679089db5207b4834aa83ff (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 f203389e5e78d1f001f68447ac2a9dd86dcfbbf6 Author: Raymond Toy <toy.raymond(a)gmail.com> Date: Tue Apr 10 00:45:13 2012 -0700 * Implement sum-big-ia. * Add series for Bessel J. (Not working yet.) diff --git a/qd-bessel.lisp b/qd-bessel.lisp index fbe9760..9759044 100644 --- a/qd-bessel.lisp +++ b/qd-bessel.lisp @@ -324,8 +324,41 @@ (incomplete-gamma-tail a (* theta z))))) (defun sum-big-ia (big-n v z) - ) + (let ((big-n-1/2 (+ big-n 1/2)) + (eps (epsilon z))) + (do* ((n 0 (1+ n)) + (term (* (big-a 0 v) + (big-i 0 big-n-1/2 z v)) + (* (big-a n v) + (big-i n big-n-1/2 z v))) + (sum term (+ sum term))) + ((<= (abs term) (* eps (abs sum))) + sum) + #+nil + (progn + (format t "n = ~D~%" n) + (format t " term = ~S~%" term) + (format t " sum = ~S~%" sum))))) +;; Series for bessel J: +;; +;; (z/2)^v*sum((-1)^k/Gamma(k+v+1)/k!*(z^2//4)^k, k, 0, inf) +(defun s-bessel-j (v z) + (with-floating-point-contagion (v z) + (let ((z2/4 (* z z 1/4)) + (eps (epsilon z))) + (do* ((k 0 (+ 1 k)) + (f (gamma (+ v 1)) + (* f (* k (+ v k)))) + (term (/ f) + (/ (* (- term) z2/4) f)) + (sum term (+ sum term))) + ((<= (abs term) (* eps (abs sum))) + (* sum (expt (* z 1/2) v))) + (format t "k = ~D~%" k) + (format t " term = ~S~%" term) + (format t " sum = ~S~%" sum))))) + (defun bessel-j (v z) (let ((vv (ftruncate v))) (cond ((= vv v) @@ -338,7 +371,8 @@ (* z (/ (sin vpi) vpi) (+ (/ -1 z) - (sum-ab big-n v z))))))))) + (sum-ab big-n v z) + (sum-big-ia big-n v z))))))))) (defun paris-series (v z n) (labels ((pochhammer (a k) ----------------------------------------------------------------------- Summary of changes: qd-bessel.lisp | 38 ++++++++++++++++++++++++++++++++++++-- 1 files changed, 36 insertions(+), 2 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. f8943af6bff60e23d679089db5207b4834aa83ff
by Raymond Toy 09 Apr '12

09 Apr '12

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 f8943af6bff60e23d679089db5207b4834aa83ff (commit) via cb1a5d41baf9d4db12a2563a230d0a1e55c8adea (commit) from 6cfb0ac4b6bcc1a25bc119e87fd2b57bfa1f4355 (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 f8943af6bff60e23d679089db5207b4834aa83ff Author: Raymond Toy <rtoy(a)google.com> Date: Mon Apr 9 08:55:27 2012 -0700 Fix two typos in the names of tests. diff --git a/rt-tests.lisp b/rt-tests.lisp index 915f808..0544ea8 100644 --- a/rt-tests.lisp +++ b/rt-tests.lisp @@ -908,7 +908,7 @@ (check-accuracy 50 val 0d0)) nil) -(rt:deftest oct.jacobi-sn.1q +(rt:deftest oct.jacobi-cn.1q (let* ((ek (elliptic-k #q.5)) (val (jacobi-cn ek #q.5))) (check-accuracy 210 val #q0)) @@ -1179,7 +1179,7 @@ for m = (random #q1) for t3 = (elliptic-theta-3 0 (elliptic-nome m)) for true = (sqrt (/ (* 2 (elliptic-k m)) (float-pi m))) - for result = (check-accuracy 206 t3 true) + for result = (check-accuracy 205.7 t3 true) when result append (list (list (list k m) result))) nil) @@ -1411,7 +1411,7 @@ (check-accuracy 49.8 p true)) nil) -(rt:deftest psi.2d +(rt:deftest psi.2q (let* ((z (float 4/3 #q1)) (p (psi z)) (true (- 3 commit cb1a5d41baf9d4db12a2563a230d0a1e55c8adea Author: Raymond Toy <rtoy(a)google.com> Date: Mon Apr 9 08:52:08 2012 -0700 Fix some mistakes in exp-integral-e when changing algorithm to use series. diff --git a/qd-gamma.lisp b/qd-gamma.lisp index fe79813..25fa38e 100644 --- a/qd-gamma.lisp +++ b/qd-gamma.lisp @@ -532,36 +532,21 @@ ;; for |arg(z)| < pi. ;; ;; - (let* ((prec (float-contagion v z)) - (v (apply-contagion v prec)) - (z (apply-contagion z prec))) + (with-floating-point-contagion (v z) (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 + ((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 (+ -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)))) + (t + ;; Use continued fraction for everything else. + (cf-exp-integral-e v z))))) ;; Series for Fresnel S ;; ----------------------------------------------------------------------- Summary of changes: qd-gamma.lisp | 33 +++++++++------------------------ rt-tests.lisp | 6 +++--- 2 files changed, 12 insertions(+), 27 deletions(-) hooks/post-receive -- OCT: A portable Lisp implementation for quad-double precision floats

1 0

[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 6cfb0ac4b6bcc1a25bc119e87fd2b57bfa1f4355
by Raymond Toy 09 Apr '12

09 Apr '12

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(a)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(a)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(a)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(a)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(a)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(a)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(a)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(a)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(a)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(a)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(a)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(a)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

1 0