GeezerVison

Background

I started this work after finding bug #153144 when building Python with Intel's math library. @skirpichev made me aware that the same bug affects Solaris.
PR #153148 is awaiting core review. In a comment in that PR, I observed:

Some math libraries (e.g., musl) don't set errno for anything, so Python cannot rely on errno for detecting overflow or invalid. I would think that errno checking can be removed everywhere....

Less that a day later, by coincidence, Mark observed on d.p.o that it's possible to remove Python's dependence on errno:

It would be nice to free the code from reliance on "global" (really thread-local) state for what should be a pure function, for all sorts of reasons - easier to reason about, opportunities for optimization, etc.

That discussion began with the OP (unintentionally) compiling Python with the flag -fno-math-errno, giving the compiler false information that leads it to perform incorrect optimizations. Thanks, again, to Sergey for bringing that to my attention and for creating feature request #156145.
Towards this goal, my PR #156551 for floatobject.c is awaiting core review.

This PR proposes the next step, i.e., updating complexobject.c.

It mostly makes #155527 obsolete, but that PR also includes tests that either should be added here or merged separately.

This PR includes Sergey's inline complex number utility classifying functions from PR #156694, which I thought had been merged, but actually was closed.

This PR should fix #156886, but tests are needed.

Anticipating pushback, there's another point I'd like to make. The documentation on C API and ABI Stability says

Names prefixed by an underscore, such as _Py_InternalState, are private API that can change without notice even in patch releases.

Also, Stable ABI Caveats warns:

Note that compiling for Stable ABI is not a complete guarantee that code will be compatible with the expected Python versions. Stable ABI prevents ABI issues, like linker errors due to missing symbols or data corruption due to changes in structure layouts or function signatures. However, other changes in Python can change the behavior of extensions.

This PR changes the behavior of _Py_c_quot and _Py_c_pow. I assume these functions are not somehow extra-stable because they have been soft deprecated. In my opinion, the C API is a C API to Python functionality, and so this PR
adds the small integer exponent optimization from Python's pow() function into the C API _Py_c_pow.

Doing so led me to question the source code comment

  // Check whether the exponent has a small integer value, and if so use
  // a faster and more accurate algorithm.

because I did not want users of _Py_c_pow (e.g., Numba) to encounter accuracy issues if this optimization were enabled for them. On the other hand, if the optimization is valuable, they should be allowed to benefit from it.

I shared my early findings with Sergey, who created issue #156695.

I'm using GCC 14.2.0 and GNU C Library 2.41.
Hardware is Intel Xeon (Sapphire Rapids).

Accuracy

I'm using the methodology described in https://inria.hal.science/hal-04714173 (though, here I've done a single, faster run, which is not as thorough as we did for the paper, so some variation is to be expected. A pseudo-random number generator is used to make test points, and it is intentionally seeded differently for each run.)

For Python's complex number division, I found examples having normwise error of over 9.0e15 ulps:

(0x1.d566ef31c4f9fp+36,-0x1.f63b2c46dc0f8p+114) / (0x1.bcaf39465046cp+1023,-0x1.d2dea988c0ef9p+1022)
(0x1.e65cd27ec1127p+121,-0x1.479364664ad5ap+40) / (0x1.c2f336389dd36p+1022,-0x1.aef0ba5f1beb1p+1023)

Note the large exponents in the denominators above.

With this PR, the worst error I found for division is less than 4.48 ulps for

(0x1.598bf798aeceep-466,-0x1.69f042a45c73fp-457) / (-0x1.6a7df2c08ffb9p-540,0x1.15d11bd015608p-554)

Note this this is not a proof of worst case error. It's just what I've been able to find.

For Python's complex power to a small negative integer exponent, I found examples having error over 1.5e15 ulps:

(-0x1.000001p+1023,-0x1.001p+1023) ** -1
(-0x1.49852f983ef92p+510,0x1.8dc42193d5c5ep+511) ** -2

and so on. Again, note the large magnitude of z.

The graph below shows the error for z**n for fixed n, with the "worst" z found by the test program:

• Red -- Python's current _Py_c_pow (which does not have the small integer exponent optimization) for integer 0 <= n <= 100.
• Black -- Python's current pow() error for integer 0 <= n <= 100.
• Green -- This PR's error for integer -100 <= n <= 100.

Performance

Using pyperf:

#!/usr/bin/env ./python
import time
import pyperf

def bench_cdiv(loops):
    z = complex(3.140625, 1.0)
    w = complex(2.718750, 4.0)
    range_it = range(loops)
    t0 = time.perf_counter()
    for _ in range_it:
        result = z/w; result = z/w; result = z/w; result = z/w;
        result = z/w; result = z/w; result = z/w; result = z/w;
        result = z/w; result = z/w; result = z/w; result = z/w;
        result = z/w; result = z/w; result = z/w; result = z/w;
        result = z/w; result = z/w; result = z/w; result = z/w;
    return time.perf_counter() - t0

def bench_cpow(loops, n):
    z = complex(3.140625, 1.0)
    range_it = range(loops)
    t0 = time.perf_counter()
    for _ in range_it:
        result = z**n; result = z**n; result = z**n; result = z**n;
        result = z**n; result = z**n; result = z**n; result = z**n;
        result = z**n; result = z**n; result = z**n; result = z**n;
        result = z**n; result = z**n; result = z**n; result = z**n;
        result = z**n; result = z**n; result = z**n; result = z**n;
    return time.perf_counter() - t0

runner = pyperf.Runner(values=10, processes=24)

runner.bench_time_func(f"cdiv(z, w)", bench_cdiv, inner_loops=20)

for n in (-5, -4, -3, -2, -1, 1, 2, 3, 4, 5):
    runner.bench_time_func(f"cpow(z, {n})", bench_cpow, n, inner_loops=20)

Benchmark hidden because not significant (1): cpow(z, -2)

Next steps

The internal function c_square doesn't seem helpful. Probably, it should just be deleted.

Consider special values (inf, nan) and add tests.
Also, some current tests are designed to verify the current implementation and need to be updated accordingly.

If a core developer is interested in working with me with the intention of then merging this, please let me know.
It's been fun, but I don't want to invest more time otherwise.

CC: @serhiy-storchaka, @eendebakpt

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