Possible undefined behavior division by zero in complex's `_Py_c_pow()`

[gpshead]

Bug report

Bug description:

https://github.com/python/cpython/blob/main/Objects/complexobject.c#L150 _Py_c_pow() contains:

        at = atan2(a.imag, a.real);
        phase = at*b.real;
        if (b.imag != 0.0) {
            len /= exp(at*b.imag);

An oss-fuzz pycompiler fuzzer identified a problem compiling the code " 9J**33J**3" within the optimizer folding the constant by doing the math in place. We haven't been able to reproduce this ourselves - but code inspection reveals a potential problem:

C exp(x) can return 0 in a couple of cases. Which would lead to the /= executing an undefined division by zero operation.

I believe the correct _Py_c_pow() code should look something like:

        if (b.imag != 0.0) {
            double tmp_exp = exp(at*b.imag);
            if (exp == 0.0 || errno == ERANGE) {
                // Detected as OverflowError by our caller.
                r.real = Py_HUGE_VAL;
                r.imag = Py_HUGE_VAL;
                return r;
            }
            len /= tmp_exp;

CPython versions tested on:

CPython main branch

Operating systems tested on:

Linux

Linked PRs

• gh-113841: fix possible undefined division by 0 in _Py_c_pow() #127211
• [3.13] gh-113841: fix possible undefined division by 0 in _Py_c_pow() (GH-127211) #127216
• [3.12] gh-113841: fix possible undefined division by 0 in _Py_c_pow() (GH-127211) (GH-127216) #127530

[serhiy-storchaka]

Why not use multiplication instead of division?

[gpshead]

The complex object code has been this way since it was added 28 years ago. f9fca92 😅 and may have come from @khinsen.

For floating point things, it is often worth figuring out how error accumulates along the way in the computation to minimize loss. I haven't given the algorithm here any thought or looked up alternatives. My proposal is a mere edit that'd keep it the same with an added check in the middle.

[gpshead]

Also of note since complexobject.c was written... C99 added a complex double type and among the functions available for use with that is cpow()...

[gpshead]

which may be less trivial to use because visual studio support appears to be partial? https://developercommunity.visualstudio.com/t/support-for-c99-complex-numbers/1409049
->
https://learn.microsoft.com/en-us/cpp/c-runtime-library/complex-math-support?view=msvc-170

[khinsen]

Yes, that was me who contributed the complex object to CPython. All the algorithms were taken from C code snippets floating around in my environment - version control hadn't arrived in computational science yet! Since then, not only the C language has evolved, but also the optimization habits of C compilers. I wouldn't waste much effort figuring out why pow was implemented this way. It's more productive to consider how best to implement pow in today's C ecosystem.

[serhiy-storchaka]

I think that the result should be NaN if exp == 0.0 && len == 0.0.

Also, should not errno == ERANGE give zero or NaN, depending on len, as result?

[serhiy-storchaka]

cc @mdickinson

[skirpichev]

We haven't been able to reproduce this ourselves

For me, it looks you example indeed trigger exp=0 condition for code

x, y = 9j, 33j
x**y**3  # here at*b.imag is -0x1.b9036a4a06cfcp+15

The C standard says for division, that "if the value of the second operand is zero, the behavior is undefined". On another hand, with Annex F - you got ERANGE from exp() and +inf (in given example). So, you will correctly get OverflowError.

The *= suggestion will work too. I doubt it's much better on non-IEEE systems. But here is a simple patch: #127211

complex double type and among the functions available for use with that is cpow()

Annex G is still optional. And e.g. M$ CRT implements complex math differently, though _Dcomplex should be compatible with Py_complex, as well as double complex type from the standard ("Each complex type has the same representation and alignment requirements as an array type containing exactly two elements of the corresponding real type"). So, in principle this could be used for _Py_c_pow() and for cmath's functions.

I'm not sure about quality of implementations. But e.g. the glibc versions of complex math functions - pass the CPython's cmath_testcases.txt.

It's more productive to consider how best to implement pow in today's C ecosystem.

It seems, glibc just do exp(e*log(b)): https://sourceware.org/git/?p=glibc.git;a=blob;f=math/s_cpow_template.c;h=92e3e986639be118ec36fb7f991ba47d8115816e;hb=HEAD

Same I see for generic case e.g. in the numpy.