[Scipy-svn] r5516 - in trunk/scipy/special: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Fri Jan 23 19:32:03 EST 2009
Author: ptvirtan
Date: 2009-01-23 18:31:39 -0600 (Fri, 23 Jan 2009)
New Revision: 5516
Modified:
trunk/scipy/special/amos_wrappers.c
trunk/scipy/special/tests/test_basic.py
Log:
Fix #853: use correct symmetries for J_v / Y_v for negative non-int. orders
Modified: trunk/scipy/special/amos_wrappers.c
===================================================================
--- trunk/scipy/special/amos_wrappers.c 2009-01-23 23:14:26 UTC (rev 5515)
+++ trunk/scipy/special/amos_wrappers.c 2009-01-24 00:31:39 UTC (rev 5516)
@@ -73,6 +73,35 @@
return w;
}
+static Py_complex
+rotate_jy(Py_complex j, Py_complex y, double v)
+{
+ Py_complex w;
+ double c = cos(v * M_PI);
+ double s = sin(v * M_PI);
+ w.real = j.real * c - y.real * s;
+ w.imag = j.imag * c - y.imag * s;
+ return w;
+}
+
+static int
+reflect_jy(Py_complex *jy, double v)
+{
+ /* NB: Y_v may be huge near negative integers -- so handle exact
+ * integers carefully
+ */
+ int i;
+ if (v != floor(v))
+ return 0;
+
+ i = v - 16384.0 * floor(v / 16384.0);
+ if (i & 1) {
+ jy->real = -jy->real;
+ jy->imag = -jy->imag;
+ }
+ return 1;
+}
+
int cairy_wrap(Py_complex z, Py_complex *ai, Py_complex *aip, Py_complex *bi, Py_complex *bip) {
int id = 0;
int ierr = 0;
@@ -144,18 +173,22 @@
int kode = 1;
int nz, ierr;
int sign = 1;
- Py_complex cy;
+ Py_complex cy_j, cy_y, cwork;
if (v < 0) {
v = -v;
sign = -1;
}
- F_FUNC(zbesj,ZBESJ)(CADDR(z), &v, &kode, &n, CADDR(cy), &nz, &ierr);
+ F_FUNC(zbesj,ZBESJ)(CADDR(z), &v, &kode, &n, CADDR(cy_j), &nz, &ierr);
DO_MTHERR("jv:");
if (sign == -1) {
- cy = rotate(cy, v);
+ if (!reflect_jy(&cy_j, v)) {
+ F_FUNC(zbesy,ZBESY)(CADDR(z), &v, &kode, &n, CADDR(cy_y), &nz, CADDR(cwork), &ierr);
+ DO_MTHERR("jv(yv):");
+ cy_j = rotate_jy(cy_j, cy_y, v);
+ }
}
- return cy;
+ return cy_j;
}
Py_complex cbesj_wrap_e( double v, Py_complex z) {
@@ -163,18 +196,22 @@
int kode = 2;
int nz, ierr;
int sign = 1;
- Py_complex cy;
+ Py_complex cy_j, cy_y, cwork;
if (v < 0) {
v = -v;
sign = -1;
}
- F_FUNC(zbesj,ZBESJ)(CADDR(z), &v, &kode, &n, CADDR(cy), &nz, &ierr);
+ F_FUNC(zbesj,ZBESJ)(CADDR(z), &v, &kode, &n, CADDR(cy_j), &nz, &ierr);
DO_MTHERR("jve:");
if (sign == -1) {
- cy = rotate(cy, v);
+ if (!reflect_jy(&cy_j, v)) {
+ F_FUNC(zbesy,ZBESY)(CADDR(z), &v, &kode, &n, CADDR(cy_y), &nz, CADDR(cwork), &ierr);
+ DO_MTHERR("jve(yve):");
+ cy_j = rotate_jy(cy_j, cy_y, v);
+ }
}
- return cy;
+ return cy_j;
}
@@ -183,19 +220,23 @@
int kode = 1;
int nz, ierr;
int sign = 1;
- Py_complex cy, cwork;
+ Py_complex cy_y, cy_j, cwork;
if (v < 0) {
v = -v;
sign = -1;
}
- F_FUNC(zbesy,ZBESY)(CADDR(z), &v, &kode, &n, CADDR(cy), &nz, CADDR(cwork), &ierr);
+ F_FUNC(zbesy,ZBESY)(CADDR(z), &v, &kode, &n, CADDR(cy_y), &nz, CADDR(cwork), &ierr);
DO_MTHERR("yv:");
if (sign == -1) {
- cy = rotate(cy, v);
+ if (!reflect_jy(&cy_y, v)) {
+ F_FUNC(zbesj,ZBESJ)(CADDR(z), &v, &kode, &n, CADDR(cy_j), &nz, &ierr);
+ DO_MTHERR("yv(jv):");
+ cy_y = rotate_jy(cy_y, cy_j, -v);
+ }
}
- return cy;
+ return cy_y;
}
Py_complex cbesy_wrap_e( double v, Py_complex z) {
@@ -203,18 +244,22 @@
int kode = 2;
int nz, ierr;
int sign = 1;
- Py_complex cy, cwork;
+ Py_complex cy_y, cy_j, cwork;
if (v < 0) {
v = -v;
sign = -1;
}
- F_FUNC(zbesy,ZBESY)(CADDR(z), &v, &kode, &n, CADDR(cy), &nz, CADDR(cwork), &ierr);
+ F_FUNC(zbesy,ZBESY)(CADDR(z), &v, &kode, &n, CADDR(cy_y), &nz, CADDR(cwork), &ierr);
DO_MTHERR("yve:");
if (sign == -1) {
- cy = rotate(cy, v);
+ if (!reflect_jy(&cy_y, v)) {
+ F_FUNC(zbesj,ZBESJ)(CADDR(z), &v, &kode, &n, CADDR(cy_j), &nz, &ierr);
+ DO_MTHERR("yv(jv):");
+ cy_y = rotate_jy(cy_y, cy_j, -v);
+ }
}
- return cy;
+ return cy_y;
}
Modified: trunk/scipy/special/tests/test_basic.py
===================================================================
--- trunk/scipy/special/tests/test_basic.py 2009-01-23 23:14:26 UTC (rev 5515)
+++ trunk/scipy/special/tests/test_basic.py 2009-01-24 00:31:39 UTC (rev 5516)
@@ -1631,12 +1631,14 @@
class TestBesselJ(object):
- def test_cephes_vs_specfun(self):
- for v in [-120, -20., -10., -1., 0., 1., 12.49, 120., 301]:
+ def test_cephes_vs_amos(self):
+ for v in [-120, -100.3, -20., -10., -1., 0., 1., 12.49, 120., 301]:
for z in [1., 10., 200.5, 400., 600.5, 700.6, 1300, 10000]:
c1, c2, c3 = jv(v, z), jv(v,z+0j), jn(int(v), z)
if np.isinf(c1):
assert np.abs(c2) >= 1e150
+ elif np.isnan(c1):
+ assert np.abs(c2.imag) > 1e-10
else:
assert_tol_equal(c1, c2, err_msg=(v, z), rtol=1e-11)
if v == int(v):
@@ -1647,6 +1649,31 @@
assert_tol_equal(jv(301, 1300), 0.0183487151115275)
assert_tol_equal(jv(301, 1296.0682), -0.0224174325312048)
+ def test_ticket_853(self):
+ # cephes
+ assert_tol_equal(jv(-1, 1 ), -0.4400505857449335)
+ assert_tol_equal(jv(-2, 1 ), 0.1149034849319005)
+ assert_tol_equal(yv(-1, 1 ), 0.7812128213002887)
+ assert_tol_equal(yv(-2, 1 ), -1.650682606816255)
+ assert_tol_equal(jv(-0.5, 1 ), 0.43109886801837607952)
+ assert_tol_equal(yv(-0.5, 1 ), 0.6713967071418031)
+ # amos
+ assert_tol_equal(jv(-1, 1+0j), -0.4400505857449335)
+ assert_tol_equal(jv(-2, 1+0j), 0.1149034849319005)
+ assert_tol_equal(yv(-1, 1+0j), 0.7812128213002887)
+ assert_tol_equal(yv(-2, 1+0j), -1.650682606816255)
+ assert_tol_equal(jv(-0.5, 1+0j), 0.43109886801837607952)
+ assert_tol_equal(jv(-0.5, 1+1j), 0.2628946385649065-0.827050182040562j)
+ assert_tol_equal(yv(-0.5, 1+0j), 0.6713967071418031)
+ assert_tol_equal(yv(-0.5, 1+1j), 0.967901282890131+0.0602046062142816j)
+
+ assert_tol_equal(jve(-0.5,1+0j), jv(-0.5, 1+0j))
+ assert_tol_equal(yve(-0.5,1+0j), yv(-0.5, 1+0j))
+
+ assert_tol_equal(hankel1(-0.5, 1+1j), jv(-0.5, 1+1j) + 1j*yv(-0.5,1+1j))
+ assert_tol_equal(hankel2(-0.5, 1+1j), jv(-0.5, 1+1j) - 1j*yv(-0.5,1+1j))
+
+
class TestBesselI(object):
def _series(self, v, z, n=200):
More information about the Scipy-svn
mailing list