gh-156695: Improve accuracy of complex div and pow

[hpkfft]

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:

perftest-div-pow.py

#!/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	main	patch
cdiv(z, w)	13.1 ns	11.1 ns: 1.19x faster
cpow(z, -5)	26.3 ns	27.0 ns: 1.03x slower
cpow(z, -4)	25.9 ns	26.6 ns: 1.03x slower
cpow(z, -3)	25.5 ns	25.1 ns: 1.01x faster
cpow(z, -1)	26.5 ns	24.7 ns: 1.07x faster
cpow(z, 1)	22.5 ns	21.3 ns: 1.06x faster
cpow(z, 2)	22.8 ns	23.0 ns: 1.01x slower
cpow(z, 3)	23.3 ns	23.5 ns: 1.01x slower
cpow(z, 4)	24.5 ns	24.8 ns: 1.01x slower
cpow(z, 5)	24.4 ns	25.5 ns: 1.04x slower
Geometric mean	(ref)	1.02x faster

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

[skirpichev]

This has merge conflict.

[skirpichev]

Few comments.

This catch-all pr seems unmanageable for me so far.

[skirpichev]

This just makes diff larger.

It's apparent for a reader, that this function is public. Then reader can take look on the header file and find this comment:

/* Operations on complex numbers (soft deprecated
   since Python 3.15). */
PyAPI_FUNC(Py_complex) _Py_c_sum(Py_complex, Py_complex);
PyAPI_FUNC(Py_complex) _Py_c_diff(Py_complex, Py_complex);
PyAPI_FUNC(Py_complex) _Py_c_neg(Py_complex);
PyAPI_FUNC(Py_complex) _Py_c_prod(Py_complex, Py_complex);
PyAPI_FUNC(Py_complex) _Py_c_quot(Py_complex, Py_complex);
PyAPI_FUNC(Py_complex) _Py_c_pow(Py_complex, Py_complex);
PyAPI_FUNC(double) _Py_c_abs(Py_complex);

IMO, this should be reverted.

[hpkfft]

I disagree it's apparent for a reader. Without the comment _Py_c_sum() and _Py_cr_sum() look the same. But the former is public, and the latter is private. It's confusing to me, at least.

I want to avoid setting errno. Look at complex_div(), which now calls c_quot(a, b) to avoid calling _Py_c_quot(a, b) because that sets errno. But complex_div() does call _Py_rc_quot(), which is private, so it doesn't set errno.

[skirpichev]

I disagree it's apparent for a reader.

I don't see any excuse for a reader, who can't take look on header files, where those functions are declared.

I want to avoid setting errno.

I don't think this is a bad, per se. #156145 is not for this. That issue suggests not rely on errno set by libm's functions. Like #156887 does, for example.

You suggest to invent essentially a new private API layer. What are arguments for?

BTW, if we will do this - there be a plenty questions. Pass in/out arguments by reference? Return error status? Rely on floating-point exceptions instead? One possible new place for using low-level private API is bytecode specialization (as we do for floats/ints). I would expect much larger speedup from this.

[hpkfft]

I don't think this is a bad, per se.

On d.p.o, Mark called it evil, and I agree with him. I respect your right to have a different opinion.

[skirpichev]

On d.p.o, Mark called it evil, and I agree with him.

I suspect someone misinterpret his words. Relevant quote is:

float_pow in particular seems unnecessarily evil in that it’s not just reading that global state, but also writing to it (code ref) to force overflow exceptions, so even on boxes whose libm doesn’t touch errno (e.g., macOS), correct operation still relies on writes and reads of errno.

Detailed explanation of the issue in the d.p.o thread is here.

I don't think this fits errno usage pattern e.g. for most of cmath's functions, like atanh(). Here we explicitly set errno, regardless on behavior of libm's functions.

Evil is some particular way of working with errno...

[skirpichev]

That, probably, can go to a separate pr. I can do this, if core devs are interested.

This pr touches many things and it's good to make it in some more manageable pieces.

[hpkfft]

Yes, I would like it if you put this in a separate PR. Then I'll merge with main, and my diff will be smaller!
Also, it's your idea and your work. (I thought it had been merged already, and I was using it locally, so I didn't want to make the effort to revert it locally before pushing this PR.)

Honestly, I've lost track of all the issues and PRs that affect complexobject.c....

[skirpichev]

That's it. I would prefer, if you reduce pr just to deal with underflows in the division. I.e. to fix referenced issue.

BTW, in CPython we assume IEEE doubles (at least as binary format), I think you can use this algorithm unconditionally.

CC @Aniketsy, perhaps you can try this approach in your pr.

[skirpichev]

You can use instead a union type, as in article.

BTW, what are reason for modification of the cdiv() algorithm? I would prefer original version, to simplify review. Unless there is a huge performance boost.

[hpkfft]

The motivation for this cdiv() algorithm is better accuracy. It does not make sense to me to leave a large worst-case error for cdiv() while changing cpow() to be more accurate. Complex division is more fundamental and more widely used than complex power. By improving the accuracy of cdiv(), one gets improved accuracy for cpow(z, -1) without any extra work.
The performance improvement (if I may make a guess) is from eliminating errno. Though, the new cdiv() algorithm has two floating-point divisions instead of three, so that may help performance a little.

[skirpichev]

The motivation for this cdiv() algorithm is better accuracy.

I understand.

My question was: why you aren't using original algorithm, presented in quoted article?

The performance improvement (if I may make a guess) is from eliminating errno.

Another good reason to leave errno stuff for a different pr.