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

[scipy-svn at scipy.org]

Author: cookedm
Date: 2006-07-05 15:02:38 -0500 (Wed, 05 Jul 2006)
New Revision: 2042

Modified:
   trunk/Lib/special/cephes/hyp2f1.c
   trunk/Lib/special/tests/test_basic.py
Log:
cephes: fix up hyp2f1 and its test case

Modified: trunk/Lib/special/cephes/hyp2f1.c
===================================================================
--- trunk/Lib/special/cephes/hyp2f1.c	2006-07-05 18:27:34 UTC (rev 2041)
+++ trunk/Lib/special/cephes/hyp2f1.c	2006-07-05 20:02:38 UTC (rev 2042)
@@ -117,28 +117,37 @@
 double d, d1, d2, e;
 double p, q, r, s, y, ax;
 double ia, ib, ic, id, err;
-int flag, i, aid;
+int i, aid;
+int neg_int_a = 0, neg_int_b = 0;
+int neg_int_ca_or_cb = 0;
 
 err = 0.0;
 ax = fabs(x);
 s = 1.0 - x;
-flag = 0;
 ia = round(a); /* nearest integer to a */
 ib = round(b);
 
-if( a <= 0 )
-	{
-	if( fabs(a-ia) < EPS )		/* a is a negative integer */
-		flag |= 1;
-	}
+if (x == 0.0) {
+    return 1.0;
+}
 
-if( b <= 0 )
-	{
-	if( fabs(b-ib) < EPS )		/* b is a negative integer */
-		flag |= 2;
-	}
+d = c - a - b;
+if (d <= -1) {
+    return pow(s, d) * hyp2f1(c-a, c-b, c, x);
+}
+if (d <= 0 && x == 1)
+    goto hypdiv;
 
-if (fabs(ax-1.0) > EPS) {               /* |x| != 1.0 */
+if (a <= 0 && fabs(a-ia) < EPS ) {	/* a is a negative integer */
+    neg_int_a = 1;
+}
+
+if (b <= 0 && fabs(b-ib) < EPS ) {	/* b is a negative integer */
+    neg_int_b = 1;
+}
+
+if (ax < 1.0 || x == -1.0) {
+	/* 2F1(a,b;b;x) = (1-x)**(-a) */
 	if( fabs(b-c) < EPS ) {		/* b = c */
 		y = pow( s, -a );	/* s to the -a power */
 		goto hypdon;
@@ -157,24 +166,28 @@
 	if( fabs(c-ic) < EPS )		/* c is a negative integer */
 		{
 		/* check if termination before explosion */
-		if( (flag & 1) && (ia > ic) )
+		if( neg_int_a && (ia > ic) )
 			goto hypok;
-		if( (flag & 2) && (ib > ic) )
+		if( neg_int_b && (ib > ic) )
 			goto hypok;
 		goto hypdiv;
 		}
 	}
 
-if( flag )			/* function is a polynomial */
+if (neg_int_a || neg_int_b)		/* function is a polynomial */
 	goto hypok;
 
 if (x < -1.0) {
     double t1;
-    r = -x;
+    t1 = fabs(b - a);
+    if (fabs(t1 - round(t1)) < EPS) {
+        /* this transformation has a pole for b-a= +-integer */
+        goto hypdiv;
+    }
     p = hyp2f1(a, 1-c+a, 1-b+a, 1.0/x);
-    q = hyp2f1(b, 1-c+b, 1-b+a, 1.0/x);
-    p *= pow(r, -a);
-    q *= pow(r, -b);
+    q = hyp2f1(b, 1-c+b, 1-a+b, 1.0/x);
+    p *= pow(-x, -a);
+    q *= pow(-x, -b);
     t1 = gamma(c);
     s = t1*gamma(b-a)/(gamma(b)*gamma(c-a));
     y = t1*gamma(a-b)/(gamma(a)*gamma(c-b));
@@ -187,41 +200,35 @@
 p = c - a;
 ia = round(p); /* nearest integer to c-a */
 if( (ia <= 0.0) && (fabs(p-ia) < EPS) )	/* negative int c - a */
-	flag |= 4;
+        neg_int_ca_or_cb = 1;
 
 r = c - b;
 ib = round(r); /* nearest integer to c-b */
 if( (ib <= 0.0) && (fabs(r-ib) < EPS) )	/* negative int c - b */
-	flag |= 8;
+        neg_int_ca_or_cb = 1;
 
-d = c - a - b;
 id = round(d); /* nearest integer to d */
 q = fabs(d-id);
 
 /* Thanks to Christian Burger <BURGER at DMRHRZ11.HRZ.Uni-Marburg.DE>
  * for reporting a bug here.  */
-if( fabs(ax-1.0) < EPS )			/* |x| == 1.0	*/
-	{
-	if( x > 0.0 )
-		{
-		if( flag & 12 ) /* negative int c-a or c-b */
-			{
+if( fabs(ax-1.0) < EPS ) {			/* |x| == 1.0	*/
+	if( x > 0.0 ) {
+                if (neg_int_ca_or_cb) {
 			if( d >= 0.0 )
 				goto hypf;
 			else
 				goto hypdiv;
-			}
+		}
 		if( d <= 0.0 )
 			goto hypdiv;
 		y = gamma(c)*gamma(d)/(gamma(p)*gamma(r));
 		goto hypdon;
-		}
-
+	}
 	if( d <= -1.0 )
 		goto hypdiv;
+}
 
-	}
-
 /* Conditionally make d > 0 by recurrence on c
  * AMS55 #15.2.27
  */
@@ -250,7 +257,7 @@
 	}
 
 
-if( flag & 12 )
+if (neg_int_ca_or_cb)
 	goto hypf; /* negative integer c-a or c-b */
 
 hypok:

Modified: trunk/Lib/special/tests/test_basic.py
===================================================================
--- trunk/Lib/special/tests/test_basic.py	2006-07-05 18:27:34 UTC (rev 2041)
+++ trunk/Lib/special/tests/test_basic.py	2006-07-05 20:02:38 UTC (rev 2042)
@@ -1166,8 +1166,10 @@
                   [3, 4, 8, 1, gamma(8)*gamma(8-4-3)/gamma(8-3)/gamma(8-4)],
                   [3, 2, 3-2+1, -1, 1./2**3*sqrt(pi)*
                       gamma(1+3-2)/gamma(1+0.5*3-2)/gamma(0.5+0.5*3)],
-                  [4, 0.5+4, 5./6+4, 1./9, (0.75)**4*sqrt(pi)*
-                      gamma(5./6+2./3*4)/gamma(0.5+4./3)*gamma(5./6+4./3)],
+                  [5, 2, 5-2+1, -1, 1./2**5*sqrt(pi)*
+                      gamma(1+5-2)/gamma(1+0.5*5-2)/gamma(0.5+0.5*5)],
+                  [4, 0.5+4, 1.5-2*4, -1./3, (8./9)**(-2*4)*gamma(4./3)*
+                      gamma(1.5-2*4)/gamma(3./2)/gamma(4./3-2*4)],
                   ]
         for i, (a, b, c, x, v) in enumerate(values):
             cv = hyp2f1(a, b, c, x)