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