[Scipy-svn] r3012 - in trunk/Lib/special: cephes tests

[scipy-svn at scipy.org]

Author: cookedm
Date: 2007-05-16 15:01:46 -0500 (Wed, 16 May 2007)
New Revision: 3012

Modified:
   trunk/Lib/special/cephes/hyperg.c
   trunk/Lib/special/tests/test_basic.py
Log:
Lib/special/cephes/hyperg.c: Use Kahan summation in power series sum for hyp1f1, instead of ad-hoc calculation of roundoff/cancellation error.


Modified: trunk/Lib/special/cephes/hyperg.c
===================================================================
--- trunk/Lib/special/cephes/hyperg.c	2007-05-16 17:13:23 UTC (rev 3011)
+++ trunk/Lib/special/cephes/hyperg.c	2007-05-16 20:01:46 UTC (rev 3012)
@@ -132,6 +132,7 @@
 {
 double n, a0, sum, t, u, temp;
 double an, bn, maxt;
+double y, c, sumc;
 
 
 /* set up for power series summation */
@@ -139,6 +140,7 @@
 bn = b;
 a0 = 1.0;
 sum = 1.0;
+c = 0.0;
 n = 1.0;
 t = 1.0;
 maxt = 0.0;
@@ -167,17 +169,14 @@
 		}
 
 	a0 *= u;
-	sum += a0;
+
+	y = a0 - c;
+	sumc = sum + y;
+	c = (sumc - sum) - y;
+	sum = sumc;
+
 	t = fabs(a0);
-	if( t > maxt )
-		maxt = t;
-/*
-	if( (maxt/fabs(sum)) > 1.0e17 )
-		{
-		pcanc = 1.0;
-		goto blowup;
-		}
-*/
+
 	an += 1.0;
 	bn += 1.0;
 	n += 1.0;
@@ -186,10 +185,11 @@
 pdone:
 
 /* estimate error due to roundoff and cancellation */
-if( sum != 0.0 )
-	maxt /= fabs(sum);
-maxt *= MACHEP; 	/* this way avoids multiply overflow */
-*err = fabs( MACHEP * n  +  maxt );
+if (sum != 0.0) {
+	*err = fabs(c / sum);
+} else {
+	*err = fabs(c);
+}
 
 return( sum );
 }

Modified: trunk/Lib/special/tests/test_basic.py
===================================================================
--- trunk/Lib/special/tests/test_basic.py	2007-05-16 17:13:23 UTC (rev 3011)
+++ trunk/Lib/special/tests/test_basic.py	2007-05-16 20:01:46 UTC (rev 3012)
@@ -193,6 +193,8 @@
         cephes.hankel2e(1,1)
 
     def check_hyp1f1(self):
+        assert_approx_equal(cephes.hyp1f1(1,1,1), exp(1.0))
+        assert_approx_equal(cephes.hyp1f1(3,4,-6), 0.026056422099537251095)
         cephes.hyp1f1(1,1,1)
     def check_hyp1f2(self):
         cephes.hyp1f2(1,1,1,1)