diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c
index 45d842f..5c5def2 100644
--- a/Modules/mathmodule.c
+++ b/Modules/mathmodule.c
@@ -362,6 +362,199 @@
 FUNC1(tanh, tanh, 0,
       "tanh(x)\n\nReturn the hyperbolic tangent of x.")
 
+/* Precision summation function as msum() by Raymond Hettinger in
+   <http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>,
+   enhanced with the exact partials sum and roundoff from Mark
+   Dickinson's post at <http://bugs.python.org/file10357/msum4.py>.
+   See those links for more details, proofs and other references.
+
+   Note 1: IEEE 754R floating point semantics are assumed,
+   but the current implementation does not re-establish special
+   value semantics across iterations (i.e. handling -Inf + Inf).
+
+   Note 2:  No provision is made for intermediate overflow handling;
+   therefore, sum([1e+308, 1e-308, 1e+308]) returns result 1e+308 while
+   sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the
+   overflow of the first partial sum.
+
+   Note 3: Aggressively optimizing compilers can potentially eliminate the
+   residual values needed for accurate summation. For instance, the statements
+   "hi = x + y; lo = y - (hi - x);" could be mis-transformed to
+   "hi = x + y; lo = 0.0;" which defeats the computation of residuals.
+
+   Note 4: A similar implementation is in Modules/cmathmodule.c.
+   Be sure to update both when making changes.
+
+   Note 5: The signature of math.sum() differs from __builtin__.sum()
+   because the start argument doesn't make sense in the context of
+   accurate summation.  Since the partials table is collapsed before
+   returning a result, sum(seq2, start=sum(seq1)) may not equal the
+   accurate result returned by sum(itertools.chain(seq1, seq2)).
+*/
+
+#define NUM_PARTIALS  32  /* initial partials array size, on stack */
+
+/* Extend the partials array p[] by doubling its size. */
+static int                          /* non-zero on error */
+_sum_realloc(double **p_ptr, Py_ssize_t  n,
+             double  *ps,    Py_ssize_t *m_ptr)
+{
+	void *v = NULL;
+	Py_ssize_t m = *m_ptr;
+
+	m += m;  /* double */
+	if (n < m && m < (PY_SSIZE_T_MAX / sizeof(double))) {
+		double *p = *p_ptr;
+		if (p == ps) {
+			v = PyMem_Malloc(sizeof(double) * m);
+			if (v != NULL)
+				memcpy(v, ps, sizeof(double) * n);
+		}
+		else
+			v = PyMem_Realloc(p, sizeof(double) * m);
+	}
+	if (v == NULL) {        /* size overflow or no memory */
+		PyErr_SetString(PyExc_MemoryError, "math sum partials");
+		return 1;
+	}
+	*p_ptr = (double*) v;
+	*m_ptr = m;
+	return 0;
+}
+
+/* Full precision summation of a sequence of floats.
+
+   def msum(iterable):
+       partials = []  # sorted, non-overlapping partial sums
+       for x in iterable:
+           i = 0
+           for y in partials:
+               if abs(x) < abs(y):
+                   x, y = y, x
+               hi = x + y
+               lo = y - (hi - x)
+               if lo:
+                   partials[i] = lo
+                   i += 1
+               x = hi
+           partials[i:] = [x]
+       return sum_exact(partials)
+
+   Rounded x+y stored in hi with the roundoff stored in lo.  Together hi+lo
+   are exactly equal to x+y.  The inner loop applies hi/lo summation to each
+   partial so that the list of partial sums remains exact.
+
+   Sum_exact() adds the partial sums exactly and correctly rounds the final
+   result (using the round-half-to-even rule).  The items in partials remain
+   non-zero, non-special, non-overlapping and strictly increasing in
+   magnitude, but possibly not all having the same sign.
+
+   Depends on IEEE 754 arithmetic guarantees and half-even rounding.
+*/
+
+static PyObject*
+math_sum(PyObject *self, PyObject *seq)
+{
+	PyObject *item, *iter, *sum = NULL;
+	Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
+	double x, y, hi, lo=0.0, ps[NUM_PARTIALS], *p = ps;
+
+	iter = PyObject_GetIter(seq);
+	if (iter == NULL)
+		return NULL;
+
+	PyFPE_START_PROTECT("sum", Py_DECREF(iter); return NULL)
+
+	for(;;) {           /* for x in iterable */
+		assert(0 <= n && n <= m);
+		assert((m == NUM_PARTIALS && p == ps) ||
+		       (m >  NUM_PARTIALS && p != NULL));
+
+		item = PyIter_Next(iter);
+		if (item == NULL) {
+			if (PyErr_Occurred())
+				goto _sum_error;
+			break;
+		}
+		x = PyFloat_AsDouble(item);
+		Py_DECREF(item);
+		if (PyErr_Occurred())
+			goto _sum_error;
+
+		for (i = j = 0; j < n; j++) {       /* for y in partials */
+			y = p[j];
+			hi = x + y;
+			lo = fabs(x) < fabs(y)
+			   ? x - (hi - y)
+			   : y - (hi - x);
+			if (lo != 0.0)
+				p[i++] = lo;
+			x = hi;
+		}
+		
+		n = i;                              /* ps[i:] = [x] */                   
+		if (x != 0.0) {
+			/* If non-finite, reset partials, effectively
+			   adding subsequent items without roundoff
+			   and yielding correct non-finite results,
+			   provided IEEE 754 rules are observed */
+			if (! Py_IS_FINITE(x))
+				n = 0;
+			else if (n >= m && _sum_realloc(&p, n, ps, &m))
+				goto _sum_error;
+			p[n++] = x;
+		}
+	}
+
+	if (n > 0) {
+		hi = p[--n];
+		if (Py_IS_FINITE(hi)) {
+			/* sum_exact(ps, hi) from the top, stop when the sum becomes inexact. */
+			while (n > 0) {
+				x = p[--n];
+				y = hi;
+				hi = x + y;
+				assert(fabs(x) < fabs(y));
+				lo = x - (hi - y);
+				if (lo != 0.0)
+					break;
+			}
+			/* Little dance to allow half-even rounding across multiple partials.
+                           Needed so that sum([1e-16, 1, 1e16]) will round-up to two instead
+                           of down to zero (the 1e16 makes the 1 slightly closer to two). */
+			if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) ||
+			              (lo > 0.0 && p[n-1] > 0.0))) {
+				y = lo * 2.0;
+				x = hi + y;
+				if (y == (x - hi))
+					hi = x;
+			}
+		}
+		else {  /* raise corresponding error */
+			errno = Py_IS_NAN(hi) ? EDOM : ERANGE;
+			if (is_error(hi))
+				goto _sum_error;
+		}
+	}
+	else  /* default */
+		hi = 0.0;
+	sum = PyFloat_FromDouble(hi);
+
+_sum_error:
+	PyFPE_END_PROTECT(hi)
+	Py_DECREF(iter);
+	if (p != ps)
+		PyMem_Free(p);
+	return sum;
+}
+
+#undef NUM_PARTIALS
+
+PyDoc_STRVAR(math_sum_doc,
+"sum(iterable)\n\n\
+Return an accurate floating point sum of values in the iterable.\n\
+Assumes IEEE-754 floating point arithmetic.");
+
 static PyObject *
 math_trunc(PyObject *self, PyObject *number)
 {
@@ -833,6 +1026,7 @@
 	{"sin",		math_sin,	METH_O,		math_sin_doc},
 	{"sinh",	math_sinh,	METH_O,		math_sinh_doc},
 	{"sqrt",	math_sqrt,	METH_O,		math_sqrt_doc},
+	{"sum",		math_sum,	METH_O,		math_sum_doc},
 	{"tan",		math_tan,	METH_O,		math_tan_doc},
 	{"tanh",	math_tanh,	METH_O,		math_tanh_doc},
  	{"trunc",	math_trunc,	METH_O,		math_trunc_doc},