Skip to content

gh-156693: correct complex powers with infinite bases #156694

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

Open
skirpichev wants to merge 2 commits into
base: main
Choose a base branch
from
Open

gh-156693: correct complex powers with infinite bases #156694

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
42391fe
gh-156693: correct complex powers with infinite bases
skirpichev Aug 31, 2026
c93746a
Drop unused _Py_ADJUST_ERANGE2 macro
skirpichev Aug 31, 2026
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
3 changes: 3 additions & 0 deletions Include/internal/pycore_complexobject.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,9 @@ PyAPI_FUNC(Py_complex) _Py_cr_prod(Py_complex, double);
PyAPI_FUNC(Py_complex) _Py_cr_quot(Py_complex, double);
PyAPI_FUNC(Py_complex) _Py_rc_quot(double, Py_complex);

PyAPI_FUNC(bool) _Py_c_isnan(Py_complex);
PyAPI_FUNC(bool) _Py_c_isinf(Py_complex);
PyAPI_FUNC(bool) _Py_c_isfinite(Py_complex);

#ifdef __cplusplus
}
Expand Down
20 changes: 2 additions & 18 deletions Include/internal/pycore_pymath.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,8 @@ extern "C" {


/* _Py_ADJUST_ERANGE1(x)
* _Py_ADJUST_ERANGE2(x, y)
* Set errno to 0 before calling a libm function, and invoke one of these
* macros after, passing the function result(s) (_Py_ADJUST_ERANGE2 is useful
* for functions returning complex results). This makes two kinds of
* Set errno to 0 before calling a libm function, and invoke macro after.
* This makes two kinds of
* adjustments to errno: (A) If it looks like the platform libm set
* errno=ERANGE due to underflow, clear errno. (B) If it looks like the
* platform libm overflowed but didn't set errno, force errno to ERANGE. In
Expand Down Expand Up @@ -42,20 +40,6 @@ static inline void _Py_ADJUST_ERANGE1(double x)
}
}

static inline void _Py_ADJUST_ERANGE2(double x, double y)
{
if (x == INFINITY || x == -INFINITY ||
y == INFINITY || y == -INFINITY)
{
if (errno == 0) {
errno = ERANGE;
}
}
else if (errno == ERANGE) {
errno = 0;
}
}


//--- HAVE_PY_SET_53BIT_PRECISION macro ------------------------------------
//
Expand Down
8 changes: 8 additions & 0 deletions Lib/test/test_complex.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -362,6 +362,14 @@ def test_pow(self):
self.assertRaises(TypeError, pow, None, 1j)
self.assertAlmostEqual(pow(1j, 0.5), 0.7071067811865476+0.7071067811865475j)

# No overflows for infinite base
r = pow(complex(INF), 2) # integer power
self.assertEqual(r.real, INF)
self.assertTrue(isnan(r.imag))
r = pow(complex(INF), 2.25) # generic algorithm
self.assertEqual(r.real, INF)
self.assertTrue(isnan(r.imag))

a = 3.33+4.43j
self.assertEqual(a ** 0j, 1)
self.assertEqual(a ** 0.+0.j, 1)
Expand Down
3 changes: 3 additions & 0 deletions Misc/NEWS.d/next/Core_and_Builtins/2026-08-31-07-33-05.gh-issue-156693.BuDcDc.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
Correct complex powers of infinite numbers to not raise
:exc:`OverflowError`'s. As proposed in :gh:`156145`, we don't use platforms
:data:`errno` to manage exceptions. Patch by Sergey B Kirpichev.
24 changes: 12 additions & 12 deletions Modules/cmathmodule.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,7 @@ special_type(double d)
}

#define SPECIAL_VALUE(z, table) \
if (!isfinite((z).real) || !isfinite((z).imag)) { \
if (!_Py_c_isfinite(z)) { \
errno = 0; \
return table[special_type((z).real)] \
[special_type((z).imag)]; \
Expand Down Expand Up @@ -446,7 +446,7 @@ cmath_cosh_impl(PyObject *module, Py_complex z)
double x_minus_one;

/* special treatment for cosh(+/-inf + iy) if y is not a NaN */
if (!isfinite(z.real) || !isfinite(z.imag)) {
if (!_Py_c_isfinite(z)) {
if (isinf(z.real) && isfinite(z.imag) &&
(z.imag != 0.)) {
if (z.real > 0) {
Expand Down Expand Up @@ -482,7 +482,7 @@ cmath_cosh_impl(PyObject *module, Py_complex z)
r.imag = sin(z.imag) * sinh(z.real);
}
/* detect overflow, and set errno accordingly */
if (isinf(r.real) || isinf(r.imag))
if (_Py_c_isinf(r))
errno = ERANGE;
else
errno = 0;
Expand Down Expand Up @@ -515,7 +515,7 @@ cmath_exp_impl(PyObject *module, Py_complex z)
Py_complex r;
double l;

if (!isfinite(z.real) || !isfinite(z.imag)) {
if (!_Py_c_isfinite(z)) {
if (isinf(z.real) && isfinite(z.imag)
&& (z.imag != 0.)) {
if (z.real > 0) {
Expand Down Expand Up @@ -552,7 +552,7 @@ cmath_exp_impl(PyObject *module, Py_complex z)
r.imag = l*sin(z.imag);
}
/* detect overflow, and set errno accordingly */
if (isinf(r.real) || isinf(r.imag))
if (_Py_c_isinf(r))
errno = ERANGE;
else
errno = 0;
Expand Down Expand Up @@ -708,7 +708,7 @@ cmath_sinh_impl(PyObject *module, Py_complex z)

/* special treatment for sinh(+/-inf + iy) if y is finite and
nonzero */
if (!isfinite(z.real) || !isfinite(z.imag)) {
if (!_Py_c_isfinite(z)) {
if (isinf(z.real) && isfinite(z.imag)
&& (z.imag != 0.)) {
if (z.real > 0) {
Expand Down Expand Up @@ -742,7 +742,7 @@ cmath_sinh_impl(PyObject *module, Py_complex z)
r.imag = sin(z.imag) * cosh(z.real);
}
/* detect overflow, and set errno accordingly */
if (isinf(r.real) || isinf(r.imag))
if (_Py_c_isinf(r))
errno = ERANGE;
else
errno = 0;
Expand Down Expand Up @@ -894,7 +894,7 @@ cmath_tanh_impl(PyObject *module, Py_complex z)

/* special treatment for tanh(+/-inf + iy) if y is finite and
nonzero */
if (!isfinite(z.real) || !isfinite(z.imag)) {
if (!_Py_c_isfinite(z)) {
if (isinf(z.real) && isfinite(z.imag)
&& (z.imag != 0.)) {
if (z.real > 0) {
Expand Down Expand Up @@ -1129,7 +1129,7 @@ static PyObject *
cmath_isfinite_impl(PyObject *module, Py_complex z)
/*[clinic end generated code: output=ac76611e2c774a36 input=e224f5c36d94f5da]*/
{
return PyBool_FromLong(isfinite(z.real) && isfinite(z.imag));
return PyBool_FromLong(_Py_c_isfinite(z));
}

/*[clinic input]
Expand All @@ -1142,7 +1142,7 @@ static PyObject *
cmath_isnan_impl(PyObject *module, Py_complex z)
/*[clinic end generated code: output=e7abf6e0b28beab7 input=71799f5d284c9baf]*/
{
return PyBool_FromLong(isnan(z.real) || isnan(z.imag));
return PyBool_FromLong(_Py_c_isnan(z));
}

/*[clinic input]
Expand All @@ -1155,7 +1155,7 @@ static PyObject *
cmath_isinf_impl(PyObject *module, Py_complex z)
/*[clinic end generated code: output=502a75a79c773469 input=363df155c7181329]*/
{
return PyBool_FromLong(isinf(z.real) || isinf(z.imag));
return PyBool_FromLong(_Py_c_isinf(z));
}

/*[clinic input]
Expand Down Expand Up @@ -1211,7 +1211,7 @@ cmath_isclose_impl(PyObject *module, Py_complex a, Py_complex b,
above.
*/

if (isinf(a.real) || isinf(a.imag) || isinf(b.real) || isinf(b.imag)) {
if (_Py_c_isinf(a) || _Py_c_isinf(b)) {
return 0;
}

Expand Down
82 changes: 56 additions & 26 deletions Objects/complexobject.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,6 @@
#include "pycore_freelist.h" // _Py_FREELIST_FREE(), _Py_FREELIST_POP()
#include "pycore_long.h" // _PyLong_GetZero()
#include "pycore_object.h" // _PyObject_Init()
#include "pycore_pymath.h" // _Py_ADJUST_ERANGE2()


#define _PyComplexObject_CAST(op) ((PyComplexObject *)(op))
Expand All @@ -28,6 +27,24 @@ class complex "PyComplexObject *" "&PyComplex_Type"

static Py_complex c_1 = {1., 0.};

bool
_Py_c_isnan(Py_complex a)
{
return isnan(a.real) || isnan(a.imag);
}

bool
_Py_c_isinf(Py_complex a)
{
return isinf(a.real) || isinf(a.imag);
}

bool
_Py_c_isfinite(Py_complex a)
{
return isfinite(a.real) && isfinite(a.imag);
}

Py_complex
_Py_c_sum(Py_complex a, Py_complex b)
{
Expand Down Expand Up @@ -98,7 +115,7 @@ _Py_c_prod(Py_complex z, Py_complex w)
if (isnan(r.real) && isnan(r.imag)) {
int recalc = 0;

if (isinf(a) || isinf(b)) { /* z is infinite */
if (_Py_c_isinf(z)) {
/* "Box" the infinity and change nans in the other factor to 0 */
a = copysign(isinf(a) ? 1.0 : 0.0, a);
b = copysign(isinf(b) ? 1.0 : 0.0, b);
Expand All @@ -110,7 +127,7 @@ _Py_c_prod(Py_complex z, Py_complex w)
}
recalc = 1;
}
if (isinf(c) || isinf(d)) { /* w is infinite */
if (_Py_c_isinf(w)) {
/* "Box" the infinity and change nans in the other factor to 0 */
c = copysign(isinf(c) ? 1.0 : 0.0, c);
d = copysign(isinf(d) ? 1.0 : 0.0, d);
Expand Down Expand Up @@ -224,17 +241,13 @@ _Py_c_quot(Py_complex a, Py_complex b)
/* Recover infinities and zeros that computed as nan+nanj. See e.g.
the C11, Annex G.5.2, routine _Cdivd(). */
if (isnan(r.real) && isnan(r.imag)) {
if ((isinf(a.real) || isinf(a.imag))
&& isfinite(b.real) && isfinite(b.imag))
{
if (_Py_c_isinf(a) && _Py_c_isfinite(b)) {
const double x = copysign(isinf(a.real) ? 1.0 : 0.0, a.real);
const double y = copysign(isinf(a.imag) ? 1.0 : 0.0, a.imag);
r.real = INFINITY * (x*b.real + y*b.imag);
r.imag = INFINITY * (y*b.real - x*b.imag);
}
else if ((isinf(abs_breal) || isinf(abs_bimag))
&& isfinite(a.real) && isfinite(a.imag))
{
else if (_Py_c_isinf(b) && _Py_c_isfinite(a)) {
const double x = copysign(isinf(b.real) ? 1.0 : 0.0, b.real);
const double y = copysign(isinf(b.imag) ? 1.0 : 0.0, b.imag);
r.real = 0.0 * (a.real*x + a.imag*y);
Expand Down Expand Up @@ -291,9 +304,7 @@ _Py_rc_quot(double a, Py_complex b)
r.real = r.imag = Py_NAN;
}

if (isnan(r.real) && isnan(r.imag) && isfinite(a)
&& (isinf(abs_breal) || isinf(abs_bimag)))
{
if (isnan(r.real) && isnan(r.imag) && isfinite(a) && _Py_c_isinf(b)) {
const double x = copysign(isinf(b.real) ? 1.0 : 0.0, b.real);
const double y = copysign(isinf(b.imag) ? 1.0 : 0.0, b.imag);
r.real = 0.0 * (a*x);
Expand All @@ -306,20 +317,25 @@ _Py_rc_quot(double a, Py_complex b)
#pragma optimize("", on)
#endif

Py_complex
_Py_c_pow(Py_complex a, Py_complex b)
static Py_complex
c_pow(Py_complex a, Py_complex b)
{
Py_complex r;
double vabs,len,at,phase;

if (b.real == 0. && b.imag == 0.) {
r.real = 1.;
r.imag = 0.;
}
else if (a.real == 0. && a.imag == 0.) {
if (b.imag != 0. || b.real < 0.)
errno = EDOM;
r.real = 0.;
r.imag = 0.;
if (b.imag != 0. || b.real < 0.) {
r.real = NAN;
r.imag = NAN;
}
else {
r.real = 0.;
r.imag = 0.;
}
}
else {
vabs = hypot(a.real,a.imag);
Expand All @@ -332,8 +348,25 @@ _Py_c_pow(Py_complex a, Py_complex b)
}
r.real = len*cos(phase);
r.imag = len*sin(phase);
Comment on lines 341 to 350

@skirpichev skirpichev Aug 31, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

BTW, this generic code should handle special cases above correctly. Maybe we should drop them?

}
return r;
}

_Py_ADJUST_ERANGE2(r.real, r.imag);
Py_complex
_Py_c_pow(Py_complex a, Py_complex b)
{
int saved_errno = errno;
Py_complex r = c_pow(a, b);

if (_Py_c_isnan(r) && a.real == 0 && a.imag == 0) {
errno = EDOM;
r.real = r.imag = 0.0; /* and set r as documented */
}
else if (_Py_c_isinf(r) && _Py_c_isfinite(a) && _Py_c_isfinite(b)) {
errno = ERANGE;
}
else {
errno = saved_errno;
}
return r;
}
Expand Down Expand Up @@ -370,7 +403,7 @@ _Py_c_abs(Py_complex z)
/* sets errno = ERANGE on overflow; otherwise errno = 0 */
double result;

if (!isfinite(z.real) || !isfinite(z.imag)) {
if (!_Py_c_isfinite(z)) {
/* C99 rules: if either the real or the imaginary part is an
infinity, return infinity, even if the other part is a
NaN. */
Expand Down Expand Up @@ -748,23 +781,20 @@ complex_pow(PyObject *v, PyObject *w, PyObject *z)
PyErr_SetString(PyExc_ValueError, "complex modulo");
return NULL;
}
errno = 0;
// Check whether the exponent has a small integer value, and if so use
// a faster and more accurate algorithm.
if (b.imag == 0.0 && b.real == floor(b.real) && fabs(b.real) <= 100.0) {
p = c_powi(a, (long)b.real);
_Py_ADJUST_ERANGE2(p.real, p.imag);
}
else {
p = _Py_c_pow(a, b);
p = c_pow(a, b);
}

if (errno == EDOM) {
if (_Py_c_isnan(p) && a.real == 0 && a.imag == 0) {
PyErr_SetString(PyExc_ZeroDivisionError,
"zero to a negative or complex power");
return NULL;
}
else if (errno == ERANGE) {
else if (_Py_c_isinf(p) && _Py_c_isfinite(a) && _Py_c_isfinite(b)) {
PyErr_SetString(PyExc_OverflowError,
"complex exponentiation");
return NULL;
Expand Down
Loading