diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py index bb307191dffcc14..e0129faaccee0dd 100644 --- a/Lib/test/test_complex.py +++ b/Lib/test/test_complex.py @@ -8,7 +8,7 @@ ) from random import random -from math import isnan, copysign +from math import isnan, copysign, ulp import operator INF = float("inf") @@ -446,6 +446,20 @@ def test_pow_with_small_integer_exponents(self): self.assertEqual(str(float_pow), str(int_pow)) self.assertEqual(str(complex_pow), str(int_pow)) + @support.requires_IEEE_754 + def test_pow_small_negative_integer_exponents(self): + z = complex(float.fromhex('0x1.47e9c711723f5p+81'), + float.fromhex('0x1.38afd1168e49fp+85')) + expected = complex(float.fromhex('0x0.4000000000000p-1022'), + float.fromhex('0x0.3ffffffffffffp-1022')) + for exponent in (-12, -12.0, complex(-12.0, 0.0)): + with self.subTest(exponent=exponent): + result = z ** exponent + self.assertLessEqual(abs(result.real - expected.real), + 4 * ulp(expected.real)) + self.assertLessEqual(abs(result.imag - expected.imag), + 4 * ulp(expected.imag)) + 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/2026-09-01-09-09-13.gh-issue-156695.-Gih-8.rst b/Misc/NEWS.d/next/Core_and_Builtins/2026-09-01-09-09-13.gh-issue-156695.-Gih-8.rst new file mode 100644 index 000000000000000..7d96b365ce114e6 --- /dev/null +++ b/Misc/NEWS.d/next/Core_and_Builtins/2026-09-01-09-09-13.gh-issue-156695.-Gih-8.rst @@ -0,0 +1,4 @@ +Improve accuracy of :class:`complex` powers with small negative integer +exponents. Previously ``z**-n`` was computed as ``1/(z**n)``; the +intermediate ``z**n`` can overflow even when the result is representable, +in which case all precision was lost. diff --git a/Objects/complexobject.c b/Objects/complexobject.c index 3612c2699a557db..76d3a23e28651cb 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -359,9 +359,29 @@ c_powi(Py_complex x, long n) { if (n > 0) return c_powu(x,n); - else - return _Py_c_quot(c_1, c_powu(x,-n)); + Py_complex r = _Py_c_quot(c_1, c_powu(x, -n)); + + /* gh-156695: x**|n| needs roughly twice the exponent range of the + result, so it can leave the range even when the result itself is + representable, leaving the quotient degenerate. Only then redo the + computation with x scaled to exponent zero; both the scaling and its + undoing are exact. The common path above is untouched. */ + if (!(isfinite(r.real) && isfinite(r.imag) + && (r.real != 0.0 || r.imag != 0.0)) + && errno != EDOM) + { + double m = fabs(x.real) > fabs(x.imag) ? fabs(x.real) : fabs(x.imag); + if (m != 0.0 && isfinite(m)) { + int e; + frexp(m, &e); + Py_complex w = {ldexp(x.real, -e), ldexp(x.imag, -e)}; + r = _Py_c_quot(c_1, c_powu(w, -n)); + r.real = ldexp(r.real, (int)(e * n)); + r.imag = ldexp(r.imag, (int)(e * n)); + } + } + return r; } double