--- Python-2.6a3/Modules/mathmoduleORIG.c 2008-04-20 18:55:50.000000000 -0700 +++ Python-2.6a3/Modules/mathmodule.c 2008-05-10 16:34:03.000000000 -0700 @@ -307,6 +307,162 @@ FUNC1(tanh, tanh, 0, "tanh(x)\n\nReturn the hyperbolic tangent of x.") +/* Summation function as function msum() in this Python recipe + + + 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 int +_math_sum_add(double x, double y, double *s_ptr) +{ + double s; + + errno = 0; + PyFPE_START_PROTECT("sum", return -1) + s = x + y; + PyFPE_END_PROTECT(s) + /* like math_2() above */ + if (Py_IS_NAN(s)) { + if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) + errno = EDOM; + else + errno = 0; + } + else if (Py_IS_INFINITY(s)) { + if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) + errno = ERANGE; + else + errno = 0; + } + if (errno != 0 && is_error(s)) + return -1; + + *s_ptr = s; + return 0; +} + +static PyObject* +math_sum(PyObject *self, PyObject *args) +{ +# define MATH_SUM_PS 128 /* static 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; + + /* for x in iterable */ + for(;;) { + item = PyIter_Next(iter); + if (item == NULL) { + if (PyErr_Occurred()) + goto math_sum_exit; + else + break; + } + x = PyFloat_AsDouble(item); + Py_DECREF(item); + if (x == -1.0 && PyErr_Occurred()) + goto math_sum_exit; + /* i = 0; for y in partials */ + for (i = j = 0; j < n; j++) { + y = p[j]; + if (_math_sum_add(x, y, &hi) < 0) + goto math_sum_exit; + 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_exit; + } + p = (double*) v; + } + p[n++] = x; + } + + /* sum(ps, 0.0) */ + s = p[0]; /* always valid */ + for (i = 1; i < n; i++) { + if (_math_sum_add(p[i], s, &s) < 0) + goto math_sum_exit; + } + sum = PyFloat_FromDouble(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 +874,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},