gh-156693: correct complex powers with infinite bases

[skirpichev]

Now we raise OverflowError's only if infinite result (i.e. cmath.isinf(result) is True) come from finite input.

This patch also adds private inline functions _Py_c_iszero(), _Py_c_isnan(), _Py_c_isinf() and _Py_c_isfinite().

• Issue: Incorrect overflows for complex powers when base is infinite #156693

[skirpichev]

CC @serhiy-storchaka

I think this will be less controversial change vs #118000 or #124243. Anyway, I think it worth for a dedicated commit.

[skirpichev]

CC @hpkfft

[bedevere-app]

A Python core developer has requested some changes be made to your pull request before we can consider merging it. If you could please address their requests along with any other requests in other reviews from core developers that would be appreciated.

Once you have made the requested changes, please leave a comment on this pull request containing the phrase I have made the requested changes; please review again. I will then notify any core developers who have left a review that you're ready for them to take another look at this pull request.

[skirpichev]

I have made the requested changes; please review again

[bedevere-app]

Thanks for making the requested changes!

@eendebakpt: please review the changes made to this pull request.

[eendebakpt]

One more caseswhere the behavior changes: (1e-200+0j)**-2 (it would not happen under #156695, but that is not merged yet).

And can you add a test for complex(inf)**1j?

[skirpichev]

One more caseswhere the behavior changes: (1e-200+0j)**-2

Indeed, I believe now all should be fixed.

I'll not add a test, though, because this behavior can be changed by #156695.

And can you add a test for complex(inf)**1j?

Done.

Last commit removes special handling for some corner cases, that's optional. Let me know if you prefer to keep this.

[eendebakpt]

With your latest iteration: what about (0j)**complex(INF), (1e-200+0j)**-2 and _Py_c_pow(nan+0j, 0+0j) (C API)?

[hpkfft]

In my opinion, anybody calling _Py_c_pow() (e.g., Numba) should get the same results as a Python developer calling pow(). The latter has an optimization for small integer exponents. Numba should be allowed to benefit from that.

You're welcome to look at my work-in-progress: hpkfft#1
I deleted c_powu(). Now, complex_pow() only calls c_pow() and raises a Python error if necessary. The C API function _Py_c_pow() also calls c_pow() and sets errno if necessary. If small integer exponent, c_pow() calls c_powi(), so the C API and the Python API give the same results.

[skirpichev]

In my opinion, anybody calling _Py_c_pow() (e.g., Numba) should get the same results as a Python developer calling pow().

Anyone, calling _Py_c_pow() should be aware, that it's a soft-deprecated interface.

We decided to not provide low-level primitives, that works same as complex arithmetics in the Python. You can use either PyNumber_* API or some external library code (e.g. like GNU GSL) and/or C types (e.g. double _Complex, if available).

You're welcome to look at my work-in-progress: hpkfft#1

Thanks, but I'll try to limit this PR just for a bugfix.

IMO, cleanest alternative to using errno is floating-point exceptions. Though, my tests shows that this come with a severe speed loss. Lets decide on this later.

[skirpichev]

Ah, it seems that added complex('inf')**1j test is still in mercy of platforms libm. Test pass for me, but breaks CI. I'll revert it.

BTW, I'm not sure that ZeroDivisionError is a right thing here. On latest gmpy2:

>>> from gmpy2 import *
>>> get_context().trap_divzero = True
>>> mpc('inf')**mpc(1j)
mpc('nan+nanj')
>>> get_context()
context(precision=53, real_prec=Default, imag_prec=Default,
        round=RoundToNearest, real_round=Default, imag_round=Default,
        emax=1073741823, emin=-1073741823,
        subnormalize=False,
        trap_underflow=False, underflow=False,
        trap_overflow=False, overflow=False,
        trap_inexact=False, inexact=False,
        trap_invalid=False, invalid=True,
        trap_erange=False, erange=False,
        trap_divzero=True, divzero=False,
        allow_complex=False,
        rational_division=False,
        allow_release_gil=False)

[hpkfft]

I'm not sure that ZeroDivisionError is a right thing here.

My preference would be ZeroDivisionError only if base is zero and exp is not a positive real number.
The following seems obviously wrong:

>>> import math
>>> c = complex(math.inf)
>>> c**1j
Traceback (most recent call last):
  File "<python-input-7>", line 1, in <module>
    c**1j
    ~^^~~
ZeroDivisionError: 0.0 to a negative or complex power

[serhiy-storchaka]

Opened #156886 for ZeroDivisionError.

[eendebakpt]

@skirpichev Do you plan to followup on this PR? Some of the added tests seem worthwhile to have in any case.