From b60242ef21611b7bf487b45065d3e9c9e45fa9b3 Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Sat, 24 Aug 2024 05:56:04 +0300 Subject: [PATCH 1/4] gh-117999: drop special case for small integer exponents --- Lib/test/test_complex.py | 5 +++ ...-08-24-10-46-35.gh-issue-117999.5K_BiA.rst | 2 + Objects/complexobject.c | 38 +------------------ 3 files changed, 8 insertions(+), 37 deletions(-) create mode 100644 Misc/NEWS.d/next/Library/2024-08-24-10-46-35.gh-issue-117999.5K_BiA.rst diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py index c5a06c53771fca..6a9e039401d3d1 100644 --- a/Lib/test/test_complex.py +++ b/Lib/test/test_complex.py @@ -364,6 +364,11 @@ def test_pow(self): except OverflowError: pass + # Check that complex numbers with special components + # are correctly handled. + self.assertComplexesAreIdentical(complex(1, +0.0)**2, complex(1, +0.0)) + self.assertComplexesAreIdentical(complex(1, -0.0)**2, complex(1, -0.0)) + def test_pow_with_small_integer_exponents(self): # Check that small integer exponents are handled identically # regardless of their type. diff --git a/Misc/NEWS.d/next/Library/2024-08-24-10-46-35.gh-issue-117999.5K_BiA.rst b/Misc/NEWS.d/next/Library/2024-08-24-10-46-35.gh-issue-117999.5K_BiA.rst new file mode 100644 index 00000000000000..d7621bf6e5b6ca --- /dev/null +++ b/Misc/NEWS.d/next/Library/2024-08-24-10-46-35.gh-issue-117999.5K_BiA.rst @@ -0,0 +1,2 @@ +Use single algorithm for complex power (case of small integer exponents was +handled before by a special method). Patch by Sergey B Kirpichev. diff --git a/Objects/complexobject.c b/Objects/complexobject.c index 4a8dac6c53f529..0acce988625651 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -23,8 +23,6 @@ class complex "PyComplexObject *" "&PyComplex_Type" /* elementary operations on complex numbers */ -static Py_complex c_1 = {1., 0.}; - Py_complex _Py_c_sum(Py_complex a, Py_complex b) { @@ -177,32 +175,6 @@ _Py_c_pow(Py_complex a, Py_complex b) return r; } -static Py_complex -c_powu(Py_complex x, long n) -{ - Py_complex r, p; - long mask = 1; - r = c_1; - p = x; - while (mask > 0 && n >= mask) { - if (n & mask) - r = _Py_c_prod(r,p); - mask <<= 1; - p = _Py_c_prod(p,p); - } - return r; -} - -static Py_complex -c_powi(Py_complex x, long n) -{ - if (n > 0) - return c_powu(x,n); - else - return _Py_c_quot(c_1, c_powu(x,-n)); - -} - double _Py_c_abs(Py_complex z) { @@ -563,15 +535,7 @@ complex_pow(PyObject *v, PyObject *w, PyObject *z) return NULL; } 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) { - p = c_powi(a, (long)b.real); - } - else { - p = _Py_c_pow(a, b); - } - + p = _Py_c_pow(a, b); _Py_ADJUST_ERANGE2(p.real, p.imag); if (errno == EDOM) { PyErr_SetString(PyExc_ZeroDivisionError, From 2d32c4d7cc728a6f28a9e8ec01c3d40ac591906b Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Sat, 24 Aug 2024 16:31:53 +0300 Subject: [PATCH 2/4] 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> --- .../Library/2024-08-24-10-46-35.gh-issue-117999.5K_BiA.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Misc/NEWS.d/next/Library/2024-08-24-10-46-35.gh-issue-117999.5K_BiA.rst b/Misc/NEWS.d/next/Library/2024-08-24-10-46-35.gh-issue-117999.5K_BiA.rst index d7621bf6e5b6ca..2e3c45ad97215f 100644 --- a/Misc/NEWS.d/next/Library/2024-08-24-10-46-35.gh-issue-117999.5K_BiA.rst +++ b/Misc/NEWS.d/next/Library/2024-08-24-10-46-35.gh-issue-117999.5K_BiA.rst @@ -1,2 +1,2 @@ -Use single algorithm for complex power (case of small integer exponents was -handled before by a special method). Patch by Sergey B Kirpichev. +Use a single algorithm for complex exponentiation (the case where the exponent +is a small integer was previously handled separately). Patch by Sergey B Kirpichev. From 8c30627eab14fffcc0b32ef635420d848f1e1ef0 Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Mon, 20 Jan 2025 15:04:18 +0300 Subject: [PATCH 3/4] Skip integer algorithm if the power base has special components --- Lib/test/test_complex.py | 38 +++++++++++++++++++++++++++---- Objects/complexobject.c | 49 +++++++++++++++++++++++++++++++++++++++- 2 files changed, 81 insertions(+), 6 deletions(-) diff --git a/Lib/test/test_complex.py b/Lib/test/test_complex.py index 38b29627a26986..0720d6ee036714 100644 --- a/Lib/test/test_complex.py +++ b/Lib/test/test_complex.py @@ -1,3 +1,4 @@ +import cmath import unittest import sys from test import support @@ -7,7 +8,9 @@ from random import random from math import isnan, copysign +from itertools import combinations_with_replacement import operator +import _testcapi INF = float("inf") NAN = float("nan") @@ -412,11 +415,6 @@ def test_pow(self): except OverflowError: pass - # Check that complex numbers with special components - # are correctly handled. - self.assertComplexesAreIdentical(complex(1, +0.0)**2, complex(1, +0.0)) - self.assertComplexesAreIdentical(complex(1, -0.0)**2, complex(1, -0.0)) - # gh-113841: possible undefined division by 0 in _Py_c_pow() x, y = 9j, 33j**3 with self.assertRaises(OverflowError): @@ -450,6 +448,36 @@ 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_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: + with self.subTest(value=z, exponent=e): + if cmath.isfinite(z) and z.real and z.imag: + continue + try: + r_pow = z**e + except OverflowError: + continue + # Use the generic complex power algorithm. + r_pro, r_pro_errno = _testcapi._py_c_pow(z, e) + self.assertEqual(r_pro_errno, 0) + 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): self.assertTrue(complex(random() + 1e-6, random() + 1e-6)) diff --git a/Objects/complexobject.c b/Objects/complexobject.c index 82e707e4bc9c6b..f9ec107ad49baa 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -26,6 +26,8 @@ class complex "PyComplexObject *" "&PyComplex_Type" /* elementary operations on complex numbers */ +static Py_complex c_1 = {1., 0.}; + Py_complex _Py_c_sum(Py_complex a, Py_complex b) { @@ -336,6 +338,38 @@ _Py_c_pow(Py_complex a, Py_complex b) return r; } +/* Switch to exponentiation by squaring if integer exponent less that this. */ +#define INT_EXP_CUTOFF 100 + +static Py_complex +c_powu(Py_complex x, long n) +{ + Py_complex r, p; + long mask = 1; + r = c_1; + p = x; + assert(0 <= n); + assert(n <= C_EXP_CUTOFF); + while (n >= mask) { + assert(mask > 0); + if (n & mask) + r = _Py_c_prod(r,p); + mask <<= 1; + p = _Py_c_prod(p,p); + } + return r; +} + +static Py_complex +c_powi(Py_complex x, long n) +{ + if (n > 0) + return c_powu(x,n); + else + return _Py_c_quot(c_1, c_powu(x,-n)); + +} + double _Py_c_abs(Py_complex z) { @@ -705,7 +739,20 @@ complex_pow(PyObject *v, PyObject *w, PyObject *z) return NULL; } errno = 0; - p = _Py_c_pow(a, b); + // Check whether the exponent has a small integer value, and if so use + // a faster and more accurate algorithm. + // Fallback on the generic code if the base has special + // components (zeros or infinities). + if (b.imag == 0.0 && b.real == floor(b.real) && fabs(b.real) <= INT_EXP_CUTOFF + && isfinite(a.real) && a.real && isfinite(a.imag) && a.imag) + { + p = c_powi(a, (long)b.real); + _Py_ADJUST_ERANGE2(p.real, p.imag); + } + else { + p = _Py_c_pow(a, b); + } + if (errno == EDOM) { PyErr_SetString(PyExc_ZeroDivisionError, "zero to a negative or complex power"); From 4f2ce53147ca71f27597f108cd88318c3e69ad8c Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Mon, 20 Jan 2025 15:08:47 +0300 Subject: [PATCH 4/4] +1 --- Objects/complexobject.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Objects/complexobject.c b/Objects/complexobject.c index f9ec107ad49baa..9d50898cf3e178 100644 --- a/Objects/complexobject.c +++ b/Objects/complexobject.c @@ -349,7 +349,7 @@ c_powu(Py_complex x, long n) r = c_1; p = x; assert(0 <= n); - assert(n <= C_EXP_CUTOFF); + assert(n <= INT_EXP_CUTOFF); while (n >= mask) { assert(mask > 0); if (n & mask)