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

[Aniketsy]

Fixes #156695

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

[eendebakpt]

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.

[Aniketsy]

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]

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]

Any example, that trigger that case?

[skirpichev]

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