Skip to content

gh-117999: Fix small integer powers of complex numbers #124243

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
serhiy-storchaka wants to merge 7 commits into
base: main
Choose a base branch
from
+88 −14
Open

gh-117999: Fix small integer powers of complex numbers #124243

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
f2280be
gh-117999: Fix small integer powers of complex numbers
serhiy-storchaka Sep 19, 2024
5b5320a
Merge branch 'main' into complex_powi
serhiy-storchaka Nov 26, 2024
1304423
Polish tests.
serhiy-storchaka Nov 26, 2024
bc4adfd
Remove unrelated change.
serhiy-storchaka Nov 26, 2024
4b8535d
Merge branch 'main' into complex_powi
serhiy-storchaka Dec 12, 2024
e876273
Add more tests.
serhiy-storchaka Dec 12, 2024
bb748c7
Merge branch 'main' into complex_powi
serhiy-storchaka Feb 4, 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
Next Next commit
gh-117999: Fix small integer powers of complex numbers 
* Fix the sign of zero components in the result. E.g. complex(1,-0.0)**2
  now evaluates to complex(1,-0.0) instead of complex(1,-0.0).
* Fix negative small integer powers of infinite complex numbers. E.g.
  complex(inf)**-1 now evaluates to complex(0,-0.0) instead of
  complex(nan,nan).
* Powers of infinite numbers no longer raise OverflowError. E.g.
  complex(inf)**1 now evaluates to complex(inf) and complex(inf)**0.5
  now evaluates to complex(inf,nan).
  • Loading branch information
@serhiy-storchaka
serhiy-storchaka committed Sep 19, 2024
commit f2280bec1b1863a5303b981537ec47f8d269bfb0
29 changes: 29 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 @@ -366,6 +366,35 @@ def test_pow_with_small_integer_exponents(self):
self.assertEqual(str(float_pow), str(int_pow))
self.assertEqual(str(complex_pow), str(int_pow))

# Check that complex numbers with special components
# are correctly handled.
for x in [2, -2, +0.0, -0.0, INF, -INF, NAN]:
for y in [2, -2, +0.0, -0.0, INF, -INF, NAN]:
c = complex(x, y)
with self.subTest(c):
self.assertComplexesAreIdentical(c**0, complex(1, +0.0))
self.assertComplexesAreIdentical(c**1, c)
self.assertComplexesAreIdentical(c**2, c*c)
self.assertComplexesAreIdentical(c**3, c*(c*c))
self.assertComplexesAreIdentical(c**4, (c*c)*(c*c))
self.assertComplexesAreIdentical(c**5, c*((c*c)*(c*c)))
Comment thread
skirpichev marked this conversation as resolved.
Outdated
for x in [+2, -2]:
for y in [+0.0, -0.0]:
with self.subTest(complex(x, y)):
self.assertComplexesAreIdentical(complex(x, y)**-1, complex(1/x, -y))
with self.subTest(complex(y, x)):
self.assertComplexesAreIdentical(complex(y, x)**-1, complex(y, -1/x))
for x in [+INF, -INF]:
for y in [+1, -1]:
c = complex(x, y)
with self.subTest(c):
self.assertComplexesAreIdentical(c**-1, complex(1/x, -0.0*y))
self.assertComplexesAreIdentical(c**-2, complex(0.0, -y/x))
c = complex(y, x)
with self.subTest(c):
self.assertComplexesAreIdentical(c**-1, complex(+0.0*y, -1/x))
self.assertComplexesAreIdentical(c**-2, complex(-0.0, -y/x))

def test_boolcontext(self):
for i in range(100):
self.assertTrue(complex(random() + 1e-6, random() + 1e-6))
Expand Down
2 changes: 2 additions & 0 deletions Misc/NEWS.d/next/Core_and_Builtins/2024-09-19-15-47-50.gh-issue-117999.Iq4jEG.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
Fix calculation of powers of complex numbers. Small integer powers now produce correct sign of zero components. Negative powers of infinite numbers now evaluate to zero instead of NaN.
Powers of infinite numbers no longer raise OverflowError.
101 changes: 85 additions & 16 deletions Objects/complexobject.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,10 +21,10 @@ class complex "PyComplexObject *" "&PyComplex_Type"

#include "clinic/complexobject.c.h"

/* elementary operations on complex numbers */

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

/* elementary operations on complex numbers */

Py_complex
_Py_c_sum(Py_complex a, Py_complex b)
{
Expand Down Expand Up @@ -143,6 +143,64 @@ _Py_c_quot(Py_complex a, Py_complex b)

return r;
}
Py_complex
_Py_dc_quot(double a, Py_complex b)
{
/* Divide a real number by a complex number.
* Same as _Py_c_quot(), but without imaginary part.
*/
Py_complex r; /* the result */
const double abs_breal = b.real < 0 ? -b.real : b.real;
const double abs_bimag = b.imag < 0 ? -b.imag : b.imag;

if (abs_breal >= abs_bimag) {
/* divide tops and bottom by b.real */
if (abs_breal == 0.0) {
errno = EDOM;
r.real = r.imag = 0.0;
}
else {
const double ratio = b.imag / b.real;
const double denom = b.real + b.imag * ratio;
r.real = a / denom;
r.imag = (- a * ratio) / denom;
}
}
else if (abs_bimag >= abs_breal) {
/* divide tops and bottom by b.imag */
const double ratio = b.real / b.imag;
const double denom = b.real * ratio + b.imag;
assert(b.imag != 0.0);
r.real = (a * ratio) / denom;
r.imag = (-a) / denom;
}
else {
/* At least one of b.real or b.imag is a NaN */
r.real = r.imag = Py_NAN;
}

/* 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)
&& isfinite(b.real) && isfinite(b.imag))
{
const double x = copysign(isinf(a) ? 1.0 : 0.0, a);
r.real = Py_INFINITY * (x*b.real);
r.imag = Py_INFINITY * (- x*b.imag);
}
else if ((isinf(abs_breal) || isinf(abs_bimag))
&& 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*x);
r.imag = 0.0 * (-a*y);
}
}

return r;
}
#ifdef _M_ARM64
#pragma optimize("", on)
#endif
Expand Down Expand Up @@ -174,23 +232,31 @@ _Py_c_pow(Py_complex a, Py_complex b)
r.real = len*cos(phase);
r.imag = len*sin(phase);

_Py_ADJUST_ERANGE2(r.real, r.imag);
if (isfinite(a.real) && isfinite(a.imag)
&& isfinite(b.real) && isfinite(b.imag))
{
_Py_ADJUST_ERANGE2(r.real, r.imag);
}
Comment thread
skirpichev marked this conversation as resolved.
}
return r;
}

static Py_complex
c_powu(Py_complex x, long n)
{
Py_complex r, p;
long mask = 1;
r = c_1;
p = x;
while (mask > 0 && n >= mask) {
if (n & mask)
r = _Py_c_prod(r,p);
mask <<= 1;
p = _Py_c_prod(p,p);
assert(n > 0);
Copy link
Copy Markdown
Member

@picnixz picnixz Aug 18, 2025

Choose a reason for hiding this comment

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

If you want to make this function a bit more natural, I would suggest having if (n == 0) { return c_1; } and remove the assertion n > 0 and just assume n is an unsigned long (or change the assertion to n >= 0). I don't think it's necessary to handle the case n == 0 in c_powi().

Copy link
Copy Markdown
Member

@skirpichev skirpichev Aug 19, 2025

Choose a reason for hiding this comment

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

See my patch above: #124243 (comment)

Copy link
Copy Markdown
Member

@picnixz picnixz Aug 19, 2025

Choose a reason for hiding this comment

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

Oh I remember this one. I think we can just keep special-casing n == 0 without actually wondering about the cutoff limit. A generic square-and-multiply algorithm is fine IMO.

Copy link
Copy Markdown
Member

@skirpichev skirpichev Aug 19, 2025

Choose a reason for hiding this comment

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

Well, it's possible to find base and exponent, when binary exponentiation is worse.

E.g. on this branch, if we omit cutoff, on my system I have:

$ ./python -m timeit -u nsec -s 'z=1e-10+1j;e=1e18+0j' 'z**e'
200000 loops, best of 5: 1.39e+03 nsec per loop
$ ./python -m timeit -u nsec -s 'z=1e-10+1j;e=1e18+0j;from _testcapi import _py_c_pow as p' 'p(z, e)'
200000 loops, best of 5: 1.12e+03 nsec per loop

(Though, both results looks to be a garbage.)

while ((n & 1) == 0) {
x = _Py_c_prod(x, x);
n >>= 1;
}
Py_complex r = x;
n >>= 1;
while (n) {
x = _Py_c_prod(x, x);
if (n & 1) {
r = _Py_c_prod(r, x);
}
n >>= 1;
Comment thread
skirpichev marked this conversation as resolved.
Comment on lines +348 to +360
Copy link
Copy Markdown
Member

@skirpichev skirpichev Aug 13, 2025

Choose a reason for hiding this comment

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

For me, your tests works fine with my patch from #118000. Here it is, with slightly adjusted formatting (and using new _Py_rc_quot):

diff --git a/Objects/complexobject.c b/Objects/complexobject.c
index b66ebe131a..9dc5d22b37 100644
--- a/Objects/complexobject.c
+++ b/Objects/complexobject.c
@@ -333,19 +333,26 @@ _Py_c_pow(Py_complex a, Py_complex b)
         r.real = len*cos(phase);
         r.imag = len*sin(phase);
 
-        _Py_ADJUST_ERANGE2(r.real, r.imag);
+        if (isfinite(a.real) && isfinite(a.imag)
+            && isfinite(b.real) && isfinite(b.imag))
+        {
+            _Py_ADJUST_ERANGE2(r.real, r.imag);
+        }
     }
     return r;
 }
 
+#define INT_EXP_CUTOFF 100
+
 static Py_complex
 c_powu(Py_complex x, long n)
 {
-    Py_complex r, p;
+    Py_complex p = x, r = n-- ? p : c_1;
     long mask = 1;
-    r = c_1;
-    p = x;
-    while (mask > 0 && n >= mask) {
+    assert(-1 <= n);
+    assert(n < INT_EXP_CUTOFF);
+    while (n >= mask) {
+        assert(mask > 0);
         if (n & mask)
             r = _Py_c_prod(r,p);
         mask <<= 1;
@@ -357,11 +364,10 @@ c_powu(Py_complex x, long n)
 static Py_complex
 c_powi(Py_complex x, long n)
 {
-    if (n > 0)
-        return c_powu(x,n);
+    if (n >= 0)
+        return c_powu(x, n);
     else
-        return _Py_c_quot(c_1, c_powu(x,-n));
-
+        return _Py_rc_quot(1.0, c_powu(x, -n));
 }
 
 double
@@ -751,9 +757,13 @@ complex_pow(PyObject *v, PyObject *w, PyObject *z)
     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) {
+    if (b.imag == 0.0 && b.real == floor(b.real)
+        && fabs(b.real) <= INT_EXP_CUTOFF)
+    {
         p = c_powi(a, (long)b.real);
-        _Py_ADJUST_ERANGE2(p.real, p.imag);
+        if (isfinite(a.real) && isfinite(a.imag)) {
+            _Py_ADJUST_ERANGE2(p.real, p.imag);
+        }
     }
     else {
         p = _Py_c_pow(a, b);

(asserts are outcome of review, I'm not sure they are useful here.)

What do you think?

}
return r;
}
Expand All @@ -199,10 +265,11 @@ static Py_complex
c_powi(Py_complex x, long n)
{
if (n > 0)
return c_powu(x,n);
return c_powu(x, n);
else if (n < 0)
return _Py_dc_quot(1.0, c_powu(x, -n));
else
return _Py_c_quot(c_1, c_powu(x,-n));

return c_1;
Copy link
Copy Markdown
Member

@picnixz picnixz Aug 18, 2025

Choose a reason for hiding this comment

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

Took me a while to remember that this was 1 + 0j (maybe rename that variable in a separate issue, e.g., ONE_AS_COMPLEX).

Copy link
Copy Markdown
Member

@skirpichev skirpichev Aug 19, 2025

Choose a reason for hiding this comment

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

This now used just in one place. If we need this branch - it's better to inline struct value here.

Copy link
Copy Markdown
Member

@picnixz picnixz Aug 19, 2025

Choose a reason for hiding this comment

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

If this is the case, then yes, let's inline it.

}

double
Expand Down Expand Up @@ -569,7 +636,9 @@ complex_pow(PyObject *v, PyObject *w, PyObject *z)
// 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);
if (isfinite(a.real) && isfinite(a.imag)) {
_Py_ADJUST_ERANGE2(p.real, p.imag);
}
Comment thread
skirpichev marked this conversation as resolved.
}
else {
p = _Py_c_pow(a, b);
Expand Down