From d4faa7bd321a7016f3e987d65962e02c778d708f Mon Sep 17 00:00:00 2001 From: Sergey B Kirpichev Date: Fri, 5 Jul 2024 18:01:05 +0300 Subject: [PATCH] gh-121149: improve accuracy of builtin sum() for complex inputs (gh-121176) --- Doc/library/functions.rst | 4 + Lib/test/test_builtin.py | 9 ++ ...-06-30-03-48-10.gh-issue-121149.lLBMKe.rst | 2 + Python/bltinmodule.c | 131 ++++++++++++++---- 4 files changed, 120 insertions(+), 26 deletions(-) create mode 100644 Misc/NEWS.d/next/Core and Builtins/2024-06-30-03-48-10.gh-issue-121149.lLBMKe.rst diff --git a/Doc/library/functions.rst b/Doc/library/functions.rst index 1d82f92ea67..17348dd907b 100644 --- a/Doc/library/functions.rst +++ b/Doc/library/functions.rst @@ -1934,6 +1934,10 @@ are always available. They are listed here in alphabetical order. .. versionchanged:: 3.12 Summation of floats switched to an algorithm that gives higher accuracy and better commutativity on most builds. + .. versionchanged:: 3.14 + Added specialization for summation of complexes, + using same algorithm as for summation of floats. + .. class:: super() super(type, object_or_type=None) diff --git a/Lib/test/test_builtin.py b/Lib/test/test_builtin.py index 9ff0f488dc4..5818e96d61f 100644 --- a/Lib/test/test_builtin.py +++ b/Lib/test/test_builtin.py @@ -1768,6 +1768,11 @@ class BuiltinTest(unittest.TestCase): sum(([x] for x in range(10)), empty) self.assertEqual(empty, []) + xs = [complex(random.random() - .5, random.random() - .5) + for _ in range(10000)] + self.assertEqual(sum(xs), complex(sum(z.real for z in xs), + sum(z.imag for z in xs))) + @requires_IEEE_754 @unittest.skipIf(HAVE_DOUBLE_ROUNDING, "sum accuracy not guaranteed on machines with double rounding") @@ -1775,6 +1780,10 @@ class BuiltinTest(unittest.TestCase): def test_sum_accuracy(self): self.assertEqual(sum([0.1] * 10), 1.0) self.assertEqual(sum([1.0, 10E100, 1.0, -10E100]), 2.0) + self.assertEqual(sum([1.0, 10E100, 1.0, -10E100, 2j]), 2+2j) + self.assertEqual(sum([2+1j, 10E100j, 1j, -10E100j]), 2+2j) + self.assertEqual(sum([1j, 1, 10E100j, 1j, 1.0, -10E100j]), 2+2j) + self.assertEqual(sum([0.1j]*10 + [fractions.Fraction(1, 10)]), 0.1+1j) def test_type(self): self.assertEqual(type(''), type('123')) diff --git a/Misc/NEWS.d/next/Core and Builtins/2024-06-30-03-48-10.gh-issue-121149.lLBMKe.rst b/Misc/NEWS.d/next/Core and Builtins/2024-06-30-03-48-10.gh-issue-121149.lLBMKe.rst new file mode 100644 index 00000000000..38d618f0609 --- /dev/null +++ b/Misc/NEWS.d/next/Core and Builtins/2024-06-30-03-48-10.gh-issue-121149.lLBMKe.rst @@ -0,0 +1,2 @@ +Added specialization for summation of complexes, this also improves accuracy +of builtin :func:`sum` for such inputs. Patch by Sergey B Kirpichev. diff --git a/Python/bltinmodule.c b/Python/bltinmodule.c index 6e50623cafa..a5b45e358d9 100644 --- a/Python/bltinmodule.c +++ b/Python/bltinmodule.c @@ -2516,6 +2516,49 @@ Without arguments, equivalent to locals().\n\ With an argument, equivalent to object.__dict__."); +/* Improved Kahan–Babuška algorithm by Arnold Neumaier + Neumaier, A. (1974), Rundungsfehleranalyse einiger Verfahren + zur Summation endlicher Summen. Z. angew. Math. Mech., + 54: 39-51. https://doi.org/10.1002/zamm.19740540106 + https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements + */ + +typedef struct { + double hi; /* high-order bits for a running sum */ + double lo; /* a running compensation for lost low-order bits */ +} CompensatedSum; + +static inline CompensatedSum +cs_from_double(double x) +{ + return (CompensatedSum) {x}; +} + +static inline CompensatedSum +cs_add(CompensatedSum total, double x) +{ + double t = total.hi + x; + if (fabs(total.hi) >= fabs(x)) { + total.lo += (total.hi - t) + x; + } + else { + total.lo += (x - t) + total.hi; + } + return (CompensatedSum) {t, total.lo}; +} + +static inline double +cs_to_double(CompensatedSum total) +{ + /* Avoid losing the sign on a negative result, + and don't let adding the compensation convert + an infinite or overflowed sum to a NaN. */ + if (total.lo && isfinite(total.lo)) { + return total.hi + total.lo; + } + return total.hi; +} + /*[clinic input] sum as builtin_sum @@ -2628,8 +2671,7 @@ builtin_sum_impl(PyObject *module, PyObject *iterable, PyObject *start) } if (PyFloat_CheckExact(result)) { - double f_result = PyFloat_AS_DOUBLE(result); - double c = 0.0; + CompensatedSum re_sum = cs_from_double(PyFloat_AS_DOUBLE(result)); Py_SETREF(result, NULL); while(result == NULL) { item = PyIter_Next(iter); @@ -2637,28 +2679,10 @@ builtin_sum_impl(PyObject *module, PyObject *iterable, PyObject *start) Py_DECREF(iter); if (PyErr_Occurred()) return NULL; - /* Avoid losing the sign on a negative result, - and don't let adding the compensation convert - an infinite or overflowed sum to a NaN. */ - if (c && isfinite(c)) { - f_result += c; - } - return PyFloat_FromDouble(f_result); + return PyFloat_FromDouble(cs_to_double(re_sum)); } if (PyFloat_CheckExact(item)) { - // Improved Kahan–Babuška algorithm by Arnold Neumaier - // Neumaier, A. (1974), Rundungsfehleranalyse einiger Verfahren - // zur Summation endlicher Summen. Z. angew. Math. Mech., - // 54: 39-51. https://doi.org/10.1002/zamm.19740540106 - // https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements - double x = PyFloat_AS_DOUBLE(item); - double t = f_result + x; - if (fabs(f_result) >= fabs(x)) { - c += (f_result - t) + x; - } else { - c += (x - t) + f_result; - } - f_result = t; + re_sum = cs_add(re_sum, PyFloat_AS_DOUBLE(item)); _Py_DECREF_SPECIALIZED(item, _PyFloat_ExactDealloc); continue; } @@ -2667,15 +2691,70 @@ builtin_sum_impl(PyObject *module, PyObject *iterable, PyObject *start) int overflow; value = PyLong_AsLongAndOverflow(item, &overflow); if (!overflow) { - f_result += (double)value; + re_sum.hi += (double)value; Py_DECREF(item); continue; } } - if (c && isfinite(c)) { - f_result += c; + result = PyFloat_FromDouble(cs_to_double(re_sum)); + if (result == NULL) { + Py_DECREF(item); + Py_DECREF(iter); + return NULL; } - result = PyFloat_FromDouble(f_result); + temp = PyNumber_Add(result, item); + Py_DECREF(result); + Py_DECREF(item); + result = temp; + if (result == NULL) { + Py_DECREF(iter); + return NULL; + } + } + } + + if (PyComplex_CheckExact(result)) { + Py_complex z = PyComplex_AsCComplex(result); + CompensatedSum re_sum = cs_from_double(z.real); + CompensatedSum im_sum = cs_from_double(z.imag); + Py_SETREF(result, NULL); + while (result == NULL) { + item = PyIter_Next(iter); + if (item == NULL) { + Py_DECREF(iter); + if (PyErr_Occurred()) { + return NULL; + } + return PyComplex_FromDoubles(cs_to_double(re_sum), + cs_to_double(im_sum)); + } + if (PyComplex_CheckExact(item)) { + z = PyComplex_AsCComplex(item); + re_sum = cs_add(re_sum, z.real); + im_sum = cs_add(im_sum, z.imag); + Py_DECREF(item); + continue; + } + if (PyLong_Check(item)) { + long value; + int overflow; + value = PyLong_AsLongAndOverflow(item, &overflow); + if (!overflow) { + re_sum.hi += (double)value; + im_sum.hi += 0.0; + Py_DECREF(item); + continue; + } + } + if (PyFloat_Check(item)) { + double value = PyFloat_AS_DOUBLE(item); + re_sum.hi += value; + im_sum.hi += 0.0; + _Py_DECREF_SPECIALIZED(item, _PyFloat_ExactDealloc); + continue; + } + result = PyComplex_FromDoubles(cs_to_double(re_sum), + cs_to_double(im_sum)); if (result == NULL) { Py_DECREF(item); Py_DECREF(iter);