From 787f2e52b94f5501bf933a54f4720ef1c4784b00 Mon Sep 17 00:00:00 2001 From: Raymond Hettinger Date: Wed, 21 Dec 2022 17:21:55 -0800 Subject: [PATCH] Kahan summation --- Lib/test/test_builtin.py | 7 ++++--- Python/bltinmodule.c | 15 ++++++--------- 2 files changed, 10 insertions(+), 12 deletions(-) diff --git a/Lib/test/test_builtin.py b/Lib/test/test_builtin.py index d4708e5a3b50b1..23e0567575e116 100644 --- a/Lib/test/test_builtin.py +++ b/Lib/test/test_builtin.py @@ -1614,9 +1614,10 @@ def test_sum(self): self.assertEqual(sum((i / 2 for i in range(10)), 1000.25), 1022.75) self.assertEqual(sum([0.5, 1]), 1.5) self.assertEqual(sum([1, 0.5]), 1.5) - # self.assertEqual(repr(sum([-0.0])), '0.0') - # self.assertEqual(repr(sum([-0.0], -0.0)), '-0.0') - # self.assertEqual(repr(sum([], -0.0)), '-0.0') + self.assertEqual(sum([0.1] * 10), 1.0) # Test accuracy of Kahan summation + self.assertEqual(repr(sum([-0.0])), '0.0') + self.assertEqual(repr(sum([-0.0], -0.0)), '-0.0') + self.assertEqual(repr(sum([], -0.0)), '-0.0') self.assertRaises(TypeError, sum) self.assertRaises(TypeError, sum, 42) diff --git a/Python/bltinmodule.c b/Python/bltinmodule.c index 1681b0f588b5b4..87feed299d8789 100644 --- a/Python/bltinmodule.c +++ b/Python/bltinmodule.c @@ -2533,7 +2533,6 @@ builtin_sum_impl(PyObject *module, PyObject *iterable, PyObject *start) if (PyFloat_CheckExact(result)) { double f_result = PyFloat_AS_DOUBLE(result); double c = 0.0; - double x, t; Py_SETREF(result, NULL); while(result == NULL) { item = PyIter_Next(iter); @@ -2541,17 +2540,15 @@ builtin_sum_impl(PyObject *module, PyObject *iterable, PyObject *start) Py_DECREF(iter); if (PyErr_Occurred()) return NULL; - return PyFloat_FromDouble(f_result + c); + return PyFloat_FromDouble(f_result); } if (PyFloat_CheckExact(item)) { - // Neumaier compensated summation + // Kahan compensated summation + double x, y, t; x = PyFloat_AS_DOUBLE(item); - t = f_result + x; - if (fabs(f_result) >= fabs(x)) { - c += (f_result - t) + x; - } else { - c += (x - t) + f_result; - } + y = x - c; + t = f_result + y; + c = (t - f_result) - y; f_result = t; _Py_DECREF_SPECIALIZED(item, _PyFloat_ExactDealloc); continue;