[Scipy-svn] r5523 - in trunk/scipy/special: . cephes tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Sun Jan 25 14:02:49 EST 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