Kahan summation

rhettinger · Dec 22, 2022 · 787f2e5 · 787f2e5
1 parent 916ab7e
commit 787f2e5
Show file tree

Hide file tree

Showing 2 changed files with 10 additions and 12 deletions.
diff --git 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
@@ -2533,25 +2533,22 @@ 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);
             if (item == NULL) {
                 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;