oct-cvs April 2012

oct-cvs@common-lisp.net

1 participants
16 discussions

[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

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

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

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

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

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

1 0

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

2009

2008

2007

oct-cvs April 2012