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).
This commit is contained in:
Tim Peters 2001-08-23 22:31:37 +00:00
parent 29a62dd6eb
commit 96685bfbf0
1 changed files with 16 additions and 45 deletions

View File

@ -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,20 +504,7 @@ 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) {
if (iv < 0.0 && iw != floor(iw)) {
PyErr_SetString(PyExc_ValueError,
"negative number cannot be raised to a fractional power");
return NULL;
@ -541,7 +513,6 @@ float_pow(PyObject *v, PyObject *w, PyObject *z)
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)
}