[Scipy-svn] r5521 - in trunk/scipy/special: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Sun Jan 25 13:28:13 EST 2009
Author: ptvirtan
Date: 2009-01-25 12:28:00 -0600 (Sun, 25 Jan 2009)
New Revision: 5521
Modified:
trunk/scipy/special/amos_wrappers.c
trunk/scipy/special/tests/test_basic.py
Log:
Fixing #853: also I_v: use correct symmetry for negative orders (complex x)
Modified: trunk/scipy/special/amos_wrappers.c
===================================================================
--- trunk/scipy/special/amos_wrappers.c 2009-01-24 10:43:46 UTC (rev 5520)
+++ trunk/scipy/special/amos_wrappers.c 2009-01-25 18:28:00 UTC (rev 5521)
@@ -102,6 +102,24 @@
return 1;
}
+static int
+reflect_i(Py_complex *ik, double v)
+{
+ if (v != floor(v))
+ return 0;
+ return 1; /* I is symmetric for integer v */
+}
+
+static Py_complex
+rotate_i(Py_complex i, Py_complex k, double v)
+{
+ Py_complex w;
+ double s = sin(v * M_PI)*(2.0/M_PI);
+ w.real = i.real + s*k.real;
+ w.imag = i.imag + s*k.imag;
+ return w;
+}
+
int cairy_wrap(Py_complex z, Py_complex *ai, Py_complex *aip, Py_complex *bi, Py_complex *bip) {
int id = 0;
int ierr = 0;
@@ -142,28 +160,57 @@
Py_complex cbesi_wrap( double v, Py_complex z) {
int n = 1;
int kode = 1;
+ int sign = 1;
int nz, ierr;
- Py_complex cy;
+ Py_complex cy, cy_k;
if (v < 0) {
v = -v;
+ sign = -1;
}
F_FUNC(zbesi,ZBESI)(CADDR(z), &v, &kode, &n, CADDR(cy), &nz, &ierr);
DO_MTHERR("iv:");
+
+ 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):");
+ cy = rotate_i(cy, cy_k, v);
+ }
+ }
+
return cy;
}
Py_complex cbesi_wrap_e( double v, Py_complex z) {
int n = 1;
int kode = 2;
+ int sign = 1;
int nz, ierr;
- Py_complex cy;
+ Py_complex cy, cy_k;
if (v < 0) {
v = -v;
+ sign = -1;
}
F_FUNC(zbesi,ZBESI)(CADDR(z), &v, &kode, &n, CADDR(cy), &nz, &ierr);
DO_MTHERR("ive:");
+
+ 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):");
+ /* adjust scaling to match zbesi */
+ cy_k = rotate(cy_k, -z.imag/M_PI);
+ if (z.real > 0) {
+ cy_k.real *= exp(-2*z.real);
+ cy_k.imag *= exp(-2*z.real);
+ }
+ /* v -> -v */
+ cy = rotate_i(cy, cy_k, v);
+ }
+ }
+
return cy;
}
@@ -181,6 +228,7 @@
}
F_FUNC(zbesj,ZBESJ)(CADDR(z), &v, &kode, &n, CADDR(cy_j), &nz, &ierr);
DO_MTHERR("jv:");
+
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);
@@ -270,6 +318,7 @@
Py_complex cy;
if (v < 0) {
+ /* K_v == K_{-v} even for non-integer v */
v = -v;
}
F_FUNC(zbesk,ZBESK)(CADDR(z), &v, &kode, &n, CADDR(cy), &nz, &ierr);
@@ -284,6 +333,7 @@
Py_complex cy;
if (v < 0) {
+ /* K_v == K_{-v} even for non-integer v */
v = -v;
}
F_FUNC(zbesk,ZBESK)(CADDR(z), &v, &kode, &n, CADDR(cy), &nz, &ierr);
Modified: trunk/scipy/special/tests/test_basic.py
===================================================================
--- trunk/scipy/special/tests/test_basic.py 2009-01-24 10:43:46 UTC (rev 5520)
+++ trunk/scipy/special/tests/test_basic.py 2009-01-25 18:28:00 UTC (rev 5521)
@@ -1650,25 +1650,45 @@
assert_tol_equal(jv(301, 1296.0682), -0.0224174325312048)
def test_ticket_853(self):
+ """Negative-order Bessels"""
# 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(iv(-1, 1 ), 0.5651591039924851)
+ assert_tol_equal(iv(-2, 1 ), 0.1357476697670383)
+ assert_tol_equal(kv(-1, 1 ), 0.6019072301972347)
+ assert_tol_equal(kv(-2, 1 ), 1.624838898635178)
assert_tol_equal(jv(-0.5, 1 ), 0.43109886801837607952)
assert_tol_equal(yv(-0.5, 1 ), 0.6713967071418031)
+ assert_tol_equal(iv(-0.5, 1 ), 1.231200214592967)
+ assert_tol_equal(kv(-0.5, 1 ), 0.4610685044478945)
# 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(iv(-1, 1+0j), 0.5651591039924851)
+ assert_tol_equal(iv(-2, 1+0j), 0.1357476697670383)
+ assert_tol_equal(kv(-1, 1+0j), 0.6019072301972347)
+ assert_tol_equal(kv(-2, 1+0j), 1.624838898635178)
+
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(iv(-0.5, 1+0j), 1.231200214592967)
+ assert_tol_equal(iv(-0.5, 1+1j), 0.77070737376928+0.39891821043561j)
+ assert_tol_equal(kv(-0.5, 1+0j), 0.4610685044478945)
+ assert_tol_equal(kv(-0.5, 1+1j), 0.06868578341999-0.38157825981268j)
- 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(jve(-0.5,1+0.3j), jv(-0.5, 1+0.3j)*exp(-0.3))
+ assert_tol_equal(yve(-0.5,1+0.3j), yv(-0.5, 1+0.3j)*exp(-0.3))
+ assert_tol_equal(ive(-0.5,0.3+1j), iv(-0.5, 0.3+1j)*exp(-0.3))
+ assert_tol_equal(kve(-0.5,0.3+1j), kv(-0.5, 0.3+1j)*exp(0.3+1j))
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))
More information about the Scipy-svn
mailing list