gh-156695: Improve accuracy of complex powers with small negative integer exponents - #156757

# 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)

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}%")