gh-156695: Improve accuracy of complex powers with small negative integer exponents by Aniketsy · Pull Request #156757 · python/cpython

Aniketsy · 2026-09-01T09:06:28Z

>>> import math
>>> z = complex(float.fromhex('0x1.47e9c711723f5p+81'),
...            float.fromhex('0x1.38afd1168e49fp+85'))
>>> ref = complex(float.fromhex('0x0.4000000000000p-1022'),
...            float.fromhex('0x0.3ffffffffffffp-1022'))
>>> z ** -12
0j
>>> abs((z**-12 - ref).real) / math.ulp(ref.real)
1125899906842624.0
>>> 0.0j ** 0
(1+0j)

Issue: Improve accuracy for complex powers with small negative integer exponents #156695

…ve integer exponents

eendebakpt · 2026-09-01T12:12:42Z

The result is improving for the example of the OP, but some results are worse as well. E.g. (3+4j)**-100 or (-2.24e-4+5.09e-5j)**-7 or (8087.7392089611985 + 8087.4504765395295j)**-2. I think we need a more extensive analysis of the overall impact on the results.

skirpichev

Yes, unfortunately this is not an easy issue.

You could compare old/new results with correctly rounded powers (using e.g. GNU MPC) to see if the net impact is positive. Take look on https://inria.hal.science/hal-04714173 for inspiration.

…ve integer exponents

Aniketsy · 2026-09-01T19:16:00Z

The result is improving for the example of the OP, but some results are worse as well. E.g. (3+4j)-100 or (-2.24e-4+5.09e-5j)-7 or (8087.7392089611985 + 8087.4504765395295j)**-2. I think we need a more extensive analysis of the overall impact on the results.

thanks for pointing out i dived into this, these are some results with these scripts i used

Details

# compare.py -- normwise error vs correctly rounded MPC results
import math, random, statistics
from gmpy2 import mpc, get_context
get_context().precision = 200

def normwise(got, ref):
    # Eq. (1) of Caprioli, Innocente & Zimmermann: e = |delta| / ulp(|z|)
    if got is None: return float('inf')
    if math.isnan(got.real) or math.isnan(got.imag): return None
    if math.isinf(got.real) or math.isinf(got.imag): return None
    return float(abs(mpc(got) - ref)) / math.ulp(float(abs(ref)))

def measure(pow_fn, samples):
    errs = []
    for z, n in samples:
        ref = mpc(z) ** n
        try: got = pow_fn(z, n)
        except (OverflowError, ZeroDivisionError): got = None
        errs.append(normwise(got, ref))
    excluded = sum(1 for e in errs if e is None)
    scored = [e for e in errs if e is not None]
    return (statistics.median(scored),
            statistics.quantiles(scored, n=100)[94],
            max(scored),
            excluded,
            errs)

random.seed(0)
samples = []
for _ in range(50000):
    n = -random.choice((1, 2, 3, 5, 7, 12, 40, 100))
    z = complex(random.uniform(-1, 1), random.uniform(-1, 1))
    z *= math.ldexp(1.0, random.randint(-1020, 1020))
    if z == 0 or not math.isfinite(abs(z)):
        continue
    ref = mpc(z) ** n
    if abs(ref) == 0 or not math.isfinite(float(abs(ref))):
        continue
    samples.append((z, n))

med, p95, mx, bad, errs = measure(lambda z, n: z ** n, samples)
print(f"n={len(samples)}  median={med:.3f}  p95={p95:.3f}  max={mx:.4g}  excluded={bad}")

import json, sys
if len(sys.argv) > 1:
    with open(sys.argv[1], "w") as f:
        json.dump([e if (e is not None and math.isfinite(e)) else None for e in errs], f)

                n        median    p95    max         excluded (NaN/Inf)

  main          32094    0.492   1.993  9.537e+14   17687
  patched       32094    0.000   1.250  63.22           0

Details

import json, statistics

a = json.load(open('errs_main.json'))
b = json.load(open('errs_patched.json'))
print(f"lengths: main={len(a)}  patched={len(b)}")

both = [(x, y) for x, y in zip(a, b) if x is not None and y is not None]
A = [x for x, _ in both]
B = [y for _, y in both]
q = lambda L: statistics.quantiles(L, n=100)[94]

print(f"common pool: {len(both)}")
print(f"  main    median={statistics.median(A):.3f}  p95={q(A):.3f}  max={max(A):.4g}")
print(f"  patched median={statistics.median(B):.3f}  p95={q(B):.3f}  max={max(B):.4g}")
print(f"  identical: {100 * sum(1 for x, y in both if x == y) / len(both):.2f}%")

common pool: 14407 samples (both builds return a finite, non-NaN result)
  main     median=0.492  p95=1.993  max=9.537e+14
  patched  median=0.492  p95=1.986  max=63.22
  identical values: 99.97%

Yes, unfortunately this is not an easy issue.

You could compare old/new results with correctly rounded powers (using e.g. GNU MPC) to see if the net impact is positive. Take look on https://inria.hal.science/hal-04714173 for inspiration.

yes it got trickier than i thought, and thanks for the reference paper

skirpichev · 2026-09-02T01:51:36Z

With original patch I've this:

ref:
n=29187  median=0.500  p95=2.358  max=9.537e+14  excluded=0
patch:
n=29187  median=1.000  p95=4.854  max=120.4  excluded=0

Details

# compare.py -- normwise error vs correctly rounded MPC results
import cmath, math, random, statistics
from gmpy2 import mpc, ieee, set_context

set_context(ieee(64))

def normwise(got, ref):
    # Eq. (1) of Caprioli, Innocente & Zimmermann: e = |delta| / ulp(|z|)
    if got is None:
        return float('inf')
    if cmath.isnan(got):
        return
    if cmath.isinf(got):
        return
    diff = abs(got - ref)
    if not math.isfinite(diff):
        return
    return diff/math.ulp(abs(ref))

def measure(pow_fn, samples):
    errs = []
    for z, n in samples:
        ref = complex(mpc(z) ** n)
        try:
            got = pow_fn(z, n)
        except (OverflowError, ZeroDivisionError):
            got = None
        errs.append(normwise(got, ref))
    excluded = sum(1 for e in errs if e is None)
    scored = [e for e in errs if e is not None]
    return (statistics.median(scored),
            statistics.quantiles(scored, n=100)[94],
            max(scored),
            excluded,
            errs)

random.seed(0)
samples = []
for _ in range(100000):
    n = -random.choice((1, 2, 3, 5, 7, 12, 40, 100))
    z = complex(random.uniform(-1, 1), random.uniform(-1, 1))
    z *= math.ldexp(1.0, random.randint(-1020, 1020))
    if z == 0 or not cmath.isfinite(z):
        continue
    ref = complex(mpc(z) ** n)
    if ref == 0 or not math.isfinite(abs(ref)):
        continue
    samples.append((z, n))

med, p95, mx, bad, errs = measure(lambda z, n: z ** n, samples)
print(f"n={len(samples)}  median={med:.3f}  p95={p95:.3f}  max={mx:.4g}  excluded={bad}")

import json, sys

if len(sys.argv) > 1:
    with open(sys.argv[1], "w") as f:
        json.dump([e if (e is not None and math.isfinite(e)) else None for e in errs], f)

skirpichev · 2026-09-02T01:53:34Z

+    if (errno == EDOM
+        || (isfinite(r.real) && isfinite(r.imag)
+            && (r.real != 0.0 || r.imag != 0.0)))
+        return r;


Any example, that trigger that case?

errno == EDOM covers x**|n| for being exactly zero: 0j ** -1, and underflow cases like (1e-200+1e-200j) ** -5 or (5e-324+0j) ** -2

skirpichev · 2026-09-02T01:55:32Z

+    /* 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. */


You recompute power and quotient again, unconditionally. I believe it will introduce a severe speed regression.

i've done some improvement in this to avoid recomputation

these are results of speed regression.

main patched z**-2 232 237 (+5) z**-100 262 262 ( 0) z**2 223 211 (-12)

also i got this on running script

n=29187 median=0.500 p95=2.062 max=62.94 excluded=0

these are results of speed regression.

Could you, please, share your benchmark?

aniket@DESKTOP-074O80J:/mnt/d/cpython/cpython$ ./python.exe -m timeit -s "z=3+4j" "z**-2" 1000000 loops, best of 5: 232 nsec per loop aniket@DESKTOP-074O80J:/mnt/d/cpython/cpython$ ./python.exe -m timeit -s "z=3+4j" "z**-100" 1000000 loops, best of 5: 262 nsec per loop aniket@DESKTOP-074O80J:/mnt/d/cpython/cpython$ ./python.exe -m timeit -s "z=3+4j" "z**2" 1000000 loops, best of 5: 223 nsec per loop

ahh sorry, after looking into the results i shared was not much clear, this is how i tested on main and on my patch

timeit is a very unreliable tool.

oh i'm not familiar with this, so i need to look into, what's the best way to test, could you please give some reference so that i can go through.

But what about z, that trigger new code?

aniket@DESKTOP-074O80J:/mnt/d/cpython/cpython$ ./python.exe -m timeit -s "z=3+4j" "z**-2" 1000000 loops, best of 5: 237 nsec per loop aniket@DESKTOP-074O80J:/mnt/d/cpython/cpython$ ./python.exe -m timeit -s "z=3+4j" "z**-100" 1000000 loops, best of 5: 262 nsec per loop aniket@DESKTOP-074O80J:/mnt/d/cpython/cpython$ ./python.exe -m timeit -s "z=3+4j" "z**2" 1000000 loops, best of 5: 211 nsec per loop

could you please give some reference so that i can go through.

Take look on benchmarking tools in other prs, for example.

aniket@DESKTOP-074O80J:/mnt/d/cpython/cpython$ ./python.exe -m timeit -s "z=3+4j" "z**-2"

Don't repeat yourself.

I meant z values, that will require to run the new code, new algorithm, not ones that utilize condition, that rule it out.

i re-ran the benchmark with pyperf using inputs (also took help from claude in running this benchmark)

Details

import pyperf, sys Z1 = "complex(float.fromhex('0x1.47e9c711723f5p+81'), float.fromhex('0x1.38afd1168e49fp+85'))" Z2 = "complex(4.221902562142989e+127, 2.821506292740584e+126)" WORKLOADS = [ ("pow2_fast", "z=3+4j", "z**2"), ("powm2_fast", "z=3+4j", "z**-2"), ("powm100_fast","z=3+4j", "z**-100"), ("powm12_slow", f"z={Z1}", "z**-12"), ("powm7_slow", f"z={Z2}", "z**-7"), ] runner = pyperf.Runner() for name, setup, stmt in WORKLOADS: runner.timeit(name=name, stmt=stmt, setup=setup)

and got these results:-

aniket@DESKTOP-074O80J:~/cpython-bench$ ./python.exe -m pyperf compare_to main.json patched.json --verbose pow2_fast ========= Mean +- std dev: [main] 62.7 ns +- 10.7 ns -> [patched] 58.5 ns +- 3.0 ns: 1.07x faster Significant (t=4.16) powm2_fast ========== Mean +- std dev: [main] 81.6 ns +- 5.5 ns -> [patched] 81.7 ns +- 14.5 ns: 1.00x slower Not significant! powm100_fast ============ Mean +- std dev: [main] 105 ns +- 8 ns -> [patched] 107 ns +- 18 ns: 1.02x slower Not significant! powm12_slow =========== Mean +- std dev: [main] 252 ns +- 17 ns -> [patched] 268 ns +- 41 ns: 1.06x slower Significant (t=-3.85) powm7_slow ========== Mean +- std dev: [main] 166 ns +- 19 ns -> [patched] 167 ns +- 23 ns: 1.01x slower Not significant! Geometric mean: 1.00x slower

Sorry, are you sure you run benchmark correctly? Looks like noise.

I run this on your patch with fixed merge conflicts. As expected, it's more than 2x slower on new code:

pow2_fast: Mean +- std dev: [ref] 184 ns +- 8 ns -> [patch] 178 ns +- 1 ns: 1.03x faster powm100_fast: Mean +- std dev: [ref] 320 ns +- 2 ns -> [patch] 293 ns +- 2 ns: 1.09x faster powm12_slow: Mean +- std dev: [ref] 263 ns +- 12 ns -> [patch] 618 ns +- 5 ns: 2.35x slower powm7_slow: Mean +- std dev: [ref] 268 ns +- 9 ns -> [patch] 633 ns +- 14 ns: 2.37x slower Benchmark hidden because not significant (1): powm2_fast Geometric mean: 1.37x slower

What if you will use simpler algorithm from the issue?

Details

diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py index 3d02bb6ec23..2fcfa02a6fd 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") @@ -517,6 +517,20 @@ def test_pow_with_small_integer_exponents(self): self.assertNotEqual(r1.imag, 0.0) self.assertTrue(r2.real == 0.0 or r2.imag == 0.0) + @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 9328baf013c..f931e6d21d9 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -369,11 +369,32 @@ static Py_complex c_powi(Py_complex x, long n) { if (n > 0) - return c_powu(x, n); - else if (n < 0) - return _Py_rc_quot(1.0, c_powu(x, -n)); - else + return c_powu(x,n); + else if (n == 0) return (Py_complex){1., 0.}; + + Py_complex r = _Py_rc_quot(1.0, 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_rc_quot(1.0, c_powu(w, -n)); + r.real = ldexp(r.real, (int)(e * n)); + r.imag = ldexp(r.imag, (int)(e * n)); + } + } + return r; } double

Edit:

I got best results with a following patch:

diff --git a/Objects/complexobject.c b/Objects/complexobject.c index 9328baf013c..096b14bdd3d 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -369,11 +369,26 @@ static Py_complex c_powi(Py_complex x, long n) { if (n > 0) - return c_powu(x, n); - else if (n < 0) - return _Py_rc_quot(1.0, c_powu(x, -n)); - else + return c_powu(x,n); + else if (n == 0) return (Py_complex){1., 0.}; + + double m = fabs(x.real) > fabs(x.imag) ? x.real : x.imag; + + if (m && isfinite(m)) { + int e; + + frexp(m, &e); + + if (-e*n > 800) { + x = (Py_complex){ldexp(x.real, -e), ldexp(x.imag, -e)}; + x = _Py_rc_quot(1.0, c_powu(x, -n)); + x.real = ldexp(x.real, (int)(e * n)); + x.imag = ldexp(x.imag, (int)(e * n)); + return x; + } + } + return _Py_rc_quot(1.0, c_powu(x, -n)); } double

pow2_fast: Mean +- std dev: [ref] 184 ns +- 8 ns -> [patch3] 176 ns +- 1 ns: 1.05x faster powm2_fast: Mean +- std dev: [ref] 242 ns +- 10 ns -> [patch3] 260 ns +- 2 ns: 1.07x slower powm100_fast: Mean +- std dev: [ref] 320 ns +- 2 ns -> [patch3] 313 ns +- 2 ns: 1.02x faster powm12_slow: Mean +- std dev: [ref] 263 ns +- 12 ns -> [patch3] 552 ns +- 7 ns: 2.10x slower powm7_slow: Mean +- std dev: [ref] 268 ns +- 9 ns -> [patch3] 536 ns +- 23 ns: 2.01x slower Geometric mean: 1.33x slower

Sorry, are you sure you run benchmark correctly? Looks like noise.

i tried to run benchmark correctly, but i had also some doubt related to noise, so you maybe right here, there would some mistakes from my ends, thank you for investing your time in this, and running benchmark as well improved patch.

the patch you shared looks good to me, in first look, ( sorry i may take time in this, with updating, as i'm not well and will be going hospital, i will catch-up with this soon )

thanks again, for the review :)

skirpichev · 2026-09-05T01:59:57Z

+        z = complex(float.fromhex('0x1.47e9c711723f5p+81'),
+                    float.fromhex('0x1.38afd1168e49fp+85'))
+        expected = complex(float.fromhex('0x0.4000000000000p-1022'),
+                           float.fromhex('0x0.3ffffffffffffp-1022'))


I think there no big sense in using hexadecimal input. Use decimal literals.

skirpichev · 2026-09-05T02:03:27Z

+                self.assertLessEqual(abs(result.real - expected.real),
+                                     4 * ulp(expected.real))
+                self.assertLessEqual(abs(result.imag - expected.imag),
+                                     4 * ulp(expected.imag))


Why not use assertAlmostEqual?

Why not use assertAlmostEqual?

make sense to me, i'll update

skirpichev · 2026-09-05T02:03:41Z

                    self.assertEqual(str(float_pow), str(int_pow))
                    self.assertEqual(str(complex_pow), str(int_pow))

+    @support.requires_IEEE_754


its just to support, to run test only on system that uses standard IEEE 754 format for floating-point numbers.

i think i have added this after looking into some test, in case if im missing something i'll verify this.

skirpichev · 2026-09-05T02:05:36Z

+          && (r.real != 0.0 || r.imag != 0.0))
+        && errno != EDOM)
+    {
+        double m = fabs(x.real) > fabs(x.imag) ? fabs(x.real) : fabs(x.imag);


There is not need to compute fabs twice. Just use components: sign doesn't matter for computing of the exponent.

You can also compute e unconditionally and then apply scaling or (1/x)**-n algorithme, based on e value.

skirpichev · 2026-09-05T02:12:35Z

+        if (m != 0.0 && isfinite(m)) {
+            int e;
+            frexp(m, &e);
+            Py_complex w = {ldexp(x.real, -e), ldexp(x.imag, -e)};


You can use r.

Aniketsy · 2026-09-09T05:36:40Z

@skirpichev i've updated with your patch as we got best results, and updated with tests ( as now our input are exact powers of two, so the expected values are exact so used assertEqual rather than assertAlmostEqual ).. please let me know if we need improvement, or you want me to verify anything from my end

skirpichev · 2026-09-09T05:47:49Z

i've updated with your patch as we got best results

This still 2x slower.

I suggested you to try Priest's algorithm as in #156968, instead of the Smith code. The article has C-coded version of it.

Aniketsy · 2026-09-09T05:57:58Z

i've updated with your patch as we got best results

This still 2x slower.

I suggested you to try Priest's algorithm as in #156968, instead of the Smith code. The article has C-coded version of it.

ah you mentioned me there, but sorry i didn't get notified ... i'll go through this, till then marking in draft

pythongh-156695: Improve accuracy of complex powers with small negati…

2a6a7cb

…ve integer exponents

bedevere-app Bot mentioned this pull request Sep 1, 2026

Improve accuracy for complex powers with small negative integer exponents #156695

Open

bedevere-app Bot added the awaiting review label Sep 1, 2026

📜🤖 Added by blurb_it.

7c611a9

skirpichev self-requested a review September 1, 2026 09:16

fix-lint

48e9342

skirpichev reviewed Sep 1, 2026

View reviewed changes

pythongh-156695: Improve accuracy of complex powers with small negati…

c0eeaee

…ve integer exponents

skirpichev reviewed Sep 2, 2026

View reviewed changes

Avoid unnecessary complex power recomputation

cdf6a07

skirpichev self-requested a review September 3, 2026 05:18

skirpichev reviewed Sep 5, 2026

View reviewed changes

Aniketsy added 3 commits September 8, 2026 17:24

resolve conflict with new patch

11c2580

improve tests

8f02c1e

imporve formatting in code

a338b40

Aniketsy marked this pull request as draft September 9, 2026 05:58

bedevere-app Bot removed the awaiting review label Sep 9, 2026

Uh oh!

Uh oh!

Conversation

Aniketsy commented Sep 1, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

eendebakpt commented Sep 1, 2026

Uh oh!

skirpichev left a comment

Choose a reason for hiding this comment

Uh oh!

Aniketsy commented Sep 1, 2026

Uh oh!

skirpichev commented Sep 2, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

skirpichev Sep 5, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Aniketsy commented Sep 9, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

skirpichev commented Sep 9, 2026

Uh oh!

Aniketsy commented Sep 9, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

Aniketsy commented Sep 1, 2026 •

edited

Loading

skirpichev Sep 5, 2026 •

edited

Loading

Aniketsy commented Sep 9, 2026 •

edited

Loading