[oct-scm] [oct-git]OCT: A portable Lisp implementation for quad-double precision floats branch master updated. 139e1571a14f8cbdaa8da16aab06bf37ad98d3a3
Raymond Toy
rtoy at common-lisp.net
Thu Apr 12 16:30:41 UTC 2012
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 at 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 at 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
More information about the oct-cvs
mailing list