[Scipy-svn] r6949 - in trunk/scipy/special: cephes specfun tests

[scipy-svn at scipy.org]

Author: ptvirtan
Date: 2010-11-27 12:40:26 -0600 (Sat, 27 Nov 2010)
New Revision: 6949

Modified:
   trunk/scipy/special/cephes/hyp2f1.c
   trunk/scipy/special/specfun/specfun.f
   trunk/scipy/special/tests/test_orthogonal_eval.py
Log:
BUG: special: avoid raising spurious div-by-zero flags in hyp2f1 and hyp1f1 (fixes #1334)

Modified: trunk/scipy/special/cephes/hyp2f1.c
===================================================================
--- trunk/scipy/special/cephes/hyp2f1.c	2010-11-27 08:10:46 UTC (rev 6948)
+++ trunk/scipy/special/cephes/hyp2f1.c	2010-11-27 18:40:26 UTC (rev 6949)
@@ -381,10 +381,14 @@
                 p *= s * (a + t + d1) / (t + 1.0);
                 p *= (b + t + d1) / (t + 1.0 + e);
                 t += 1.0;
+                if (t > 10000) {      /* should never happen */
+                    mtherr("hyp2f1", TOOMANY);
+                    *loss = 1.0;
+                    return NPY_NAN;
+                }
             }
-            while (fabs(q / y) > EPS);
+            while (y == 0 || fabs(q / y) > EPS);
 
-
             if (id == 0.0) {
                 y *= gamma(c) / (gamma(a) * gamma(b));
                 goto psidon;
@@ -498,7 +502,7 @@
             return (s);
         }
     }
-    while (fabs(u / s) > MACHEP);
+    while (s == 0 || fabs(u / s) > MACHEP);
 
     /* return estimated relative error */
     *loss = (MACHEP * umax) / fabs(s) + (MACHEP * i);

Modified: trunk/scipy/special/specfun/specfun.f
===================================================================
--- trunk/scipy/special/specfun/specfun.f	2010-11-27 08:10:46 UTC (rev 6948)
+++ trunk/scipy/special/specfun/specfun.f	2010-11-27 18:40:26 UTC (rev 6949)
@@ -5675,7 +5675,7 @@
               DO 15 J=1,500
                  RG=RG*(A+J-1.0D0)/(J*(B+J-1.0D0))*X
                  HG=HG+RG
-                 IF (DABS(RG/HG).LT.1.0D-15) GO TO 25
+                 IF (HG.NE.0D0.AND.DABS(RG/HG).LT.1.0D-15) GO TO 25
 15            CONTINUE
            ELSE
               CALL GAMMA2(A,TA)

Modified: trunk/scipy/special/tests/test_orthogonal_eval.py
===================================================================
--- trunk/scipy/special/tests/test_orthogonal_eval.py	2010-11-27 08:10:46 UTC (rev 6948)
+++ trunk/scipy/special/tests/test_orthogonal_eval.py	2010-11-27 18:40:26 UTC (rev 6949)
@@ -1,5 +1,3 @@
-#from numpy.testing import 
-
 import numpy as np
 from numpy.testing import assert_
 import scipy.special.orthogonal as orth
@@ -15,6 +13,18 @@
     assert_(np.allclose(v1, v2, rtol=1e-15))
 
 
+def test_warnings():
+    # ticket 1334
+    olderr = np.seterr(all='raise')
+    try:
+        # these should raise no fp warnings
+        orth.eval_legendre(1, 0)
+        orth.eval_laguerre(1, 1)
+        orth.eval_gegenbauer(1, 1, 0)
+    finally:
+        np.seterr(**olderr)
+
+
 class TestPolys(object):
     """
     Check that the eval_* functions agree with the constructed polynomials
@@ -47,9 +57,13 @@
             p = (p[0].astype(int),) + p[1:]
             return func(*p)
 
-        ds = FuncData(polyfunc, dataset, range(len(param_ranges)+2), -1,
-                      rtol=rtol)
-        ds.check()
+        olderr = np.seterr(all='raise')
+        try:
+            ds = FuncData(polyfunc, dataset, range(len(param_ranges)+2), -1,
+                          rtol=rtol)
+            ds.check()
+        finally:
+            np.seterr(**olderr)
 
     def test_jacobi(self):
         self.check_poly(orth.eval_jacobi, orth.jacobi,