from 120010: fix invalid (nan+nanj) results in _Py_c_prod()

In some cases, previously computed as (nan+nanj), we could recover
meaningful component values in the result, see e.g. the C11, Annex
G.5.2, routine _Cmultd():

>>> z = 1e300+1j
>>> z*(nan+infj)  # was (nan+nanj)
(-inf+infj)

That also fix some complex powers for small integer exponents, computed
with optimised algorithm (by squaring):

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

Author: skirpichev

Lib/test/test_complex.py

@@ -348,6 +348,43 @@ def test_mul(self):
         self.assertComplexesAreIdentical((1+0j) * float(1.0), complex(1, 0))
         self.assertComplexesAreIdentical(float(1.0) * (1+0j), complex(1, 0))
 
+        self.assertComplexesAreIdentical((1e300+1j) * complex(INF, INF),
+                                         complex(NAN, INF))
+        self.assertComplexesAreIdentical(complex(INF, INF) * (1e300+1j),
+                                         complex(NAN, INF))
+        self.assertComplexesAreIdentical((1e300+1j) * complex(NAN, INF),
+                                         complex(-INF, INF))
+        self.assertComplexesAreIdentical(complex(NAN, INF) * (1e300+1j),
+                                         complex(-INF, INF))
+        self.assertComplexesAreIdentical((1e300+1j) * complex(INF, NAN),
+                                         complex(INF, INF))
+        self.assertComplexesAreIdentical(complex(INF, NAN) * (1e300+1j),
+                                         complex(INF, INF))
+        self.assertComplexesAreIdentical(complex(INF, 1) * complex(NAN, INF),
+                                         complex(NAN, INF))
+        self.assertComplexesAreIdentical(complex(INF, 1) * complex(INF, NAN),
+                                         complex(INF, NAN))
+        self.assertComplexesAreIdentical(complex(NAN, INF) * complex(INF, 1),
+                                         complex(NAN, INF))
+        self.assertComplexesAreIdentical(complex(INF, NAN) * complex(INF, 1),
+                                         complex(INF, NAN))
+        self.assertComplexesAreIdentical(complex(NAN, 1) * complex(1, INF),
+                                         complex(-INF, NAN))
+        self.assertComplexesAreIdentical(complex(1, NAN) * complex(1, INF),
+                                         complex(NAN, INF))
+
+        self.assertComplexesAreIdentical(complex(1e200, NAN) * complex(1e200, NAN),
+                                         complex(INF, NAN))
+        self.assertComplexesAreIdentical(complex(1e200, NAN) * complex(NAN, 1e200),
+                                         complex(NAN, INF))
+        self.assertComplexesAreIdentical(complex(NAN, 1e200) * complex(1e200, NAN),
+                                         complex(NAN, INF))
+        self.assertComplexesAreIdentical(complex(NAN, 1e200) * complex(NAN, 1e200),
+                                         complex(-INF, NAN))
+
+        self.assertComplexesAreIdentical(complex(NAN, NAN) * complex(NAN, NAN),
+                                         complex(NAN, NAN))
+
     def test_mod(self):
         # % is no longer supported on complex numbers
         with self.assertRaises(TypeError):
@@ -389,6 +426,7 @@ def test_pow(self):
         self.assertAlmostEqual(pow(1j, 200), 1)
         self.assertRaises(ValueError, pow, 1+1j, 1+1j, 1+1j)
         self.assertRaises(OverflowError, pow, 1e200+1j, 1e200+1j)
+        self.assertRaises(OverflowError, pow, 1e200+1j, 5)
         self.assertRaises(TypeError, pow, 1j, None)
         self.assertRaises(TypeError, pow, None, 1j)
         self.assertAlmostEqual(pow(1j, 0.5), 0.7071067811865476+0.7071067811865475j)

Objects/complexobject.c

@@ -55,11 +55,66 @@ _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;
+
+    /* Recover infinities that computed as nan+nanj.  See e.g. the C11,
+       Annex G.5.2, routine _Cmultd(). */
+    if (isnan(r.real) && isnan(r.imag)) {
+        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;
 }