Skip to content

Drop workarounds for buggy libm's pow() #8

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 2 commits into from
Closed

Drop workarounds for buggy libm's pow() #8

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
81eab50
Drop workarounds for buggy libm's pow()
skirpichev May 2, 2025
7fb691e
Merge branch 'main' into no-pow-wa/133274
skirpichev May 2, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
58 changes: 8 additions & 50 deletions Modules/mathmodule.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2980,64 +2980,22 @@ static PyObject *
math_pow_impl(PyObject *module, double x, double y)
/*[clinic end generated code: output=fff93e65abccd6b0 input=c26f1f6075088bfd]*/
{
double r;
int odd_y;
errno = 0;

/* deal directly with IEEE specials, to cope with problems on various
platforms whose semantics don't exactly match C99 */
r = 0.; /* silence compiler warning */
if (!isfinite(x) || !isfinite(y)) {
errno = 0;
if (isnan(x))
r = y == 0. ? 1. : x; /* NaN**0 = 1 */
else if (isnan(y))
r = x == 1. ? 1. : y; /* 1**NaN = 1 */
else if (isinf(x)) {
odd_y = isfinite(y) && fmod(fabs(y), 2.0) == 1.0;
if (y > 0.)
r = odd_y ? x : fabs(x);
else if (y == 0.)
r = 1.;
else /* y < 0. */
r = odd_y ? copysign(0., x) : 0.;
}
else {
assert(isinf(y));
if (fabs(x) == 1.0)
r = 1.;
else if (y > 0. && fabs(x) > 1.0)
r = y;
else if (y < 0. && fabs(x) < 1.0) {
r = -y; /* result is +inf */
}
else
r = 0.;
}
}
else {
/* let libm handle finite**finite */
errno = 0;
r = pow(x, y);
double r = pow(x, y);

if (isfinite(x) && isfinite(y)) {
/* a NaN result should arise only from (-ve)**(finite
non-integer); in this case we want to raise ValueError. */
if (!isfinite(r)) {
if (isnan(r)) {
errno = EDOM;
}
/*
an infinite result here arises either from:
(A) (+/-0.)**negative (-> divide-by-zero)
(B) overflow of x**y with x and y finite
/* an infinite result here arises from:
(+/-0.)**negative (-> divide-by-zero)
*/
else if (isinf(r)) {
if (x == 0.)
errno = EDOM;
else
errno = ERANGE;
if (!x && isinf(r)) {
errno = EDOM;
}
}
}

if (errno && is_error(r, 1))
return NULL;
else
Expand Down
116 changes: 18 additions & 98 deletions Objects/floatobject.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -708,7 +708,6 @@ static PyObject *
float_pow(PyObject *v, PyObject *w, PyObject *z)
{
double iv, iw, ix;
int negate_result = 0;

if ((PyObject *)z != Py_None) {
PyErr_SetString(PyExc_TypeError, "pow() 3rd argument not "
Expand All @@ -719,105 +718,26 @@ float_pow(PyObject *v, PyObject *w, PyObject *z)
CONVERT_TO_DOUBLE(v, iv);
CONVERT_TO_DOUBLE(w, iw);

/* Sort out special cases here instead of relying on pow() */
if (iw == 0) { /* v**0 is 1, even 0**0 */
return PyFloat_FromDouble(1.0);
}
if (isnan(iv)) { /* nan**w = nan, unless w == 0 */
return PyFloat_FromDouble(iv);
}
if (isnan(iw)) { /* v**nan = nan, unless v == 1; 1**nan = 1 */
return PyFloat_FromDouble(iv == 1.0 ? 1.0 : iw);
}
if (isinf(iw)) {
/* v**inf is: 0.0 if abs(v) < 1; 1.0 if abs(v) == 1; inf if
* abs(v) > 1 (including case where v infinite)
*
* v**-inf is: inf if abs(v) < 1; 1.0 if abs(v) == 1; 0.0 if
* abs(v) > 1 (including case where v infinite)
*/
iv = fabs(iv);
if (iv == 1.0)
return PyFloat_FromDouble(1.0);
else if ((iw > 0.0) == (iv > 1.0))
return PyFloat_FromDouble(fabs(iw)); /* return inf */
else
return PyFloat_FromDouble(0.0);
}
if (isinf(iv)) {
/* (+-inf)**w is: inf for w positive, 0 for w negative; in
* both cases, we need to add the appropriate sign if w is
* an odd integer.
*/
int iw_is_odd = DOUBLE_IS_ODD_INTEGER(iw);
if (iw > 0.0)
return PyFloat_FromDouble(iw_is_odd ? iv : fabs(iv));
else
return PyFloat_FromDouble(iw_is_odd ?
copysign(0.0, iv) : 0.0);
}
if (iv == 0.0) { /* 0**w is: 0 for w positive, 1 for w zero
(already dealt with above), and an error
if w is negative. */
int iw_is_odd = DOUBLE_IS_ODD_INTEGER(iw);
if (iw < 0.0) {
PyErr_SetString(PyExc_ZeroDivisionError,
"zero to a negative power");
return NULL;
}
/* use correct sign if iw is odd */
return PyFloat_FromDouble(iw_is_odd ? iv : 0.0);
}

if (iv < 0.0) {
/* Whether this is an error is a mess, and bumps into libm
* bugs so we have to figure it out ourselves.
*/
if (iw != floor(iw)) {
/* Negative numbers raised to fractional powers
* become complex.
*/
return PyComplex_Type.tp_as_number->nb_power(v, w, z);
}
/* iw is an exact integer, albeit perhaps a very large
* one. Replace iv by its absolute value and remember
* to negate the pow result if iw is odd.
*/
iv = -iv;
negate_result = DOUBLE_IS_ODD_INTEGER(iw);
}

if (iv == 1.0) { /* 1**w is 1, even 1**inf and 1**nan */
/* (-1) ** large_integer also ends up here. Here's an
* extract from the comments for the previous
* implementation explaining why this special case is
* necessary:
*
* -1 raised to an exact integer should never be exceptional.
* Alas, some libms (chiefly glibc as of early 2003) return
* NaN and set EDOM on pow(-1, large_int) if the int doesn't
* happen to be representable in a *C* integer. That's a
* bug.
*/
return PyFloat_FromDouble(negate_result ? -1.0 : 1.0);
}

/* Now iv and iw are finite, iw is nonzero, and iv is
* positive and not equal to 1.0. We finally allow
* the platform pow to step in and do the rest.
*/
errno = 0;
ix = pow(iv, iw);
_Py_ADJUST_ERANGE1(ix);
if (negate_result)
ix = -ix;

if (errno != 0) {
/* We don't expect any errno value other than ERANGE, but
* the range of libm bugs appears unbounded.
*/
PyErr_SetFromErrno(errno == ERANGE ? PyExc_OverflowError :
PyExc_ValueError);
if (isfinite(iv) && isfinite(iw)) {
_Py_ADJUST_ERANGE1(ix);
}
if (errno) {
if (iv <= 0.0) {
if (iw != floor(iw)) {
/* Negative numbers raised to fractional powers
* become complex.
*/
return PyComplex_Type.tp_as_number->nb_power(v, w, z);
}
else {
PyErr_SetString(PyExc_ZeroDivisionError,
"zero to a negative power");
return NULL;
}
}
PyErr_SetFromErrno(PyExc_OverflowError);
return NULL;
}
return PyFloat_FromDouble(ix);
Expand Down
Loading