Merged revisions 63542-63544,63546,63553,63563-63564,63567,63569,63576 via svnmerge from
svn+ssh://pythondev@svn.python.org/python/trunk
........
r63542 | mark.dickinson | 2008-05-22 20:35:30 -0500 (Thu, 22 May 2008) | 5 lines
Issue #2819: Add math.sum, a function that sums a sequence of floats
efficiently but with no intermediate loss of precision. Based on
Raymond Hettinger's ASPN recipe. Thanks Jean Brouwers for the patch.
........
r63543 | mark.dickinson | 2008-05-22 21:36:48 -0500 (Thu, 22 May 2008) | 2 lines
Add tests for math.sum (Issue #2819)
........
r63544 | mark.dickinson | 2008-05-22 22:30:01 -0500 (Thu, 22 May 2008) | 2 lines
Better error reporting in test_math.py
........
r63546 | raymond.hettinger | 2008-05-22 23:32:43 -0500 (Thu, 22 May 2008) | 1 line
Tweak the comments and formatting.
........
r63553 | mark.dickinson | 2008-05-23 07:07:36 -0500 (Fri, 23 May 2008) | 3 lines
Skip math.sum tests on non IEEE 754 platforms, and on IEEE 754 platforms
that exhibit the problem described in issue #2937.
........
r63563 | martin.v.loewis | 2008-05-23 10:18:28 -0500 (Fri, 23 May 2008) | 3 lines
Issue #1390: Raise ValueError in toxml when an invalid comment would
otherwise be produced.
........
r63564 | raymond.hettinger | 2008-05-23 12:21:44 -0500 (Fri, 23 May 2008) | 1 line
Issue 2909: show how to name unpacked fields.
........
r63567 | raymond.hettinger | 2008-05-23 12:34:34 -0500 (Fri, 23 May 2008) | 1 line
Fix typo
........
r63569 | martin.v.loewis | 2008-05-23 14:33:13 -0500 (Fri, 23 May 2008) | 3 lines
Mention that the leaking of variables from list comprehensions
is fixed in 3.0.
........
r63576 | martin.v.loewis | 2008-05-24 04:36:45 -0500 (Sat, 24 May 2008) | 3 lines
Don't try to get the window size if it was never set before.
Fixes the test failure on Solaris.
........
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},