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

[scipy-svn at scipy.org]

Author: ptvirtan
Date: 2009-01-15 18:04:38 -0600 (Thu, 15 Jan 2009)
New Revision: 5466

Modified:
   trunk/scipy/special/specfun_wrappers.c
   trunk/scipy/special/tests/test_basic.py
Log:
Ensure that STVHV is never called with x<0

Use symmetries for integer v, and return NaN for non-integer v (when
the result would be complex). This works around bad accuracy in
Cephes and Specfun.

Modified: trunk/scipy/special/specfun_wrappers.c
===================================================================
--- trunk/scipy/special/specfun_wrappers.c	2009-01-16 00:04:07 UTC (rev 5465)
+++ trunk/scipy/special/specfun_wrappers.c	2009-01-16 00:04:38 UTC (rev 5466)
@@ -212,27 +212,40 @@
 
 double struve_wrap(double v, double x) {
   double out;
+  double rem;
   int flag=0;
 
+  if (x < 0) {
+      rem = fmod(v, 2.0);
+      if (rem == 0) {
+          x = -x;
+          flag = 1;
+      } else if (rem == 1 || rem == -1) {
+          x = -x;
+          flag = 0;
+      } else {
+          /* non-integer v and x < 0 => complex-valued */
+          return NAN;
+      }
+  }
+
   if ((v<-8.0) || (v>12.5)) {
-    return cephes_struve(v, x);  /* from cephes */
+    out = cephes_struve(v, x);  /* from cephes */
   }
-  if (v==0.0) {
-    if (x < 0) {x = -x; flag=1;}
+  else if (v==0.0) {
     F_FUNC(stvh0,STVH0)(&x,&out);
     CONVINF(out);
-    if (flag) out = -out;
-    return out;
   }
-  if (v==1.0) {
-    if (x < 0) x=-x;
+  else if (v==1.0) {
     F_FUNC(stvh1,STVH1)(&x,&out);
     CONVINF(out);
-    return out;
   }
-  F_FUNC(stvhv,STVHV)(&v,&x,&out);
-  CONVINF(out);
-  return out;  
+  else {
+    F_FUNC(stvhv,STVHV)(&v,&x,&out);
+    CONVINF(out);
+  }
+  if (flag) out = -out;
+  return out;
 }
 
 double modstruve_wrap(double v, double x) {

Modified: trunk/scipy/special/tests/test_basic.py
===================================================================
--- trunk/scipy/special/tests/test_basic.py	2009-01-16 00:04:07 UTC (rev 5465)
+++ trunk/scipy/special/tests/test_basic.py	2009-01-16 00:04:38 UTC (rev 5466)
@@ -2014,7 +2014,7 @@
 
     def test_vs_series(self):
         """Check Struve function versus its power series"""
-        for v in [-10, -7.99, -3.4, -1, 0, 1, 3.4, 12.49, 16]:
+        for v in [-20, -10, -7.99, -3.4, -1, 0, 1, 3.4, 12.49, 16]:
             for z in [1, 10, 19, 21, 30]:
                 value, err = self._series(v, z)
                 assert allclose(struve(v, z), value, atol=err), (v, z)
@@ -2023,7 +2023,15 @@
         assert_almost_equal(struve(-7.99, 21), 0.0467547614113, decimal=8)
         assert_almost_equal(struve(-8.01, 21), 0.0398716951023, decimal=9)
         assert_almost_equal(struve(-3.0, 200), 0.0142134427432, decimal=13)
+        assert_almost_equal(struve(-8.0, -41), 0.0192469727846, decimal=9)
+        assert_equal(struve(-12, -41), -struve(-12, 41))
+        assert_equal(struve(+12, -41), -struve(+12, 41))
+        assert_equal(struve(-11, -41), +struve(-11, 41))
+        assert_equal(struve(+11, -41), +struve(+11, 41))
 
+        assert isnan(struve(-7.1, -1))
+        assert isnan(struve(-10.1, -1))
+
     def test_regression_679(self):
         """Regression test for #679"""
         assert_almost_equal(struve(-1.0, 20 - 1e-8), struve(-1.0, 20 + 1e-8))