From 399e3e37f4fdd144e6d2a8068fb811103ba4bc84 Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Tue, 4 Jun 2024 06:13:58 +0300 Subject: [PATCH 1/7] gh-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 "", line 1, in z**5 ~^^~ OverflowError: complex exponentiation --- Lib/test/test_complex.py | 42 +++++++++++++ ...-06-04-08-26-25.gh-issue-120010._z-AWz.rst | 2 + Objects/complexobject.c | 61 ++++++++++++++++++- 3 files changed, 102 insertions(+), 3 deletions(-) create mode 100644 Misc/NEWS.d/next/Core and Builtins/2024-06-04-08-26-25.gh-issue-120010._z-AWz.rst diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py index fb510ca9b70902..91642aaad4f9bd 100644 --- a/Lib/test/test_complex.py +++ b/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) @@ -227,6 +231,43 @@ def test_mul(self): self.assertRaises(TypeError, operator.mul, 1j, None) self.assertRaises(TypeError, operator.mul, None, 1j) + 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): @@ -268,6 +309,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) diff --git a/Misc/NEWS.d/next/Core and Builtins/2024-06-04-08-26-25.gh-issue-120010._z-AWz.rst b/Misc/NEWS.d/next/Core and Builtins/2024-06-04-08-26-25.gh-issue-120010._z-AWz.rst new file mode 100644 index 00000000000000..45e08b98c6883c --- /dev/null +++ b/Misc/NEWS.d/next/Core and Builtins/2024-06-04-08-26-25.gh-issue-120010._z-AWz.rst @@ -0,0 +1,2 @@ +Correct invalid corner cases (resulted in ``(nan+nanj)`` output) in complex +multiplication, e.g. ``(1e300+1j)*(nan+infj)``. Patch by Sergey B Kirpichev. diff --git a/Objects/complexobject.c b/Objects/complexobject.c index a8be266970afd0..c5f40b096552b8 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -53,11 +53,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; } From 324c6022d4719cec78bcf4a7003f2988cb17987f Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Mon, 10 Jun 2024 03:05:59 +0300 Subject: [PATCH 2/7] address review: declarations --- Objects/complexobject.c | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/Objects/complexobject.c b/Objects/complexobject.c index c5f40b096552b8..e201cd06c31daf 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -55,12 +55,9 @@ _Py_c_neg(Py_complex a) Py_complex _Py_c_prod(Py_complex z, Py_complex w) { - Py_complex r; 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; + Py_complex r = {ac - bd, ad + bc}; /* Recover infinities that computed as nan+nanj. See e.g. the C11, Annex G.5.2, routine _Cmultd(). */ From 5509a2d209206b9481794015f6a9109fbb70dee0 Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Fri, 27 Sep 2024 13:30:01 +0300 Subject: [PATCH 3/7] address review: use subTest --- Lib/test/test_complex.py | 51 ++++++++++++---------------------------- 1 file changed, 15 insertions(+), 36 deletions(-) diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py index cd499a49993eaa..a43c02351fa57b 100644 --- a/Lib/test/test_complex.py +++ b/Lib/test/test_complex.py @@ -242,42 +242,21 @@ def test_mul(self): self.assertRaises(TypeError, operator.mul, 1j, None) self.assertRaises(TypeError, operator.mul, None, 1j) - 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)) + for z, w, r in [(1e300+1j, complex(INF, INF), complex(NAN, INF)), + (1e300+1j, complex(NAN, INF), complex(-INF, INF)), + (1e300+1j, complex(INF, NAN), complex(INF, INF)), + (complex(INF, 1), complex(NAN, INF), complex(NAN, INF)), + (complex(INF, 1), complex(INF, NAN), complex(INF, NAN)), + (complex(NAN, 1), complex(1, INF), complex(-INF, NAN)), + (complex(1, NAN), complex(1, INF), complex(NAN, INF)), + (complex(1e200, NAN), complex(1e200, NAN), complex(INF, NAN)), + (complex(1e200, NAN), complex(NAN, 1e200), complex(NAN, INF)), + (complex(NAN, 1e200), complex(1e200, NAN), complex(NAN, INF)), + (complex(NAN, 1e200), complex(NAN, 1e200), complex(-INF, NAN)), + (complex(NAN, NAN), complex(NAN, NAN), complex(NAN, NAN))]: + with self.subTest(z=z, w=w, r=r): + self.assertComplexesAreIdentical(z * w, s) + self.assertComplexesAreIdentical(w * z, s) def test_mod(self): # % is no longer supported on complex numbers From 0e9811677417d4339fab0d4d01013e0c5b33d81a Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Fri, 27 Sep 2024 14:43:56 +0300 Subject: [PATCH 4/7] Update Misc/NEWS.d/next/Core and Builtins/2024-06-04-08-26-25.gh-issue-120010._z-AWz.rst MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Bénédikt Tran <10796600+picnixz@users.noreply.github.com> --- .../2024-06-04-08-26-25.gh-issue-120010._z-AWz.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Misc/NEWS.d/next/Core and Builtins/2024-06-04-08-26-25.gh-issue-120010._z-AWz.rst b/Misc/NEWS.d/next/Core and Builtins/2024-06-04-08-26-25.gh-issue-120010._z-AWz.rst index 45e08b98c6883c..7954c7f5927397 100644 --- a/Misc/NEWS.d/next/Core and Builtins/2024-06-04-08-26-25.gh-issue-120010._z-AWz.rst +++ b/Misc/NEWS.d/next/Core and Builtins/2024-06-04-08-26-25.gh-issue-120010._z-AWz.rst @@ -1,2 +1,2 @@ -Correct invalid corner cases (resulted in ``(nan+nanj)`` output) in complex -multiplication, e.g. ``(1e300+1j)*(nan+infj)``. Patch by Sergey B Kirpichev. +Correct invalid corner cases which resulted in ``(nan+nanj)`` output in complex +multiplication, e.g., ``(1e300+1j)*(nan+infj)``. Patch by Sergey B Kirpichev. From 338176d66160bbbd021f6b11d3d5299e8358bb34 Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Fri, 27 Sep 2024 15:21:55 +0300 Subject: [PATCH 5/7] + fix typo --- Lib/test/test_complex.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py index a43c02351fa57b..df80a7cda7f65b 100644 --- a/Lib/test/test_complex.py +++ b/Lib/test/test_complex.py @@ -255,8 +255,8 @@ def test_mul(self): (complex(NAN, 1e200), complex(NAN, 1e200), complex(-INF, NAN)), (complex(NAN, NAN), complex(NAN, NAN), complex(NAN, NAN))]: with self.subTest(z=z, w=w, r=r): - self.assertComplexesAreIdentical(z * w, s) - self.assertComplexesAreIdentical(w * z, s) + self.assertComplexesAreIdentical(z * w, r) + self.assertComplexesAreIdentical(w * z, r) def test_mod(self): # % is no longer supported on complex numbers From 219c0cb331c84f95bf175a2a4e1ac1f9d6c50c9e Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Wed, 27 Nov 2024 05:22:52 +0300 Subject: [PATCH 6/7] + rename --- .../2024-06-04-08-26-25.gh-issue-120010._z-AWz.rst | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename Misc/NEWS.d/next/{Core and Builtins => Core_and_Builtins}/2024-06-04-08-26-25.gh-issue-120010._z-AWz.rst (100%) diff --git a/Misc/NEWS.d/next/Core and Builtins/2024-06-04-08-26-25.gh-issue-120010._z-AWz.rst b/Misc/NEWS.d/next/Core_and_Builtins/2024-06-04-08-26-25.gh-issue-120010._z-AWz.rst similarity index 100% rename from Misc/NEWS.d/next/Core and Builtins/2024-06-04-08-26-25.gh-issue-120010._z-AWz.rst rename to Misc/NEWS.d/next/Core_and_Builtins/2024-06-04-08-26-25.gh-issue-120010._z-AWz.rst From e368f23bb16e4f88e14048068f56efcf015dc340 Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Fri, 6 Dec 2024 12:27:07 +0300 Subject: [PATCH 7/7] Update Objects/complexobject.c Co-authored-by: Serhiy Storchaka --- Objects/complexobject.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Objects/complexobject.c b/Objects/complexobject.c index e299344af6980c..bf6187efac941f 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -92,7 +92,7 @@ _Py_c_prod(Py_complex z, Py_complex w) Py_complex r = {ac - bd, ad + bc}; /* Recover infinities that computed as nan+nanj. See e.g. the C11, - Annex G.5.2, routine _Cmultd(). */ + Annex G.5.1, routine _Cmultd(). */ if (isnan(r.real) && isnan(r.imag)) { int recalc = 0;