gh-117999: Fix small integer powers of complex numbers (GH-124243)

* 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).

Author: serhiy-storchaka

Lib/test/test_complex.py

@@ -84,6 +84,10 @@ def assertClose(self, x, y, eps=1e-9):
         # check that relative difference < eps
         self.assertTrue(abs(x-y)/abs(y) < eps)
 
+    def assertSameSign(self, x, y):
+        if copysign(1., x) != copysign(1., y):
+            self.fail(f'{x!r} and {y!r} have different signs')
+
     def check_div(self, x, y):
         """Compute complex z=x*y, and check that z/x==y and z/y==x."""
         z = x * y
@@ -446,6 +450,63 @@ 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.
+        values = [complex(x, y)
+                  for x in [5, -5, +0.0, -0.0, INF, -INF, NAN]
+                  for y in [12, -12, +0.0, -0.0, INF, -INF, NAN]]
+        for c in values:
+            with self.subTest(value=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**3, (c*c)*c)
+                if not c:
+                    continue
+                for n in range(1, 9):
+                    with self.subTest(exponent=-n):
+                        self.assertComplexesAreIdentical(c**-n, 1/(c**n))
+
+        # Special cases for complex division.
+        for x in [+2, -2]:
+            for y in [+0.0, -0.0]:
+                c = complex(x, y)
+                with self.subTest(value=c):
+                    self.assertComplexesAreIdentical(c**-1, complex(1/x, -y))
+                c = complex(y, x)
+                with self.subTest(value=c):
+                    self.assertComplexesAreIdentical(c**-1, complex(y, -1/x))
+        for x in [+INF, -INF]:
+            for y in [+1, -1]:
+                c = complex(x, y)
+                with self.subTest(value=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(value=c):
+                    self.assertComplexesAreIdentical(c**-1, complex(+0.0*y, -1/x))
+                    self.assertComplexesAreIdentical(c**-2, complex(-0.0, -y/x))
+
+        # Test that zeroes has the same sign as small non-zero values.
+        eps = 1e-11
+        pairs = [(complex(x, y), complex(x, copysign(0.0, y)))
+                 for x in [+1, -1] for y in [+eps, -eps]]
+        pairs += [(complex(y, x), complex(copysign(0.0, y), x))
+                  for x in [+1, -1] for y in [+eps, -eps]]
+        for c1, c2 in pairs:
+            for n in exponents:
+                with self.subTest(value=c1, exponent=n):
+                    r1 = c1**n
+                    r2 = c2**n
+                    self.assertClose(r1, r2)
+                    self.assertSameSign(r1.real, r2.real)
+                    self.assertSameSign(r1.imag, r2.imag)
+                    self.assertNotEqual(r1.real, 0.0)
+                    if n != 0:
+                        self.assertNotEqual(r1.imag, 0.0)
+                    self.assertTrue(r2.real == 0.0 or r2.imag == 0.0)
+
     def test_boolcontext(self):
         for i in range(100):
             self.assertTrue(complex(random() + 1e-6, random() + 1e-6))

Misc/NEWS.d/next/Core_and_Builtins/2024-09-19-15-47-50.gh-issue-117999.Iq4jEG.rst

@@ -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.

Objects/complexobject.c

@@ -26,8 +26,6 @@ class complex "PyComplexObject *" "&PyComplex_Type"
 
 /* elementary operations on complex numbers */
 
-static Py_complex c_1 = {1., 0.};
-
 Py_complex
 _Py_c_sum(Py_complex a, Py_complex b)
 {
@@ -333,23 +331,33 @@ _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;
-    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(0 < n && n <= INT_EXP_CUTOFF);
+    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;
     }
     return r;
 }
@@ -358,10 +366,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_rc_quot(1.0, c_powu(x, -n));
     else
-        return _Py_c_quot(c_1, c_powu(x,-n));
-
+        return (Py_complex){1., 0.};
 }
 
 double
@@ -751,9 +760,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);