float_pow: Put *all* of the burden on the libm pow in normal

cases.
powu:  Deleted.

This started with a nonsensical error msg:

>>> x = -1.
>>> import sys
>>> x**(-sys.maxint-1L)
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
ValueError: negative number cannot be raised to a fractional power
>>>

The special-casing in float_pow was simply wrong in this case (there's
not even anything peculiar about these inputs), and I don't see any point
to it in *any* case:  a decent libm pow should have worst-case error under
1 ULP, so in particular should deliver the exact result whenever the exact
result is representable (else its error is at least 1 ULP).  Thus our
special fiddling for integral values "shouldn't" buy anything in accuracy,
and, to the contrary, repeated multiplication is less accurate than a
decent pow when the true result isn't exactly representable.  So just
letting pow() do its job here (we may not be able to trust libm x-platform
in exceptional cases, but these are normal cases).

1 files changed, 16 insertions, 45 deletions

diff --git a/Objects/floatobject.c b/Objects/floatobject.c
index 044d1d3..69aede4 100644
--- a/Objects/floatobject.c
+++ b/Objects/floatobject.c
@@ -469,25 +469,10 @@ float_divmod(PyObject *v, PyObject *w)
 	return Py_BuildValue("(dd)", floordiv, mod);
 }
 
-static double powu(double x, long n)
-{
-	double r = 1.;
-	double p = x;
-	long mask = 1;
-	while (mask > 0 && n >= mask) {
-		if (n & mask)
-			r *= p;
-		mask <<= 1;
-		p *= p;
-	}
-	return r;
-}
-
 static PyObject *
 float_pow(PyObject *v, PyObject *w, PyObject *z)
 {
 	double iv, iw, ix;
-	long intw;
  /* XXX Doesn't handle overflows if z!=None yet; it may never do so :(
   * The z parameter is really only going to be useful for integers and
   * long integers.  Maybe something clever with logarithms could be done.
@@ -495,23 +480,23 @@ float_pow(PyObject *v, PyObject *w, PyObject *z)
   */
 	CONVERT_TO_DOUBLE(v, iv);
 	CONVERT_TO_DOUBLE(w, iw);
-	intw = (long)iw;
 
 	/* Sort out special cases here instead of relying on pow() */
-	if (iw == 0) { 		/* x**0 is 1, even 0**0 */
+	if (iw == 0) { 		/* v**0 is 1, even 0**0 */
 		PyFPE_START_PROTECT("pow", return NULL)
 		if ((PyObject *)z != Py_None) {
 			double iz;
 			CONVERT_TO_DOUBLE(z, iz);
-			ix=fmod(1.0, iz);
-			if (ix!=0 && iz<0) ix+=iz;
+			ix = fmod(1.0, iz);
+			if (ix != 0 && iz < 0)
+				ix += iz;
 		}
 		else
 			ix = 1.0;
 		PyFPE_END_PROTECT(ix)
 		return PyFloat_FromDouble(ix); 
 	}
-	if (iv == 0.0) {
+	if (iv == 0.0) {  /* 0**w is error if w<0, else 1 */
 		if (iw < 0.0) {
 			PyErr_SetString(PyExc_ZeroDivisionError,
 					"0.0 cannot be raised to a negative power");
@@ -519,29 +504,15 @@ float_pow(PyObject *v, PyObject *w, PyObject *z)
 		}
 		return PyFloat_FromDouble(0.0);
 	}
-
-	if (iw == intw && intw > LONG_MIN) {
-		/* ruled out LONG_MIN because -LONG_MIN isn't representable */
-		errno = 0;
-		PyFPE_START_PROTECT("pow", return NULL)
-		if (intw > 0)
-			ix = powu(iv, intw);
-		else
-			ix = 1./powu(iv, -intw);
-		PyFPE_END_PROTECT(ix)
-	}
-	else {
-		/* Sort out special cases here instead of relying on pow() */
-		if (iv < 0.0) {
-			PyErr_SetString(PyExc_ValueError,
-					"negative number cannot be raised to a fractional power");
-			return NULL;
-		}
-		errno = 0;
-		PyFPE_START_PROTECT("pow", return NULL)
-		ix = pow(iv, iw);
-		PyFPE_END_PROTECT(ix)
+	if (iv < 0.0 && iw != floor(iw)) {
+		PyErr_SetString(PyExc_ValueError,
+				"negative number cannot be raised to a fractional power");
+		return NULL;
 	}
+	errno = 0;
+	PyFPE_START_PROTECT("pow", return NULL)
+	ix = pow(iv, iw);
+	PyFPE_END_PROTECT(ix)
 	CHECK(ix);
 	if (errno != 0) {
 		/* XXX could it be another type of error? */
@@ -552,9 +523,9 @@ float_pow(PyObject *v, PyObject *w, PyObject *z)
 		double iz;
 		CONVERT_TO_DOUBLE(z, iz);
 		PyFPE_START_PROTECT("pow", return 0)
-	 	ix=fmod(ix, iz);	/* XXX To Be Rewritten */
-	 	if (ix!=0 && ((iv<0 && iz>0) || (iv>0 && iz<0) )) {
-		     ix+=iz;
+	 	ix = fmod(ix, iz);	/* XXX To Be Rewritten */
+	 	if (ix != 0 && ((iv < 0 && iz > 0) || (iv > 0 && iz < 0) )) {
+		     ix += iz;
 		}
   		PyFPE_END_PROTECT(ix)
 	}