From 99d1b8c017823036fe2c6bba460110cb480d38ac Mon Sep 17 00:00:00 2001
From: Sergey B Kirpichev <skirpichev@gmail.com>
Date: Tue, 21 May 2024 18:25:09 +0300
Subject: [PATCH] gh-119372: Recover inf's and zeros in
 _Py_c_quot/_Py_c_quot_real

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

>>> from cmath import inf, infj, nanj
>>> (1+1j)/(inf+infj)  # was (nan+nanj)
0j
>>> (inf+nanj)/complex(2**1000, 2**-1000)
(inf-infj)
---
 Lib/test/test_complex.py | 27 +++++++++++++++++++++++++++
 Objects/complexobject.c  | 40 ++++++++++++++++++++++++++++++++++++++--
 2 files changed, 65 insertions(+), 2 deletions(-)

diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py
index fb9d54a78fac25c..3b4577c1846b06f 100644
--- a/Lib/test/test_complex.py
+++ b/Lib/test/test_complex.py
@@ -142,6 +142,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))
+
         self.assertEqual((1+1j)/float(2), 0.5+0.5j)
         self.assertRaises(TypeError, operator.truediv, None, 1+1j)
         self.assertRaises(TypeError, operator.truediv, 1+1j, None)
diff --git a/Objects/complexobject.c b/Objects/complexobject.c
index da657229d946422..bdf2798b389a607 100644
--- a/Objects/complexobject.c
+++ b/Objects/complexobject.c
@@ -90,8 +90,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;
@@ -122,6 +121,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 ((fmax(abs_breal, abs_bimag) == Py_INFINITY)
+                 && 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;
 }
 
@@ -155,6 +176,21 @@ _Py_c_quot_real(double a, Py_complex b)
     else {
         r.real = r.imag = Py_NAN;
     }
+
+    if (isnan(r.real) && isnan(r.imag)) {
+        if (isinf(a) && isfinite(b.real) && isfinite(b.imag)) {
+            const double x = copysign(1.0, a);
+            r.real = Py_INFINITY * (x*b.real);
+            r.imag = Py_INFINITY * (-x*b.imag);
+        }
+        else if ((fmax(abs_breal, abs_bimag) == Py_INFINITY) && isfinite(a)) {
+            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*x);
+            r.imag = 0.0 * (-a*y);
+        }
+    }
+
     return r;
 }
 #ifdef _M_ARM64