diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py index fd002fb00ac338..3facd38e0f732f 100644 --- a/Lib/test/test_complex.py +++ b/Lib/test/test_complex.py @@ -87,6 +87,10 @@ def assertClose(self, x, y, eps=1e-9): self.assertCloseAbs(x.real, y.real, eps) self.assertCloseAbs(x.imag, y.imag, 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 @@ -445,6 +449,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-8 + 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)) diff --git a/Misc/NEWS.d/next/Core_and_Builtins/2024-09-19-15-47-50.gh-issue-117999.Iq4jEG.rst b/Misc/NEWS.d/next/Core_and_Builtins/2024-09-19-15-47-50.gh-issue-117999.Iq4jEG.rst new file mode 100644 index 00000000000000..7c3a9b38372cad --- /dev/null +++ b/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. diff --git a/Objects/complexobject.c b/Objects/complexobject.c index 5d9b3c9f0e3e76..5a82bf3778120e 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -333,7 +333,11 @@ _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; } @@ -341,15 +345,19 @@ _Py_c_pow(Py_complex a, Py_complex b) 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); + 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 c_1; } double @@ -737,7 +746,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); + } } else { p = _Py_c_pow(a, b);