gh-119372: Recover inf's and zeros in _Py_c_quot (GH-119457)

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 _Cdivd().

Author: skirpichev

Lib/test/test_complex.py

@@ -94,6 +94,10 @@ def assertFloatsAreIdentical(self, x, y):
                 msg += ': zeros have different signs'
         self.fail(msg.format(x, y))
 
+    def assertComplexesAreIdentical(self, x, y):
+        self.assertFloatsAreIdentical(x.real, y.real)
+        self.assertFloatsAreIdentical(x.imag, y.imag)
+
     def assertClose(self, x, y, eps=1e-9):
         """Return true iff complexes x and y "are close"."""
         self.assertCloseAbs(x.real, y.real, eps)
@@ -139,6 +143,33 @@ def test_truediv(self):
             self.assertTrue(isnan(z.real))
             self.assertTrue(isnan(z.imag))
 
+        self.assertComplexesAreIdentical(complex(INF, 1)/(0.0+1j),
+                                         complex(NAN, -INF))
+
+        # test recover of infs if numerator has infs and denominator is finite
+        self.assertComplexesAreIdentical(complex(INF, -INF)/(1+0j),
+                                         complex(INF, -INF))
+        self.assertComplexesAreIdentical(complex(INF, INF)/(0.0+1j),
+                                         complex(INF, -INF))
+        self.assertComplexesAreIdentical(complex(NAN, INF)/complex(2**1000, 2**-1000),
+                                         complex(INF, INF))
+        self.assertComplexesAreIdentical(complex(INF, NAN)/complex(2**1000, 2**-1000),
+                                         complex(INF, -INF))
+
+        # test recover of zeros if denominator is infinite
+        self.assertComplexesAreIdentical((1+1j)/complex(INF, INF), (0.0+0j))
+        self.assertComplexesAreIdentical((1+1j)/complex(INF, -INF), (0.0+0j))
+        self.assertComplexesAreIdentical((1+1j)/complex(-INF, INF),
+                                         complex(0.0, -0.0))
+        self.assertComplexesAreIdentical((1+1j)/complex(-INF, -INF),
+                                         complex(-0.0, 0))
+        self.assertComplexesAreIdentical((INF+1j)/complex(INF, INF),
+                                         complex(NAN, NAN))
+        self.assertComplexesAreIdentical(complex(1, INF)/complex(INF, INF),
+                                         complex(NAN, NAN))
+        self.assertComplexesAreIdentical(complex(INF, 1)/complex(1, INF),
+                                         complex(NAN, NAN))
+
     def test_truediv_zero_division(self):
         for a, b in ZERO_DIVISION:
             with self.assertRaises(ZeroDivisionError):

Misc/NEWS.d/next/Core and Builtins/2024-05-22-12-49-03.gh-issue-119372.PXig1R.rst

@@ -0,0 +1,2 @@
+Correct invalid corner cases in complex division (resulted in ``(nan+nanj)``
+output), e.g.  ``1/complex('(inf+infj)')``.  Patch by Sergey B Kirpichev.

Objects/complexobject.c

@@ -88,8 +88,7 @@ _Py_c_quot(Py_complex a, Py_complex b)
      * numerators and denominator by whichever of {b.real, b.imag} has
      * larger magnitude.  The earliest reference I found was to CACM
      * Algorithm 116 (Complex Division, Robert L. Smith, Stanford
-     * University).  As usual, though, we're still ignoring all IEEE
-     * endcases.
+     * University).
      */
      Py_complex r;      /* the result */
      const double abs_breal = b.real < 0 ? -b.real : b.real;
@@ -120,6 +119,28 @@ _Py_c_quot(Py_complex a, Py_complex b)
         /* At least one of b.real or b.imag is a NaN */
         r.real = r.imag = Py_NAN;
     }
+
+    /* Recover infinities and zeros that computed as nan+nanj.  See e.g.
+       the C11, Annex G.5.2, routine _Cdivd(). */
+    if (isnan(r.real) && isnan(r.imag)) {
+        if ((isinf(a.real) || isinf(a.imag))
+            && isfinite(b.real) && isfinite(b.imag))
+        {
+            const double x = copysign(isinf(a.real) ? 1.0 : 0.0, a.real);
+            const double y = copysign(isinf(a.imag) ? 1.0 : 0.0, a.imag);
+            r.real = Py_INFINITY * (x*b.real + y*b.imag);
+            r.imag = Py_INFINITY * (y*b.real - x*b.imag);
+        }
+        else if ((isinf(abs_breal) || isinf(abs_bimag))
+                 && isfinite(a.real) && isfinite(a.imag))
+        {
+            const double x = copysign(isinf(b.real) ? 1.0 : 0.0, b.real);
+            const double y = copysign(isinf(b.imag) ? 1.0 : 0.0, b.imag);
+            r.real = 0.0 * (a.real*x + a.imag*y);
+            r.imag = 0.0 * (a.imag*x - a.real*y);
+        }
+    }
+
     return r;
 }
 #ifdef _M_ARM64