From 4db428d47538b0948c79d0ba90ab29f19c3f12bd Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Tue, 4 Jun 2024 06:13:58 +0300 Subject: [PATCH] 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 "", line 1, in z**5 ~^^~ OverflowError: complex exponentiation --- Lib/test/test_complex.py | 38 +++++++++++++++++++++++++ Objects/complexobject.c | 61 ++++++++++++++++++++++++++++++++++++++-- 2 files changed, 96 insertions(+), 3 deletions(-) diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py index a3ccc9abeca09a0..26101b725cd4b7a 100644 --- a/Lib/test/test_complex.py +++ b/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) diff --git a/Objects/complexobject.c b/Objects/complexobject.c index b6703f37ce5060a..47c28b51cf6c6b6 100644 --- a/Objects/complexobject.c +++ b/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; }