--- Python-2.6a3/Modules/mathmoduleORIG.c	2008-04-20 18:55:50.000000000 -0700
+++ Python-2.6a3/Modules/mathmodule.c	2008-05-12 08:51:15.000000000 -0700
@@ -307,6 +307,186 @@
 FUNC1(tanh, tanh, 0,
       "tanh(x)\n\nReturn the hyperbolic tangent of x.")
 
+/* Summation function as function msum() in this Python recipe
+   <http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>
+
+   def msum(iterable):
+       "Full precision summation using multiple floats for intermediate values"
+
+       # Rounded x+y stored in hi with the round-off 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.
+
+       # Depends on IEEE-754 arithmetic guarantees.  See proof of correctness at:
+       # www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps
+
+       ps = []  # sorted, non-overlapping partial sums
+       for x in iterable:
+           i = 0
+           for y in ps:
+               if abs(x) < abs(y):
+                   x, y = y, x
+               hi = x + y
+               lo = y - (hi - x)
+               if lo:
+                   ps[i] = lo
+                   i += 1
+               x = hi
+           ps[i:] = [x]
+       return sum(ps, 0.0)
+
+   Note, a duplicate of this code is in file ./cmathmodule.c.
+   Make any changes in both places.
+*/
+
+static double
+_math_sum_add(double x, double y)
+{
+	double s;
+
+	errno = 0;
+	s = x + y;
+	/* following the IEEE 754 rules for summation:
+
+	1. if the summands include a NaN, return a NaN
+
+	2. else if the summands include infinities of
+	   both signs, raise ValueError,
+
+	3. else if the summands include infinities of
+	   only one sign, return infinity with that sign,
+
+	4. else (all summands are finite) if the result
+	   is infinite, raise OverflowError.  The result
+	   can never be a NaN if all summands are finite.
+	*/
+	if (Py_IS_NAN(s)) {
+		if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
+			errno = EDOM;
+		else
+			errno = 0;  /* rule 1 */
+	}
+	else if (Py_IS_INFINITY(s)) {
+		if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
+			errno = ERANGE;  /* rule 4 */
+		else if (Py_IS_INFINITY(x) && Py_IS_INFINITY(y) && x != y)
+			errno = EDOM;  /* rule 2 */
+		else
+			errno = 0;  /* rule 3 */
+	}
+	return s;
+}
+
+static PyObject*
+math_sum(PyObject *self, PyObject *args)
+{
+#	define MATH_SUM_PS  128  /* stack partials */
+
+	Py_ssize_t i, j, n = 0, m = MATH_SUM_PS;
+	double x, y, s, hi, lo, ps[MATH_SUM_PS], *p = ps;
+	PyObject *seq, *iter, *item, *sum = NULL;
+
+	if (!PyArg_UnpackTuple(args, "sum", 1, 2, &seq, &sum))
+		return NULL;
+
+	iter = PyObject_GetIter(seq);
+	if (iter == NULL)
+		return NULL;
+
+	if (sum != NULL) { 
+		s = PyFloat_AsDouble(sum);
+		sum = NULL;
+		if (s == -1.0 && PyErr_Occurred())
+			goto math_sum_exit;
+		p[n++] = s;
+	} else
+		p[0] = 0.0;
+
+	PyFPE_START_PROTECT("sum", goto math_sum_exit)
+
+	/* for x in iterable */
+	for(;;) {
+		/* some invariants */
+		assert(sum == NULL);
+		assert(n >= 0);
+		assert(m >= n);
+		assert((m == MATH_SUM_PS && p == ps) ||
+		       (m >  MATH_SUM_PS && p != NULL));
+
+		item = PyIter_Next(iter);
+		if (item == NULL) {
+			if (PyErr_Occurred())
+				goto math_sum_error;
+			else
+				break;
+		}
+		x = PyFloat_AsDouble(item);
+		Py_DECREF(item);
+		if (x == -1.0 && PyErr_Occurred())
+			goto math_sum_error;
+		/* i = 0; for y in partials */
+		for (i = j = 0; j < n; j++) {
+			y = p[j];
+			hi = _math_sum_add(x, y);
+			if (errno && is_error(hi))
+				goto math_sum_error;
+			/* lo = min_x_y - (hi - max_x_y) */
+			lo = fabs(x) < fabs(y)
+			   ? x - (hi - y)
+			   : y - (hi - x);
+			if (lo != 0.0)
+				p[i++] = lo;
+			x = hi;
+		}
+		/* ps[i:] = [x] */
+		n = i;
+		if (n >= m) {  /* extend partials */
+			void *v = NULL;
+			m += m / 2;  /* plus 50% */
+			if (n < m && m < (PY_SSIZE_T_MAX / sizeof(double))) {
+				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) {  /* no memory or size overflow */
+				PyErr_SetString(PyExc_MemoryError, "math sum error");
+				goto math_sum_error;
+			}
+			p = (double*) v;
+		}
+		p[n++] = x;
+	}
+
+	/* sum(ps, 0.0) */
+	s = p[0];  /* always valid */
+	for (i = 1; i < n; i++) {
+		s = _math_sum_add(p[i], s);
+		if (errno && is_error(s))
+			goto math_sum_error;
+	}
+	sum = PyFloat_FromDouble(s);
+
+math_sum_error:
+	PyFPE_END_PROTECT(s)
+
+math_sum_exit:
+	Py_DECREF(iter);
+	if (p != ps)
+		PyMem_Free(p);
+	return sum;
+
+#	undef MATH_SUM_PS
+}
+
+PyDoc_STRVAR(math_sum_doc,
+"sum(sequence[, start=0])\n"
+"\n"
+"Return the full precision sum of a sequence of numbers plus the value\n"
+"of parameter 'start'.  When the sequence is empty, return start.");
+
 static PyObject *
 math_trunc(PyObject *self, PyObject *number)
 {
@@ -718,6 +898,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_VARARGS,	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},