From 70d21c765f00adec71ea12913a38d6bd833339c6 Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Wed, 17 Apr 2024 19:39:42 +0300 Subject: [PATCH 1/6] gh-117999: fixed small nonnegative integer powers of complex numbers Before, handling of numbers with special values in components (infinities, nans, signed zero) was invalid. Simple example: >>> z = complex(1, -0.0) >>> z*z (1-0j) >>> z**2 (1+0j) Now: >>> z**2 (1-0j) --- Lib/test/test_complex.py | 28 +++++++++++++++++++ ...-04-17-19-39-30.gh-issue-117999.oCJkud.rst | 2 ++ Objects/complexobject.c | 11 +++----- 3 files changed, 34 insertions(+), 7 deletions(-) create mode 100644 Misc/NEWS.d/next/Core and Builtins/2024-04-17-19-39-30.gh-issue-117999.oCJkud.rst diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py index fa3017b24e16c8..d092336552a76d 100644 --- a/Lib/test/test_complex.py +++ b/Lib/test/test_complex.py @@ -6,6 +6,9 @@ from random import random from math import atan2, isnan, copysign +from cmath import log, exp, isclose, isnan as cisnan +from functools import reduce +from itertools import combinations import operator INF = float("inf") @@ -330,6 +333,31 @@ def test_pow_with_small_integer_exponents(self): self.assertEqual(str(float_pow), str(int_pow)) self.assertEqual(str(complex_pow), str(int_pow)) + # Check that complex numbers with special components + # are correctly handled. + values = [complex(*_) for _ in combinations([1, -1, 0.0, -0.0, 2, + -3, INF, -INF, NAN], 2)] + exponents = [0, 1, 2, 3, 4, 5, 6, 19] + for z in values: + for e in exponents: + try: + r_pow = z**e + except OverflowError: + continue + r_pro = reduce(lambda x, y: x*y, [z]*e) if e else 1+0j + test = str(r_pow) == str(r_pro) + if not test: + # We might fail here, because associativity of multiplication + # is broken already for floats. + # Consider z = 1-1j. Then z*z*z*z = ((z*z)*z)*z = -4+0j, + # while in the algorithm for pow() a diffenent grouping + # of operations is used: z**4 = (z*z)*(z*z) = -4-0j. + r_pro = exp(e*log(z)) + self.assertTrue(test or isclose(r_pow, r_pro)) + if not cisnan(r_pow): + self.assertEqual(copysign(1, r_pow.real), copysign(1, r_pro.real)) + self.assertEqual(copysign(1, r_pow.imag), copysign(1, r_pro.imag)) + def test_boolcontext(self): for i in range(100): self.assertTrue(complex(random() + 1e-6, random() + 1e-6)) diff --git a/Misc/NEWS.d/next/Core and Builtins/2024-04-17-19-39-30.gh-issue-117999.oCJkud.rst b/Misc/NEWS.d/next/Core and Builtins/2024-04-17-19-39-30.gh-issue-117999.oCJkud.rst new file mode 100644 index 00000000000000..0d1fc9955411b1 --- /dev/null +++ b/Misc/NEWS.d/next/Core and Builtins/2024-04-17-19-39-30.gh-issue-117999.oCJkud.rst @@ -0,0 +1,2 @@ +Fixed small nonnegative integer powers for complex numbers with special values +(``±0.0``, infinities or nan's) in components. Patch by Sergey B Kirpichev. diff --git a/Objects/complexobject.c b/Objects/complexobject.c index 943c5ccabfd5c4..52d271c39fb655 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -159,10 +159,8 @@ _Py_c_pow(Py_complex a, Py_complex b) static Py_complex c_powu(Py_complex x, long n) { - Py_complex r, p; + Py_complex p = x, r = n-- ? p : c_1; long mask = 1; - r = c_1; - p = x; while (mask > 0 && n >= mask) { if (n & mask) r = _Py_c_prod(r,p); @@ -175,11 +173,10 @@ c_powu(Py_complex x, long n) static Py_complex c_powi(Py_complex x, long n) { - if (n > 0) - return c_powu(x,n); - else + if (n < 0) return _Py_c_quot(c_1, c_powu(x,-n)); - + else + return c_powu(x,n); } double From 8e0c48264b1acba0b2d35b0639c6510e1fc63404 Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Wed, 29 May 2024 10:58:05 +0300 Subject: [PATCH 2/6] Polishing of c_powu(): exclude unreachible cases and magic numbers --- Objects/complexobject.c | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/Objects/complexobject.c b/Objects/complexobject.c index 52d271c39fb655..655a9c7001eeed 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -156,12 +156,16 @@ _Py_c_pow(Py_complex a, Py_complex b) return r; } +#define INT_EXP_CUTOFF 100 + static Py_complex c_powu(Py_complex x, long n) { Py_complex p = x, r = n-- ? p : c_1; long mask = 1; - while (mask > 0 && n >= mask) { + assert(-1 <= n && n <= INT_EXP_CUTOFF); + while (n >= mask) { + assert(mask>0); if (n & mask) r = _Py_c_prod(r,p); mask <<= 1; @@ -541,7 +545,7 @@ complex_pow(PyObject *v, PyObject *w, PyObject *z) errno = 0; // Check whether the exponent has a small integer value, and if so use // a faster and more accurate algorithm. - if (b.imag == 0.0 && b.real == floor(b.real) && fabs(b.real) <= 100.0) { + if (b.imag == 0.0 && b.real == floor(b.real) && fabs(b.real) <= INT_EXP_CUTOFF) { p = c_powi(a, (long)b.real); } else { From 6fb1d12541e2eeeb366c6d67f9b84d3d09f4ed3c Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Thu, 30 May 2024 09:50:57 +0300 Subject: [PATCH 3/6] Use _Py_c_pow() to test implementation --- Lib/test/test_complex.py | 40 +++++++++++++++++++++++++++------------- 1 file changed, 27 insertions(+), 13 deletions(-) diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py index d092336552a76d..4c2db4261c0147 100644 --- a/Lib/test/test_complex.py +++ b/Lib/test/test_complex.py @@ -6,10 +6,10 @@ from random import random from math import atan2, isnan, copysign -from cmath import log, exp, isclose, isnan as cisnan from functools import reduce from itertools import combinations import operator +import _testcapi INF = float("inf") NAN = float("nan") @@ -345,18 +345,32 @@ def test_pow_with_small_integer_exponents(self): except OverflowError: continue r_pro = reduce(lambda x, y: x*y, [z]*e) if e else 1+0j - test = str(r_pow) == str(r_pro) - if not test: - # We might fail here, because associativity of multiplication - # is broken already for floats. - # Consider z = 1-1j. Then z*z*z*z = ((z*z)*z)*z = -4+0j, - # while in the algorithm for pow() a diffenent grouping - # of operations is used: z**4 = (z*z)*(z*z) = -4-0j. - r_pro = exp(e*log(z)) - self.assertTrue(test or isclose(r_pow, r_pro)) - if not cisnan(r_pow): - self.assertEqual(copysign(1, r_pow.real), copysign(1, r_pro.real)) - self.assertEqual(copysign(1, r_pow.imag), copysign(1, r_pro.imag)) + if str(r_pow) == str(r_pro): + continue + + self.assertNotIn(z.real, {0.0, -0.0, INF, -INF, NAN}) + self.assertNotIn(z.imag, {0.0, -0.0, INF, -INF, NAN}) + + # We might fail here, because associativity of multiplication + # is broken already for floats. + # Consider z = 1-1j. Then z*z*z*z = ((z*z)*z)*z = -4+0j, + # while in the algorithm for pow() a diffenent grouping + # of operations is used: z**4 = (z*z)*(z*z) = -4-0j. + # + # Fallback to the generic complex power algorithm. + r_pro, r_pro_errno = _testcapi._py_c_pow(z, e) + self.assertEqual(r_pro_errno, 0) + self.assertClose(r_pow, r_pro) + if not isnan(r_pow.real): + self.assertEqual(copysign(1, r_pow.real), + copysign(1, r_pro.real)) + else: + self.assertTrue(isnan(r_pro.real)) + if not isnan(r_pow.imag): + self.assertEqual(copysign(1, r_pow.imag), + copysign(1, r_pro.imag)) + else: + self.assertTrue(isnan(r_pro.imag)) def test_boolcontext(self): for i in range(100): From 389441d99eff56d03efa13b969edc7e25429fb89 Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Mon, 19 Aug 2024 11:55:02 +0300 Subject: [PATCH 4/6] Apply suggestions from code review 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> --- Lib/test/test_complex.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py index bc4d94632f8f38..3e6b17d699f005 100644 --- a/Lib/test/test_complex.py +++ b/Lib/test/test_complex.py @@ -382,11 +382,11 @@ def test_pow_with_small_integer_exponents(self): r_pro, r_pro_errno = _testcapi._py_c_pow(z, e) self.assertEqual(r_pro_errno, 0) self.assertClose(r_pow, r_pro) - if not isnan(r_pow.real): + if isnan(r_pow.real): + self.assertTrue(isnan(r_pro.real)) + else: self.assertEqual(copysign(1, r_pow.real), copysign(1, r_pro.real)) - else: - self.assertTrue(isnan(r_pro.real)) if not isnan(r_pow.imag): self.assertEqual(copysign(1, r_pow.imag), copysign(1, r_pro.imag)) From 9f24ad1f9af5f68ef2cb67d2a2809a825df81dac Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Mon, 19 Aug 2024 11:57:00 +0300 Subject: [PATCH 5/6] Apply suggestions from code review 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> --- Lib/test/test_complex.py | 15 ++++++++------- Objects/complexobject.c | 2 +- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py index 3e6b17d699f005..ce08b3fe59a2f9 100644 --- a/Lib/test/test_complex.py +++ b/Lib/test/test_complex.py @@ -5,9 +5,9 @@ INVALID_UNDERSCORE_LITERALS) from random import random -from math import atan2, isnan, copysign +from math import isnan, copysign from functools import reduce -from itertools import combinations +from itertools import combinations_with_replacement import operator import _testcapi @@ -356,8 +356,9 @@ def test_pow_with_small_integer_exponents(self): # Check that complex numbers with special components # are correctly handled. - values = [complex(*_) for _ in combinations([1, -1, 0.0, -0.0, 2, - -3, INF, -INF, NAN], 2)] + values = [complex(*_) + for _ in combinations_with_replacement([1, -1, 0.0, 0, -0.0, 2, + -3, INF, -INF, NAN], 2)] exponents = [0, 1, 2, 3, 4, 5, 6, 19] for z in values: for e in exponents: @@ -387,11 +388,11 @@ def test_pow_with_small_integer_exponents(self): else: self.assertEqual(copysign(1, r_pow.real), copysign(1, r_pro.real)) - if not isnan(r_pow.imag): + if isnan(r_pow.imag): + self.assertTrue(isnan(r_pro.imag)) + else: self.assertEqual(copysign(1, r_pow.imag), copysign(1, r_pro.imag)) - else: - self.assertTrue(isnan(r_pro.imag)) def test_boolcontext(self): for i in range(100): diff --git a/Objects/complexobject.c b/Objects/complexobject.c index 3b76f94d46d65b..508ea49fefaaa5 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -163,7 +163,7 @@ c_powu(Py_complex x, long n) { Py_complex p = x, r = n-- ? p : c_1; long mask = 1; - assert(-1 <= n && n <= INT_EXP_CUTOFF); + assert(-1 <= n && n < INT_EXP_CUTOFF); while (n >= mask) { assert(mask>0); if (n & mask) From 807bcc6184b74a0ac179d06cf7f3b703bc6240c0 Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Mon, 19 Aug 2024 12:42:06 +0300 Subject: [PATCH 6/6] Address review: subTest and asserts --- Lib/test/test_complex.py | 63 ++++++++++++++++++++-------------------- Objects/complexobject.c | 3 +- 2 files changed, 34 insertions(+), 32 deletions(-) diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py index ce08b3fe59a2f9..dc5536f29d3178 100644 --- a/Lib/test/test_complex.py +++ b/Lib/test/test_complex.py @@ -362,37 +362,38 @@ def test_pow_with_small_integer_exponents(self): exponents = [0, 1, 2, 3, 4, 5, 6, 19] for z in values: for e in exponents: - try: - r_pow = z**e - except OverflowError: - continue - r_pro = reduce(lambda x, y: x*y, [z]*e) if e else 1+0j - if str(r_pow) == str(r_pro): - continue - - self.assertNotIn(z.real, {0.0, -0.0, INF, -INF, NAN}) - self.assertNotIn(z.imag, {0.0, -0.0, INF, -INF, NAN}) - - # We might fail here, because associativity of multiplication - # is broken already for floats. - # Consider z = 1-1j. Then z*z*z*z = ((z*z)*z)*z = -4+0j, - # while in the algorithm for pow() a diffenent grouping - # of operations is used: z**4 = (z*z)*(z*z) = -4-0j. - # - # Fallback to the generic complex power algorithm. - r_pro, r_pro_errno = _testcapi._py_c_pow(z, e) - self.assertEqual(r_pro_errno, 0) - self.assertClose(r_pow, r_pro) - if isnan(r_pow.real): - self.assertTrue(isnan(r_pro.real)) - else: - self.assertEqual(copysign(1, r_pow.real), - copysign(1, r_pro.real)) - if isnan(r_pow.imag): - self.assertTrue(isnan(r_pro.imag)) - else: - self.assertEqual(copysign(1, r_pow.imag), - copysign(1, r_pro.imag)) + with self.subTest(value=z, exponent=e): + try: + r_pow = z**e + except OverflowError: + continue + r_pro = reduce(lambda x, y: x*y, [z]*e) if e else 1+0j + if str(r_pow) == str(r_pro): + continue + + self.assertNotIn(z.real, {0.0, -0.0, INF, -INF, NAN}) + self.assertNotIn(z.imag, {0.0, -0.0, INF, -INF, NAN}) + + # We might fail here, because associativity of multiplication + # is broken already for floats. + # Consider z = 1-1j. Then z*z*z*z = ((z*z)*z)*z = -4+0j, + # while in the algorithm for pow() a diffenent grouping + # of operations is used: z**4 = (z*z)*(z*z) = -4-0j. + # + # Fallback to the generic complex power algorithm. + r_pro, r_pro_errno = _testcapi._py_c_pow(z, e) + self.assertEqual(r_pro_errno, 0) + self.assertClose(r_pow, r_pro) + if isnan(r_pow.real): + self.assertTrue(isnan(r_pro.real)) + else: + self.assertEqual(copysign(1, r_pow.real), + copysign(1, r_pro.real)) + if isnan(r_pow.imag): + self.assertTrue(isnan(r_pro.imag)) + else: + self.assertEqual(copysign(1, r_pow.imag), + copysign(1, r_pro.imag)) def test_boolcontext(self): for i in range(100): diff --git a/Objects/complexobject.c b/Objects/complexobject.c index 508ea49fefaaa5..43b390a7ff9013 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -163,7 +163,8 @@ c_powu(Py_complex x, long n) { Py_complex p = x, r = n-- ? p : c_1; long mask = 1; - assert(-1 <= n && n < INT_EXP_CUTOFF); + assert(-1 <= n); + assert(n < INT_EXP_CUTOFF); while (n >= mask) { assert(mask>0); if (n & mask)