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

[scipy-svn at scipy.org]

Author: ptvirtan
Date: 2009-01-17 13:46:56 -0600 (Sat, 17 Jan 2009)
New Revision: 5478

Modified:
   trunk/scipy/special/cephes/hyperg.c
   trunk/scipy/special/cephes/iv.c
   trunk/scipy/special/tests/test_basic.py
Log:
Fixing #640: Make Cephes iv and hyperg behave better for large |z|

Bump iteration limits and handle infs and nans more carefully.

However, since iv is implemented in Cephes using 1F1, it overflows to inf
at smaller |z| than necessary, as intermediate results are larger than
floating point max. Fixing this would likely require a complete rewrite.

Modified: trunk/scipy/special/cephes/hyperg.c
===================================================================
--- trunk/scipy/special/cephes/hyperg.c	2009-01-17 19:45:37 UTC (rev 5477)
+++ trunk/scipy/special/cephes/hyperg.c	2009-01-17 19:46:56 UTC (rev 5478)
@@ -95,15 +95,22 @@
 	return( exp(x) * hyperg( temp, b, -x )  );
 
 
-psum = hy1f1p( a, b, x, &pcanc );
-if( pcanc < 1.0e-15 )
-	goto done;
+/* Try power & asymptotic series, starting from the one that is likely OK */
+if (fabs(x) < 10 + fabs(a) + fabs(b))
+        {
+        psum = hy1f1p( a, b, x, &pcanc );
+        if( pcanc < 1.0e-15 )
+ 	        goto done;
+        asum = hy1f1a( a, b, x, &acanc );
+        }
+else
+        {
+        psum = hy1f1a( a, b, x, &pcanc );
+        if( pcanc < 1.0e-15 )
+ 	        goto done;
+        asum = hy1f1p( a, b, x, &acanc );
+        }
 
-
-/* try asymptotic series */
-
-asum = hy1f1a( a, b, x, &acanc );
-
 /* Pick the result with less estimated error */
 
 if( acanc < pcanc )
@@ -129,7 +136,7 @@
 double a, b, x;
 double *err;
 {
-double n, a0, sum, t, u, temp;
+double n, a0, sum, t, u, temp, maxn;
 double an, bn, maxt;
 double y, c, sumc;
 
@@ -145,6 +152,7 @@
 maxt = 0.0;
 *err = 1.0;
 
+maxn = 200.0 + 2*fabs(a) + 2*fabs(b);
 
 while( t > MACHEP )
 	{
@@ -155,11 +163,11 @@
 		}
 	if( an == 0 )			/* a singularity		*/
 		return( sum );
-	if( n > 200 )
+	if( n > maxn )
                 {
-                /* did not converge: estimate 100% error */
-                *err = 1.0;
-                return sum;
+                /* too many terms; take the last one as error estimate */
+                c = fabs(c) + fabs(t)*50.0;
+		goto pdone;
                 }
 	u = x * ( an / (bn * n) );
 
@@ -194,6 +202,11 @@
 	*err = fabs(c);
 }
 
+if (*err != *err) {
+        /* nan */
+        *err = 1.0;
+}
+
 return( sum );
 }
 
@@ -261,7 +274,6 @@
 
 acanc = fabs(err1) + fabs(err2);
 
-
 if( b < 0 )
 	{
 	temp = gamma(b);
@@ -273,6 +285,10 @@
 if( asum != 0.0 )
 	acanc /= fabs(asum);
 
+if( asum*2 == asum )
+        /* infinity */
+        acanc = 0;
+
 acanc *= 30.0;	/* fudge factor, since error of asymptotic formula
 		 * often seems this much larger than advertised */
 
@@ -321,8 +337,8 @@
 	t = fabs(a0);
 
 	/* terminating condition for asymptotic series */
-	if( t > tlast )
-		goto ndone;
+	if( t < tlast && t < MACHEP*fabs(sum) )
+                goto ndone;
 
 	tlast = t;
 	sum += alast;	/* the sum is one term behind */
@@ -373,7 +389,6 @@
 /* estimate error due to roundoff, cancellation, and nonconvergence */
 *err = MACHEP * (n + maxt)  +  fabs ( a0 );
 
-
 done:
 sum += alast;
 return( sum );

Modified: trunk/scipy/special/cephes/iv.c
===================================================================
--- trunk/scipy/special/cephes/iv.c	2009-01-17 19:45:37 UTC (rev 5477)
+++ trunk/scipy/special/cephes/iv.c	2009-01-17 19:46:56 UTC (rev 5478)
@@ -69,7 +69,7 @@
 double v, x;
 {
 int sign;
-double t, ax;
+double t, ax, res;
 
 /* If v is a negative integer, invoke symmetry */
 t = floor(v);
@@ -112,5 +112,8 @@
 t = v * log( 0.5 * ax )  -  x;
 t = sign * exp(t) / gamma( v + 1.0 );
 ax = v + 0.5;
-return( t * hyperg( ax,  2.0 * ax,  2.0 * x ) );
+res = hyperg( ax,  2.0 * ax,  2.0 * x );
+if (res*2 == res)
+    return sign*res / gamma( v + 1.0 ); /* inf */
+return( t * res );
 }

Modified: trunk/scipy/special/tests/test_basic.py
===================================================================
--- trunk/scipy/special/tests/test_basic.py	2009-01-17 19:45:37 UTC (rev 5477)
+++ trunk/scipy/special/tests/test_basic.py	2009-01-17 19:46:56 UTC (rev 5478)
@@ -1591,9 +1591,15 @@
 
 class TestBesselI(object):
 
+    def _lggamma(self, q):
+        res = zeros_like(q)
+        res[q>=2] = gammaln(q[q>=2])
+        res[q<2] = log(gamma(q[q<2]))
+        return res
+
     def _series(self, v, z, n=200):
         k = arange(0, n).astype(float_)
-        r = (v+2*k)*log(.5*z) - log(gamma(k+1)) - log(gamma(v+k+1))
+        r = (v+2*k)*log(.5*z) - self._lggamma(k+1) - self._lggamma(v+k+1)
         r[isnan(r)] = inf
         r = exp(r)
         err = abs(r).max() * finfo(float_).eps * n + abs(r[-1])*10
@@ -1615,6 +1621,15 @@
                 value, err = self._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],