mirror of
https://github.com/python/cpython
synced 2026-09-29 12:10:30 +03:00
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. ........
This commit is contained in:
Benjamin Peterson
parent
1b466f25e1
commit
2b7411df5c
7 changed files with 368 additions and 4 deletions
|
|
@ -362,6 +362,199 @@ FUNC1(tan, tan, 0,
|
|||
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 @@ static PyMethodDef math_methods[] = {
|
|||
{"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},
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue