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).

diff --git a/Objects/floatobject.c b/Objects/floatobject.c
index 044d1d3af8c2a126ee90d4366e42fe52cd8b26e9..69aede4cc1950417c2ef334278487b1b56fcf027 100644 (file)
--- 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)
        }