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

[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