diff Objects/floatobject.c @ 19457:6c2b7b624ef2 legacy-trunk

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).
author Tim Peters <tim.peters@gmail.com>
date Thu, 23 Aug 2001 22:31:37 +0000
parents 7257e87e0720
children 6c30d20d27f2
line wrap: on
line diff
--- a/Objects/floatobject.c	Thu Aug 23 21:40:38 2001 +0000
+++ b/Objects/floatobject.c	Thu Aug 23 22:31:37 2001 +0000
@@ -469,25 +469,10 @@
 	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 @@
   */
 	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 @@
 		}
 		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)
+	if (iv < 0.0 && iw != floor(iw)) {
+		PyErr_SetString(PyExc_ValueError,
+				"negative number cannot be raised to a fractional power");
+		return NULL;
 	}
-	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)
-	}
+	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 @@
 		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)
 	}