From 2a6a7cb29160af77bfe69ac209ed5df3b7b229c0 Mon Sep 17 00:00:00 2001 From: Aniket Singh Yadav Date: Tue, 1 Sep 2026 09:04:18 +0000 Subject: [PATCH 1/5] gh-156695: Improve accuracy of complex powers with small negative integer exponents --- Lib/test/test_complex.py | 16 +++++++++++++++- Objects/complexobject.c | 4 ++-- 2 files changed, 17 insertions(+), 3 deletions(-) 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/Objects/complexobject.c b/Objects/complexobject.c index 3612c2699a557db..1deccb387acb000 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -357,10 +357,10 @@ c_powu(Py_complex x, long n) static Py_complex c_powi(Py_complex x, long n) { - if (n > 0) + if (n >= 0) return c_powu(x,n); else - return _Py_c_quot(c_1, c_powu(x,-n)); + return c_powu(_Py_c_quot(c_1, x), -n); } From 7c611a9d1af1b2fc65bafadd3dafbcfcb7542f63 Mon Sep 17 00:00:00 2001 From: "blurb-it[bot]" <43283697+blurb-it[bot]@users.noreply.github.com> Date: Tue, 1 Sep 2026 09:09:16 +0000 Subject: [PATCH 2/5] =?UTF-8?q?=F0=9F=93=9C=F0=9F=A4=96=20Added=20by=20blu?= =?UTF-8?q?rb=5Fit.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../2026-09-01-09-09-13.gh-issue-156695.-Gih-8.rst | 4 ++++ 1 file changed, 4 insertions(+) create mode 100644 Misc/NEWS.d/next/Core_and_Builtins/2026-09-01-09-09-13.gh-issue-156695.-Gih-8.rst 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..0b3e88258fbb44c --- /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. From 48e9342be066d27f0c0e0fc42e23d0f99e9723e9 Mon Sep 17 00:00:00 2001 From: Aniket Singh Yadav Date: Tue, 1 Sep 2026 09:20:22 +0000 Subject: [PATCH 3/5] fix-lint --- .../2026-09-01-09-09-13.gh-issue-156695.-Gih-8.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) 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 index 0b3e88258fbb44c..7d96b365ce114e6 100644 --- 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 @@ -1,4 +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, +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. From c0eeaeecd98f31aee308cd30f9fd4e7d6bb73e07 Mon Sep 17 00:00:00 2001 From: Aniket Singh Yadav Date: Tue, 1 Sep 2026 19:13:12 +0000 Subject: [PATCH 4/5] gh-156695: Improve accuracy of complex powers with small negative integer exponents --- Objects/complexobject.c | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/Objects/complexobject.c b/Objects/complexobject.c index 1deccb387acb000..1ce9a6c82d72d3b 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -357,11 +357,29 @@ c_powu(Py_complex x, long n) static Py_complex c_powi(Py_complex x, long n) { - if (n >= 0) + if (n > 0) return c_powu(x,n); - else - return c_powu(_Py_c_quot(c_1, x), -n); + Py_complex r = _Py_c_quot(c_1, c_powu(x, -n)); + if (errno == EDOM + || (isfinite(r.real) && isfinite(r.imag) + && (r.real != 0.0 || r.imag != 0.0))) + return r; + + /* gh-156695: x**|n| left the exponent range although the result is + representable. Redo it with x scaled to exponent zero; both the + scaling and its undoing are exact. */ + double m = fabs(x.real) > fabs(x.imag) ? fabs(x.real) : fabs(x.imag); + if (m == 0.0 || !isfinite(m)) + return r; + + 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 From cdf6a07aede9c36b72ed8d6b325e0bb442134501 Mon Sep 17 00:00:00 2001 From: Aniket Singh Yadav Date: Wed, 2 Sep 2026 08:52:37 +0000 Subject: [PATCH 5/5] Avoid unnecessary complex power recomputation --- Objects/complexobject.c | 38 ++++++++++++++++++++------------------ 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/Objects/complexobject.c b/Objects/complexobject.c index 1ce9a6c82d72d3b..76d3a23e28651cb 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -361,24 +361,26 @@ c_powi(Py_complex x, long n) return c_powu(x,n); Py_complex r = _Py_c_quot(c_1, c_powu(x, -n)); - if (errno == EDOM - || (isfinite(r.real) && isfinite(r.imag) - && (r.real != 0.0 || r.imag != 0.0))) - return r; - - /* gh-156695: x**|n| left the exponent range although the result is - representable. Redo it with x scaled to exponent zero; both the - scaling and its undoing are exact. */ - double m = fabs(x.real) > fabs(x.imag) ? fabs(x.real) : fabs(x.imag); - if (m == 0.0 || !isfinite(m)) - return r; - - 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)); + + /* 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; }