Invalid corner cases (resulting in nan+nanj) in _Py_c_prod()

[skirpichev]

Bug report

Bug description:

Reproducer:

>>> z = 1e300+1j
>>> z*complex('(inf+infj)') # should be (nan+infj)
(nan+nanj)

c.f. C code

#include <stdio.h>
#include <complex.h>
#include <math.h>
int main (void)
{
 complex double z = CMPLX(1e300, 1);
 z = CMPLX(INFINITY, INFINITY)*z;
 printf("%e %e\n", creal(z), cimag(z));
 return 0;
}

$ cc a.c && ./a.out
-nan inf

That quickly leads to problems in optimized algorithm for exponentiation (by squaring), that claims to be "more accurate":

>>> z**5
(nan+nanj)
>>> z**5.0000001 # generic algorithm
Traceback (most recent call last):
 File "<python-input-5>", line 1, in <module>
 z**5.0000001
 ~^^~~~~~~~~~
OverflowError: complex exponentiation
>>> _testcapi._py_c_pow(z, 5) # for integer exponent it also signals overflow
((inf+infj), 34)

Following patch (mostly a literal translation of the _Cmultd() routine from the C11 Annex G.5.2) should solve the issue.

a patch

diff --git a/Objects/complexobject.c b/Objects/complexobject.c
index 59c84f1359..1d9707895c 100644
--- a/Objects/complexobject.c
+++ b/Objects/complexobject.c
@@ -53,11 +53,48 @@ _Py_c_neg(Py_complex a)
 }

 Py_complex
-_Py_c_prod(Py_complex a, Py_complex b)
+_Py_c_prod(Py_complex z, Py_complex w)
{
 Py_complex r;
- r.real = a.real*b.real - a.imag*b.imag;
- r.imag = a.real*b.imag + a.imag*b.real;
+ double a = z.real, b = z.imag, c = w.real, d = w.imag;
+ double ac = a*c, bd = b*d, ad = a*d, bc = b*c;
+
+ r.real = ac - bd;
+ r.imag = ad + bc;
+
+ if (isnan(r.real) && isnan(r.imag)){
+ /* Recover infinities that computed as (nan+nanj) */
+ int recalc = 0;
+ if (isinf(a) || isinf(b)){/* z is infinite */
+ /* "Box" the infinity and change nans in the other factor to 0 */
+ a = copysign(isinf(a) ? 1.0 : 0.0, a);
+ b = copysign(isinf(b) ? 1.0 : 0.0, b);
+ if (isnan(c)) c = copysign(0.0, c);
+ if (isnan(d)) d = copysign(0.0, d);
+ recalc = 1;
+ }
+ if (isinf(c) || isinf(d)){/* w is infinite */
+ /* "Box" the infinity and change nans in the other factor to 0 */
+ c = copysign(isinf(c) ? 1.0 : 0.0, c);
+ d = copysign(isinf(d) ? 1.0 : 0.0, d);
+ if (isnan(a)) a = copysign(0.0, a);
+ if (isnan(b)) b = copysign(0.0, b);
+ recalc = 1;
+ }
+ if (!recalc && (isinf(ac) || isinf(bd) || isinf(ad) || isinf(bc))){
+ /* Recover infinities from overflow by changing nans to 0 */
+ if (isnan(a)) a = copysign(0.0, a);
+ if (isnan(b)) b = copysign(0.0, b);
+ if (isnan(c)) c = copysign(0.0, c);
+ if (isnan(d)) d = copysign(0.0, d);
+ recalc = 1;
+ }
+ if (recalc){
+ r.real = Py_INFINITY*(a*c - b*d);
+ r.imag = Py_INFINITY*(a*d + b*c);
+ }
+ }
+
 return r;
 }

c.f. clang's version:
https://github.com/llvm/llvm-project/blob/4973ad47181710d2a69292018cad7bc6f95a6c1a/libcxx/include/complex#L712-L792

Also, maybe _Py_c_prod() code should set errno on overflows, just as _Py_c_abs(). E.g. with above version we have:

>>> z**5
Traceback (most recent call last):
 File "<python-input-1>", line 1, in <module>
 z**5
 ~^^~
OverflowError: complex exponentiation
>>> z*z*z*z*z
(-inf+infj)

If this issue does make sense, I'll provide a patch.

CPython versions tested on:

CPython main branch

Operating systems tested on:

No response

Linked PRs

• gh-120010: fix invalid (nan+nanj) results in _Py_c_prod() #120287

[skirpichev]

See also #119372 and #117999 (another issue for optimized algorithm for exponentiation).

[serhiy-storchaka]

cc @mdickinson

[skirpichev]

PR is ready for review: #120287

[serhiy-storchaka]

Thank you for your contribution @skirpichev.

I do not think that this change should be backported.