https://github.com/python/cpython/commit/d4faa7bd321a7016f3e987d65962e02c778d708f
commit: d4faa7bd321a7016f3e987d65962e02c778d708f
branch: main
author: Sergey B Kirpichev <[email protected]>
committer: rhettinger <[email protected]>
date: 2024-07-05T10:01:05-05:00
summary:
gh-121149: improve accuracy of builtin sum() for complex inputs (gh-121176)
files:
A Misc/NEWS.d/next/Core and
Builtins/2024-06-30-03-48-10.gh-issue-121149.lLBMKe.rst
M Doc/library/functions.rst
M Lib/test/test_builtin.py
M Python/bltinmodule.c
diff --git a/Doc/library/functions.rst b/Doc/library/functions.rst
index 1d82f92ea67857..17348dd907bf67 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 9ff0f488dc4fa9..5818e96d61f480 100644
--- a/Lib/test/test_builtin.py
+++ b/Lib/test/test_builtin.py
@@ -1768,6 +1768,11 @@ def __getitem__(self, index):
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 @@ def __getitem__(self, index):
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 00000000000000..38d618f06090fd
--- /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 6e50623cafa4ed..a5b45e358d9efb 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;
+ }
+ }
+ result = PyFloat_FromDouble(cs_to_double(re_sum));
+ if (result == NULL) {
+ Py_DECREF(item);
+ Py_DECREF(iter);
+ return NULL;
+ }
+ 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 (c && isfinite(c)) {
- f_result += c;
+ 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 = PyFloat_FromDouble(f_result);
+ result = PyComplex_FromDoubles(cs_to_double(re_sum),
+ cs_to_double(im_sum));
if (result == NULL) {
Py_DECREF(item);
Py_DECREF(iter);
_______________________________________________
Python-checkins mailing list -- [email protected]
To unsubscribe send an email to [email protected]
https://mail.python.org/mailman3/lists/python-checkins.python.org/
Member address: [email protected]
