[Scipy-svn] r5516 - in trunk/scipy/special: . tests

scipy-svn@scip... scipy-svn@scip...
Fri Jan 23 18:32:03 CST 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