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

scipy-svn@scip... scipy-svn@scip...
Sun Jan 25 13:02:49 CST 2009


Author: ptvirtan
Date: 2009-01-25 13:02:17 -0600 (Sun, 25 Jan 2009)
New Revision: 5523

Modified:
   trunk/scipy/special/_cephesmodule.c
   trunk/scipy/special/amos_wrappers.c
   trunk/scipy/special/amos_wrappers.h
   trunk/scipy/special/cephes/iv.c
   trunk/scipy/special/cephes/jv.c
   trunk/scipy/special/tests/test_basic.py
Log:
Fixing #854: make iv, jv, etc. return NaN when out-of-domain or failure

- Fix Cephes iv, jv return values; on DOMAIN return NaN, not zero

- Fix amos_wrappers return value handling: return NaNs when Amos signals
  that it didn't do any computation

- Fix _cephesmodule signatures: all real-valued functions should return
  NaN when the result would be complex, instead of just returning the
  real part.

Modified: trunk/scipy/special/_cephesmodule.c
===================================================================
--- trunk/scipy/special/_cephesmodule.c	2009-01-25 18:50:54 UTC (rev 5522)
+++ trunk/scipy/special/_cephesmodule.c	2009-01-25 19:02:17 UTC (rev 5523)
@@ -30,14 +30,12 @@
 static PyUFuncGenericFunction cephes1_2_functions[] = { NULL, NULL, NULL, NULL,};
 static PyUFuncGenericFunction cephes1_2c_functions[] = { NULL, NULL,};
 static PyUFuncGenericFunction cephes1c_4_functions[] = { NULL, NULL, NULL, NULL };
-static PyUFuncGenericFunction cephes1cp_4_functions[] = { NULL, NULL, NULL, NULL};
 static PyUFuncGenericFunction cephes1cpb_4_functions[] = { NULL, NULL,};
 static PyUFuncGenericFunction cephes2_functions[] = { NULL, NULL, };
 static PyUFuncGenericFunction cephes2_2_functions[] = { NULL, NULL, };
 static PyUFuncGenericFunction cephes2_4_functions[] = { NULL, NULL, };
 static PyUFuncGenericFunction cephes2a_functions[] = { NULL, NULL, };
 static PyUFuncGenericFunction cephes2c_functions[] = { NULL, NULL, NULL, NULL };
-static PyUFuncGenericFunction cephes2cp_functions[] = { NULL, NULL, NULL, NULL, };
 static PyUFuncGenericFunction cephes2cpp_functions[] = { NULL, NULL, };
 static PyUFuncGenericFunction cephes3_functions[] = { NULL, NULL, NULL, NULL};
 static PyUFuncGenericFunction cephes3a_functions[] = { NULL, NULL, };
@@ -50,7 +48,7 @@
 static PyUFuncGenericFunction cephes1c_functions[] = { NULL, NULL, };
 
 static void * airy_data[] = { (void *)airy, (void *)airy, (void *)cairy_wrap, (void *)cairy_wrap,};
-static void * airye_data[] = { (void *)cairy_wrap_e, (void *)cairy_wrap_e, (void *)cairy_wrap_e, (void *)cairy_wrap_e, };
+static void * airye_data[] = { (void *)cairy_wrap_e_real, (void *)cairy_wrap_e_real, (void *)cairy_wrap_e, (void *)cairy_wrap_e, };
 static void * itairy_data[] = { (void *)itairy_wrap, (void *)itairy_wrap, };
 
 static void * kelvin_data[] = { (void *)kelvin_wrap, (void *)kelvin_wrap,};
@@ -148,22 +146,22 @@
                                      (void *)gammaincinv, };
 
 static void * iv_data[] = { (void *)iv, (void *)iv, (void *)cbesi_wrap, (void *)cbesi_wrap,};
-static void * ive_data[] = { (void *)cbesi_wrap_e, (void *)cbesi_wrap_e, (void *)cbesi_wrap_e, (void *)cbesi_wrap_e, };
+static void * ive_data[] = { (void *)cbesi_wrap_e_real, (void *)cbesi_wrap_e_real, (void *)cbesi_wrap_e, (void *)cbesi_wrap_e, };
 static void * j0_data[] = { (void *)j0,  (void *)j0,  };
 static void * y0_data[] = { (void *)y0, (void *)y0, };
 static void * j1_data[] = { (void *)j1,  (void *)j1,  };
 static void * y1_data[] = { (void *)y1, (void *)y1, };
 static void * jv_data[] = { (void *)jv, (void *)jv, (void *)cbesj_wrap, (void *)cbesj_wrap,};
-static void * jve_data[] = { (void *)cbesj_wrap_e, (void *)cbesj_wrap_e, (void *)cbesj_wrap_e, (void *)cbesj_wrap_e, };
+static void * jve_data[] = { (void *)cbesj_wrap_e_real, (void *)cbesj_wrap_e_real, (void *)cbesj_wrap_e, (void *)cbesj_wrap_e, };
 static void * yv_data[] = { (void *)yv, (void *)yv, (void *)cbesy_wrap, (void *)cbesy_wrap,};
-static void * yve_data[] = { (void *)cbesy_wrap_e, (void *)cbesy_wrap_e, (void *)cbesy_wrap_e, (void *)cbesy_wrap_e, };
+static void * yve_data[] = { (void *)cbesy_wrap_e_real, (void *)cbesy_wrap_e_real, (void *)cbesy_wrap_e, (void *)cbesy_wrap_e, };
 
 static void * k0_data[] = { (void *)k0, (void *)k0, };
 static void * k0e_data[] = { (void *)k0e, (void *)k0e, };
 static void * k1_data[] = { (void *)k1, (void *)k1, };
 static void * k1e_data[] = { (void *)k1e, (void *)k1e, };
-static void * kv_data[] = { (void *)cbesk_wrap, (void *)cbesk_wrap, (void *)cbesk_wrap, (void *)cbesk_wrap,};
-static void * kve_data[] = { (void *)cbesk_wrap_e, (void *)cbesk_wrap_e, (void *)cbesk_wrap_e, (void *)cbesk_wrap_e,};
+static void * kv_data[] = { (void *)cbesk_wrap_real, (void *)cbesk_wrap_real, (void *)cbesk_wrap, (void *)cbesk_wrap,};
+static void * kve_data[] = { (void *)cbesk_wrap_e_real, (void *)cbesk_wrap_e_real, (void *)cbesk_wrap_e, (void *)cbesk_wrap_e,};
 static void * hankel1_data[] = { (void *)cbesh_wrap1, (void *)cbesh_wrap1,};
 static void * hankel1e_data[] = { (void *)cbesh_wrap1_e, (void *)cbesh_wrap1_e,};
 static void * hankel2_data[] = { (void *)cbesh_wrap2, (void *)cbesh_wrap2,};
@@ -356,10 +354,6 @@
         cephes1c_4_functions[1] = PyUFunc_d_dddd;
         cephes1c_4_functions[2] = PyUFunc_F_FFFF_As_D_DDDD;
         cephes1c_4_functions[3] = PyUFunc_D_DDDD;
-        cephes1cp_4_functions[0] = PyUFunc_f_ffff_As_D_DDDD;
-        cephes1cp_4_functions[1] = PyUFunc_d_dddd_As_D_DDDD;
-        cephes1cp_4_functions[2] = PyUFunc_F_FFFF_As_D_DDDD;
-        cephes1cp_4_functions[3] = PyUFunc_D_DDDD;
         cephes1cpb_4_functions[0] = PyUFunc_f_FFFF_As_d_DDDD;
         cephes1cpb_4_functions[1] = PyUFunc_d_DDDD;
         cephes2_functions[0] = PyUFunc_ff_f_As_dd_d;
@@ -372,10 +366,6 @@
         cephes2c_functions[1] = PyUFunc_dd_d;
         cephes2c_functions[2] = PyUFunc_fF_F_As_dD_D;
         cephes2c_functions[3] = PyUFunc_dD_D;
-        cephes2cp_functions[0] = PyUFunc_ff_f_As_dD_D;
-        cephes2cp_functions[1] = PyUFunc_dd_d_As_dD_D;
-        cephes2cp_functions[2] = PyUFunc_fF_F_As_dD_D;
-        cephes2cp_functions[3] = PyUFunc_dD_D;
         cephes2cpp_functions[0] = PyUFunc_fF_F_As_dD_D;
         cephes2cpp_functions[1] = PyUFunc_dD_D;
         cephes2_4_functions[0] = PyUFunc_ff_ffff_As_dd_dddd;
@@ -555,7 +545,7 @@
         f = PyUFunc_FromFuncAndData(cephes2c_functions, iv_data, cephes_3c_types, 4, 2, 1, PyUFunc_None, "iv", iv_doc, 0);
 	PyDict_SetItemString(dictionary, "iv", f);
 	Py_DECREF(f);
-	f = PyUFunc_FromFuncAndData(cephes2cp_functions, ive_data, cephes_3c_types, 4, 2, 1, PyUFunc_None, "ive", ive_doc, 0);
+	f = PyUFunc_FromFuncAndData(cephes2c_functions, ive_data, cephes_3c_types, 4, 2, 1, PyUFunc_None, "ive", ive_doc, 0);
 	PyDict_SetItemString(dictionary, "ive", f);
 	Py_DECREF(f);
 
@@ -614,7 +604,7 @@
 	Py_DECREF(f);
 
 
-	f = PyUFunc_FromFuncAndData(cephes1cp_4_functions, airye_data, cephes_5c_types, 4, 1, 4, PyUFunc_None, "airye", airye_doc, 0);
+	f = PyUFunc_FromFuncAndData(cephes1c_4_functions, airye_data, cephes_5c_types, 4, 1, 4, PyUFunc_None, "airye", airye_doc, 0);
 	PyDict_SetItemString(dictionary, "airye", f);
 	Py_DECREF(f);
 
@@ -662,13 +652,13 @@
            accurate. So we alias jv to jn */
         PyDict_SetItemString(dictionary, "jn", f);
 	Py_DECREF(f);
-	f = PyUFunc_FromFuncAndData(cephes2cp_functions, jve_data, cephes_3c_types, 4, 2, 1, PyUFunc_None, "jve", jve_doc, 0);
+	f = PyUFunc_FromFuncAndData(cephes2c_functions, jve_data, cephes_3c_types, 4, 2, 1, PyUFunc_None, "jve", jve_doc, 0);
 	PyDict_SetItemString(dictionary, "jve", f);
 	Py_DECREF(f);
 	f = PyUFunc_FromFuncAndData(cephes2c_functions, yv_data, cephes_3c_types, 4, 2, 1, PyUFunc_None, "yv", yv_doc, 0);
 	PyDict_SetItemString(dictionary, "yv", f);
 	Py_DECREF(f);
-	f = PyUFunc_FromFuncAndData(cephes2cp_functions, yve_data, cephes_3c_types, 4, 2, 1, PyUFunc_None, "yve", yve_doc, 0);
+	f = PyUFunc_FromFuncAndData(cephes2c_functions, yve_data, cephes_3c_types, 4, 2, 1, PyUFunc_None, "yve", yve_doc, 0);
 	PyDict_SetItemString(dictionary, "yve", f);
 	Py_DECREF(f);
 
@@ -685,10 +675,10 @@
 	f = PyUFunc_FromFuncAndData(cephes1_functions, k1e_data, cephes_2_types, 2, 1, 1, PyUFunc_None, "k1e", k1e_doc, 0);
 	PyDict_SetItemString(dictionary, "k1e", f);
 	Py_DECREF(f);
-	f = PyUFunc_FromFuncAndData(cephes2cp_functions, kv_data, cephes_3c_types, 4, 2, 1, PyUFunc_None, "kv", kv_doc, 0);
+	f = PyUFunc_FromFuncAndData(cephes2c_functions, kv_data, cephes_3c_types, 4, 2, 1, PyUFunc_None, "kv", kv_doc, 0);
 	PyDict_SetItemString(dictionary, "kv", f);
 	Py_DECREF(f);
-	f = PyUFunc_FromFuncAndData(cephes2cp_functions, kve_data, cephes_3c_types, 4, 2, 1, PyUFunc_None, "kve", kve_doc, 0);
+	f = PyUFunc_FromFuncAndData(cephes2c_functions, kve_data, cephes_3c_types, 4, 2, 1, PyUFunc_None, "kve", kve_doc, 0);
 	PyDict_SetItemString(dictionary, "kve", f);
 	Py_DECREF(f);
 
@@ -1064,7 +1054,7 @@
     va_end(ap);
 
     NPY_ALLOW_C_API
-    PyErr_WarnEx(scipy_special_SpecialFunctionWarning, msg, 2);
+    PyErr_WarnEx(scipy_special_SpecialFunctionWarning, msg, 1);
     NPY_DISABLE_C_API
 }
 

Modified: trunk/scipy/special/amos_wrappers.c
===================================================================
--- trunk/scipy/special/amos_wrappers.c	2009-01-25 18:50:54 UTC (rev 5522)
+++ trunk/scipy/special/amos_wrappers.c	2009-01-25 19:02:17 UTC (rev 5523)
@@ -62,6 +62,13 @@
   return -1;
 }
 
+void set_nan_if_no_computation_done(Py_complex *v, int ierr) {
+  if (v != NULL && (ierr == 1 || ierr == 2 || ierr == 4 || ierr == 5)) {
+    v->real = NAN;
+    v->imag = NAN;
+  }
+}
+
 static Py_complex
 rotate(Py_complex z, double v)
 {
@@ -127,15 +134,15 @@
   int nz;
 
   F_FUNC(zairy,ZAIRY)(CADDR(z), &id, &kode, F2C_CST(ai), &nz, &ierr);
-  DO_MTHERR("airy:");
+  DO_MTHERR("airy:", ai);
   F_FUNC(zbiry,ZBIRY)(CADDR(z), &id, &kode, F2C_CST(bi), &nz, &ierr);
-  DO_MTHERR("airy:");
+  DO_MTHERR("airy:", bi);
   
   id = 1;
   F_FUNC(zairy,ZAIRY)(CADDR(z), &id, &kode, F2C_CST(aip), &nz, &ierr);
-  DO_MTHERR("airy:");
+  DO_MTHERR("airy:", aip);
   F_FUNC(zbiry,ZBIRY)(CADDR(z), &id, &kode, F2C_CST(bip), &nz, &ierr);
-  DO_MTHERR("airy:");
+  DO_MTHERR("airy:", bip);
   return 0;
 }
 
@@ -145,18 +152,52 @@
   int nz, ierr;
 
   F_FUNC(zairy,ZAIRY)(CADDR(z), &id, &kode, F2C_CST(ai), &nz, &ierr);
-  DO_MTHERR("airye:");
+  DO_MTHERR("airye:", ai);
   F_FUNC(zbiry,ZBIRY)(CADDR(z), &id, &kode, F2C_CST(bi), &nz, &ierr);
-  DO_MTHERR("airye:");
+  DO_MTHERR("airye:", bi);
   
   id = 1;
   F_FUNC(zairy,ZAIRY)(CADDR(z), &id, &kode, F2C_CST(aip), &nz, &ierr);
-  DO_MTHERR("airye:");
+  DO_MTHERR("airye:", aip);
   F_FUNC(zbiry,ZBIRY)(CADDR(z), &id, &kode, F2C_CST(bip), &nz, &ierr);
-  DO_MTHERR("airye:");
+  DO_MTHERR("airye:", bip);
   return 0;
 }
 
+int cairy_wrap_e_real(double z, double *ai, double *aip, double *bi, double *bip) {
+  int id = 0;
+  int kode = 2;        /* Exponential scaling */
+  int nz, ierr;
+  Py_complex cz, cai, caip, cbi, cbip;
+
+  cz.real = z;
+  cz.imag = 0;
+
+  if (z < 0) {
+      *ai = NAN;
+  } else {
+      F_FUNC(zairy,ZAIRY)(CADDR(cz), &id, &kode, CADDR(cai), &nz, &ierr);
+      DO_MTHERR("airye:", &cai);
+      *ai = cai.real;
+  }
+  F_FUNC(zbiry,ZBIRY)(CADDR(cz), &id, &kode, CADDR(cbi), &nz, &ierr);
+  DO_MTHERR("airye:", &cbi);
+  *bi = cbi.real;
+  
+  id = 1;
+  if (z < 0) {
+      *aip = NAN;
+  } else {
+      F_FUNC(zairy,ZAIRY)(CADDR(cz), &id, &kode, CADDR(caip), &nz, &ierr);
+      DO_MTHERR("airye:", &caip);
+      *aip = caip.real;
+  }
+  F_FUNC(zbiry,ZBIRY)(CADDR(cz), &id, &kode, CADDR(cbip), &nz, &ierr);
+  DO_MTHERR("airye:", &cbip);
+  *bip = cbip.real;
+  return 0;
+}
+
 Py_complex cbesi_wrap( double v, Py_complex z) {
   int n = 1;
   int kode = 1;
@@ -169,16 +210,16 @@
     sign = -1;
   }
   F_FUNC(zbesi,ZBESI)(CADDR(z), &v,  &kode, &n, CADDR(cy), &nz, &ierr);
-  DO_MTHERR("iv:");
+  DO_MTHERR("iv:", &cy);
 
   if (sign == -1) {
     if (!reflect_i(&cy, v)) {
       F_FUNC(zbesk,ZBESK)(CADDR(z), &v,  &kode, &n, CADDR(cy_k), &nz, &ierr);
-      DO_MTHERR("iv(kv):");
+      DO_MTHERR("iv(kv):", &cy_k);
       cy = rotate_i(cy, cy_k, v);
     }
   }
-      
+
   return cy;
 }
 
@@ -194,12 +235,12 @@
     sign = -1;
   }
   F_FUNC(zbesi,ZBESI)(CADDR(z), &v,  &kode, &n, CADDR(cy), &nz, &ierr);
-  DO_MTHERR("ive:");
+  DO_MTHERR("ive:", &cy);
 
   if (sign == -1) {
     if (!reflect_i(&cy, v)) {
       F_FUNC(zbesk,ZBESK)(CADDR(z), &v,  &kode, &n, CADDR(cy_k), &nz, &ierr);
-      DO_MTHERR("ive(kv):");
+      DO_MTHERR("ive(kv):", &cy_k);
       /* adjust scaling to match zbesi */
       cy_k = rotate(cy_k, -z.imag/M_PI);
       if (z.real > 0) {
@@ -214,6 +255,17 @@
   return cy;
 }
 
+double cbesi_wrap_e_real(double v, double z) {
+  Py_complex cy, w;
+  if (v != floor(v) && z < 0) {
+    return NAN;
+  } else {
+    w.real = z;
+    w.imag = 0;
+    cy = cbesi_wrap_e(v, w);
+    return cy.real;
+  }
+}
   
 Py_complex cbesj_wrap( double v, Py_complex z) {
   int n = 1;
@@ -227,12 +279,12 @@
     sign = -1;
   }
   F_FUNC(zbesj,ZBESJ)(CADDR(z), &v,  &kode, &n, CADDR(cy_j), &nz, &ierr);
-  DO_MTHERR("jv:");
+  DO_MTHERR("jv:", &cy_j);
 
   if (sign == -1) {
     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):");
+      DO_MTHERR("jv(yv):", &cy_y);
       cy_j = rotate_jy(cy_j, cy_y, v);
     }
   }
@@ -251,17 +303,28 @@
     sign = -1;
   }
   F_FUNC(zbesj,ZBESJ)(CADDR(z), &v, &kode, &n, CADDR(cy_j), &nz, &ierr);
-  DO_MTHERR("jve:");
+  DO_MTHERR("jve:", &cy_j);
   if (sign == -1) {
     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):");
+      DO_MTHERR("jve(yve):", &cy_y);
       cy_j = rotate_jy(cy_j, cy_y, v);
     }
   }
   return cy_j;
 }
 
+double cbesj_wrap_e_real(double v, double z) {
+  Py_complex cy, w;
+  if (v != floor(v) && z < 0) {
+    return NAN;
+  } else {
+    w.real = z;
+    w.imag = 0;
+    cy = cbesj_wrap_e(v, w);
+    return cy.real;
+  }
+}
   
 Py_complex cbesy_wrap( double v, Py_complex z) {
   int n = 1;
@@ -276,11 +339,11 @@
   }
   F_FUNC(zbesy,ZBESY)(CADDR(z), &v,  &kode, &n, CADDR(cy_y), &nz, CADDR(cwork), &ierr);
 
-  DO_MTHERR("yv:");
+  DO_MTHERR("yv:", &cy_y);
   if (sign == -1) {
     if (!reflect_jy(&cy_y, v)) {
       F_FUNC(zbesj,ZBESJ)(CADDR(z), &v,  &kode, &n, CADDR(cy_j), &nz, &ierr);
-      DO_MTHERR("yv(jv):");
+      DO_MTHERR("yv(jv):", &cy_j);
       cy_y = rotate_jy(cy_y, cy_j, -v);
     }
   }
@@ -299,17 +362,28 @@
     sign = -1;
   }
   F_FUNC(zbesy,ZBESY)(CADDR(z), &v, &kode, &n, CADDR(cy_y), &nz, CADDR(cwork), &ierr);
-  DO_MTHERR("yve:");
+  DO_MTHERR("yve:", &cy_y);
   if (sign == -1) {
     if (!reflect_jy(&cy_y, v)) {
       F_FUNC(zbesj,ZBESJ)(CADDR(z), &v,  &kode, &n, CADDR(cy_j), &nz, &ierr);
-      DO_MTHERR("yv(jv):");
+      DO_MTHERR("yv(jv):", &cy_j);
       cy_y = rotate_jy(cy_y, cy_j, -v);
     }
   }
   return cy_y;
 }
 
+double cbesy_wrap_e_real(double v, double z) {
+  Py_complex cy, w;
+  if (z < 0) {
+    return NAN;
+  } else {
+    w.real = z;
+    w.imag = 0;
+    cy = cbesy_wrap_e(v, w);
+    return cy.real;
+  }
+}
   
 Py_complex cbesk_wrap( double v, Py_complex z) {
   int n = 1;
@@ -322,7 +396,7 @@
     v = -v;
   }
   F_FUNC(zbesk,ZBESK)(CADDR(z), &v,  &kode, &n, CADDR(cy), &nz, &ierr);
-  DO_MTHERR("kv:");
+  DO_MTHERR("kv:", &cy);
   return cy;
 }
 
@@ -337,10 +411,34 @@
     v = -v;
   }
   F_FUNC(zbesk,ZBESK)(CADDR(z), &v, &kode, &n, CADDR(cy), &nz, &ierr);
-  DO_MTHERR("kve:");
+  DO_MTHERR("kve:", &cy);
   return cy;
 }
   
+double cbesk_wrap_real( double v, double z) {
+  Py_complex cy, w;
+  if (z < 0) {
+    return NAN;
+  } else {
+    w.real = z;
+    w.imag = 0;
+    cy = cbesk_wrap(v, w);
+    return cy.real;
+  }
+}
+
+double cbesk_wrap_e_real( double v, double z) {
+  Py_complex cy, w;
+  if (z < 0) {
+    return NAN;
+  } else {
+    w.real = z;
+    w.imag = 0;
+    cy = cbesk_wrap_e(v, w);
+    return cy.real;
+  }
+}
+  
 Py_complex cbesh_wrap1( double v, Py_complex z) {
   int n = 1;
   int kode = 1;
@@ -354,7 +452,7 @@
     sign = -1;
   }
   F_FUNC(zbesh,ZBESH)(CADDR(z), &v,  &kode, &m, &n, CADDR(cy), &nz, &ierr);
-  DO_MTHERR("hankel1:");
+  DO_MTHERR("hankel1:", &cy);
   if (sign == -1) {
     cy = rotate(cy, v);
   }
@@ -374,7 +472,7 @@
     sign = -1;
   }
   F_FUNC(zbesh,ZBESH)(CADDR(z), &v, &kode, &m, &n, CADDR(cy), &nz, &ierr);
-  DO_MTHERR("hankel1e:");
+  DO_MTHERR("hankel1e:", &cy);
   if (sign == -1) {
     cy = rotate(cy, v);
   }
@@ -394,8 +492,7 @@
     sign = -1;
   }
   F_FUNC(zbesh,ZBESH)(CADDR(z), &v, &kode, &m, &n, CADDR(cy), &nz, &ierr);
-  F_FUNC(zbesh,ZBESH)(CADDR(z), &v,  &kode, &m, &n, CADDR(cy), &nz, &ierr);
-  DO_MTHERR("hankel2:");
+  DO_MTHERR("hankel2:", &cy);
   if (sign == -1) {
     cy = rotate(cy, -v);
   }
@@ -415,8 +512,7 @@
     sign = -1;
   }
   F_FUNC(zbesh,ZBESH)(CADDR(z), &v, &kode, &m, &n, CADDR(cy), &nz, &ierr);
-  F_FUNC(zbesh,ZBESH)(CADDR(z), &v, &kode, &m, &n, CADDR(cy), &nz, &ierr);
-  DO_MTHERR("hankel2e:");
+  DO_MTHERR("hankel2e:", &cy);
   if (sign == -1) {
     cy = rotate(cy, -v);
   }

Modified: trunk/scipy/special/amos_wrappers.h
===================================================================
--- trunk/scipy/special/amos_wrappers.h	2009-01-25 18:50:54 UTC (rev 5522)
+++ trunk/scipy/special/amos_wrappers.h	2009-01-25 19:02:17 UTC (rev 5523)
@@ -10,18 +10,32 @@
 #include "Python.h"
 #include "cephes/mconf.h"
 
-#define DO_MTHERR(name) if (nz !=0 || ierr !=0) mtherr(name, ierr_to_mtherr(nz, ierr))
-int ierr_to_mtherr( int nz, int ierr); 
+#define DO_MTHERR(name, varp)                         \
+    do {                                              \
+      if (nz !=0 || ierr != 0) {                      \
+        mtherr(name, ierr_to_mtherr(nz, ierr));       \
+        set_nan_if_no_computation_done(varp, ierr);   \
+      }                                               \
+    } while (0)
+
+int ierr_to_mtherr( int nz, int ierr);
+void set_nan_if_no_computation_done(Py_complex *var, int ierr);
 int cairy_wrap(Py_complex z, Py_complex *ai, Py_complex *aip, Py_complex *bi, Py_complex *bip);
 int cairy_wrap_e(Py_complex z, Py_complex *ai, Py_complex *aip, Py_complex *bi, Py_complex *bip);
+int cairy_wrap_e_real(double z, double *ai, double *aip, double *bi, double *bip);
 Py_complex cbesi_wrap( double v, Py_complex z);
 Py_complex cbesi_wrap_e( double v, Py_complex z);
+double cbesi_wrap_e_real( double v, double z);
 Py_complex cbesj_wrap( double v, Py_complex z);
 Py_complex cbesj_wrap_e( double v, Py_complex z);
+double cbesj_wrap_e_real( double v, double z);
 Py_complex cbesy_wrap( double v, Py_complex z);
 Py_complex cbesy_wrap_e( double v, Py_complex z);
+double cbesy_wrap_e_real( double v, double z);
 Py_complex cbesk_wrap( double v, Py_complex z);
 Py_complex cbesk_wrap_e( double v, Py_complex z);  
+double cbesk_wrap_real( double v, double z);
+double cbesk_wrap_e_real( double v, double z);
 Py_complex cbesh_wrap1( double v, Py_complex z);
 Py_complex cbesh_wrap1_e( double v, Py_complex z);  
 Py_complex cbesh_wrap2( double v, Py_complex z);

Modified: trunk/scipy/special/cephes/iv.c
===================================================================
--- trunk/scipy/special/cephes/iv.c	2009-01-25 18:50:54 UTC (rev 5522)
+++ trunk/scipy/special/cephes/iv.c	2009-01-25 19:02:17 UTC (rev 5523)
@@ -63,7 +63,7 @@
 #else
 double hyperg(), exp(), gamma(), log(), fabs(), floor();
 #endif
-extern double MACHEP, MAXNUM;
+extern double MACHEP, MAXNUM, NAN;
 
 double iv( v, x )
 double v, x;
@@ -88,7 +88,7 @@
 	if( t != v )
 		{
 		mtherr( "iv", DOMAIN );
-		return( 0.0 );
+		return( NAN );
 		}
 	if( v != 2.0 * floor(v/2.0) )
 		sign = -1;

Modified: trunk/scipy/special/cephes/jv.c
===================================================================
--- trunk/scipy/special/cephes/jv.c	2009-01-25 18:50:54 UTC (rev 5522)
+++ trunk/scipy/special/cephes/jv.c	2009-01-25 19:02:17 UTC (rev 5523)
@@ -90,7 +90,7 @@
 static double recur(), jvs(), hankel(), jnx(), jnt();
 #endif
 
-extern double MAXNUM, MACHEP, MINLOG, MAXLOG;
+extern double MAXNUM, MACHEP, MINLOG, MAXLOG, INFINITY, NAN;
 #define BIG  1.44115188075855872E+17
 
 double jv(double n, double x)
@@ -123,7 +123,7 @@
 
     if ((x < 0.0) && (y != an)) {
 	mtherr("Jv", DOMAIN);
-	y = 0.0;
+	y = NAN;
 	goto done;
     }
 
@@ -231,7 +231,7 @@
         */
 	if (n < 0.0) {
 	    mtherr("Jv", TLOSS);
-	    y = 0.0;
+	    y = NAN;
 	    goto done;
 	}
 	t = x / n;

Modified: trunk/scipy/special/tests/test_basic.py
===================================================================
--- trunk/scipy/special/tests/test_basic.py	2009-01-25 18:50:54 UTC (rev 5522)
+++ trunk/scipy/special/tests/test_basic.py	2009-01-25 19:02:17 UTC (rev 5523)
@@ -1629,21 +1629,32 @@
         yvp1 = yvp(2,.2)
         assert_array_almost_equal(yvp1,yvpr,10)
 
-class TestBesselJ(object):
-
-    def test_cephes_vs_amos(self):
+    def check_cephes_vs_amos(self, f1, f2, rtol=1e-11, atol=0):
         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)
+            for z in [-1300, -11, -10, -1, 1., 10., 200.5, 401., 600.5, 700.6,
+                      1300, 10003]:
+                c1, c2, c3 = f1(v, z), f1(v,z+0j), f2(int(v), z)
                 if np.isinf(c1):
-                    assert np.abs(c2) >= 1e150
+                    assert np.abs(c2) >= 1e150, (v, z)
                 elif np.isnan(c1):
-                    assert np.abs(c2.imag) > 1e-10
+                    assert c2.imag != 0, (v, z)
                 else:
-                    assert_tol_equal(c1, c2, err_msg=(v, z), rtol=1e-11)
+                    assert_tol_equal(c1, c2, err_msg=(v, z), rtol=rtol, atol=atol)
                     if v == int(v):
-                        assert_tol_equal(c2, c3, err_msg=(v, z), rtol=1e-11)
+                        assert_tol_equal(c3, c2, err_msg=(v, z), rtol=rtol, atol=atol)
 
+    def test_jv_cephes_vs_amos(self):
+        self.check_cephes_vs_amos(jv, jn, rtol=1e-10, atol=1e-305)
+
+    def test_yv_cephes_vs_amos(self):
+        self.check_cephes_vs_amos(yv, yn, rtol=1e-11, atol=1e-305)
+
+    def test_iv_cephes_vs_amos(self):
+        self.check_cephes_vs_amos(iv, iv, rtol=1e-8, atol=1e-305)
+
+    def test_kv_cephes_vs_amos(self):
+        self.check_cephes_vs_amos(kv, kn, rtol=1e-9, atol=1e-305)
+
     def test_ticket_623(self):
         assert_tol_equal(jv(3, 4), 0.43017147387562193)
         assert_tol_equal(jv(301, 1300), 0.0183487151115275)
@@ -1693,10 +1704,24 @@
         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))
 
+    def test_ticket_854(self):
+        """Real-valued Bessel domains"""
+        assert isnan(jv(0.5, -1))
+        assert isnan(iv(0.5, -1))
+        assert isnan(yv(0.5, -1))
+        assert isnan(yv(1, -1))
+        assert isnan(kv(0.5, -1))
+        assert isnan(kv(1, -1))
+        assert isnan(jve(0.5, -1))
+        assert isnan(ive(0.5, -1))
+        assert isnan(yve(0.5, -1))
+        assert isnan(yve(1, -1))
+        assert isnan(kve(0.5, -1))
+        assert isnan(kve(1, -1))
+        assert isnan(airye(-1)[0:2]).all(), airye(-1)
+        assert not isnan(airye(-1)[2:4]).any(), airye(-1)
 
-class TestBesselI(object):
-
-    def _series(self, v, z, n=200):
+    def iv_series(self, v, z, n=200):
         k = arange(0, n).astype(float_)
         r = (v+2*k)*log(.5*z) - gammaln(k+1) - gammaln(v+k+1)
         r[isnan(r)] = inf
@@ -1706,29 +1731,20 @@
 
     def test_i0_series(self):
         for z in [1., 10., 200.5]:
-            value, err = self._series(0, z)
+            value, err = self.iv_series(0, z)
             assert_tol_equal(i0(z), value, atol=err, err_msg=z)
 
     def test_i1_series(self):
             for z in [1., 10., 200.5]:
-                value, err = self._series(1, z)
+                value, err = self.iv_series(1, z)
                 assert_tol_equal(i1(z), value, atol=err, err_msg=z)
 
     def test_iv_series(self):
         for v in [-20., -10., -1., 0., 1., 12.49, 120.]:
             for z in [1., 10., 200.5, -1+2j]:
-                value, err = self._series(v, z)
+                value, err = self.iv_series(v, z)
                 assert_tol_equal(iv(v, z), value, atol=err, err_msg=(v, z))
 
-    def test_cephes_vs_specfun(self):
-        for v in [-120, -20., -10., -1., 0., 1., 12.49, 120.]:
-            for z in [1., 10., 200.5, 400., 600.5, 700.6]:
-                c1, c2 = iv(v, z), iv(v,z+0j)
-                if np.isinf(c1):
-                    assert np.abs(c2) >= 1e150
-                else:
-                    assert_tol_equal(c1, c2, err_msg=(v, z), rtol=1e-11)
-
     def test_i0(self):
         values = [[0.0, 1.0],
                   [1e-10, 1.0],



More information about the Scipy-svn mailing list