From 96685bfbf0e954febd4f873721d3db5d03cb2d99 Mon Sep 17 00:00:00 2001 From: Tim Peters Date: Thu, 23 Aug 2001 22:31:37 +0000 Subject: [PATCH] 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 "", 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). --- Objects/floatobject.c | 61 ++++++++++++------------------------------- 1 file changed, 16 insertions(+), 45 deletions(-) diff --git a/Objects/floatobject.c b/Objects/floatobject.c index 044d1d3af8c..69aede4cc19 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) }