--- Python-2.6a3/Modules/mathmoduleORIG.c 2008-04-20 18:55:50.000000000 -0700 +++ Python-2.6a3/Modules/mathmodule.c 2008-05-18 14:07:34.000000000 -0700 @@ -55,6 +55,8 @@ #include "Python.h" #include "longintrepr.h" /* just for SHIFT */ +#include /* for FLT_RADIX and DBL_MAX_EXP */ + #ifdef _OSF_SOURCE /* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */ extern double copysign(double, double); @@ -307,6 +309,297 @@ FUNC1(tanh, tanh, 0, "tanh(x)\n\nReturn the hyperbolic tangent of x.") +/* Summation function enhanced to handle intermediate overflow and + rounding errors like function msum() version written by Mark + Dickinson posted at . + + The orginal msum() function by Raymond Hettinger is this recipe + . + + See both of those for more details, proofs and other references. + + Assumes IEEE 754 floating point format and semantics, specifically + the following rules: + + 1. if the summands include a NaN, return a NaN, + + 2. if the summands include infinities of both signs, raise ValueError, + + 3. if the summands include infinities of only one sign, return + infinity with that sign, + + 4. otherwise (all summands are finite) if the result is infinite, + raise OverflowError. The result can never be a NaN if all + summands are finite. + + The same summation function ais also in ./cmathmodule.c, except for + 2 lines marked *** below. Make sure to updated both when changing. +*/ + + +#define DO_SUM_PS 128 /* stack partials */ +#define DO_SUM_SIGNOF(x) ((x) < 0 ? 0 : 1) + +static double _do_sum_pow_2 = 0.0; /* FLT_RADIX ** DBL_MAX_EXP / 2 */ + +static double +_do_sum_add2(double x, double y, double *lo_ptr) +{ + double hi; + + errno = 0; + hi = x + y; + /* follow the IEEE 754 rules for summation */ + if (Py_IS_NAN(hi)) { + if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) + errno = EDOM; + else + errno = 0; /* rule 1 */ + } + else if (Py_IS_INFINITY(hi)) { + /* *** see is_error() for hi != 0.0 */ + if (hi != 0.0 && Py_IS_FINITE(x) && Py_IS_FINITE(y)) + errno = ERANGE; /* rule 4 */ + else if (x != y && Py_IS_INFINITY(x) && Py_IS_INFINITY(y)) + errno = EDOM; /* rule 2 */ + else + errno = 0; /* rule 3 */ + } + if (lo_ptr) + *lo_ptr = fabs(x) < fabs(y) + ? x - (hi - y) + : y - (hi - x); + return hi; +} + +static double +_do_sum_adjust(double *x_ptr) +{ + double x; + + /* get sys.float_info.radix ** sys.float_info.max_exp / 2.0 */ + if (_do_sum_pow_2 == 0.0) +#if FLT_RADIX == 2 + _do_sum_pow_2 = pow((double)FLT_RADIX, (double)(DBL_MAX_EXP - 1)); +#else + _do_sum_pow_2 = pow((double)FLT_RADIX, (double)DBL_MAX_EXP) / 2.0; +#endif + + if ((x = *x_ptr) > 0.0) { + *x_ptr = x - _do_sum_pow_2 - _do_sum_pow_2; + return 1.0; + } else { + *x_ptr = x + _do_sum_pow_2 + _do_sum_pow_2; + return -1.0; + } +} + +static double +_do_sum_exact(double top, double *p, Py_ssize_t *n_ptr) +{ + double lo, hi = top; + Py_ssize_t n = *n_ptr; + + /* sum from the top (most significant), + stop as soon as sum becomes inexact */ + while (n-- > 0) { + hi = _do_sum_add2(hi, p[n], &lo); + if (lo != 0.0 || errno) { + p[n++] = lo; + *n_ptr = n; + break; + } + } + /* round correctly if necessary */ + if (n-- > 1) { + lo = p[n]; + if (DO_SUM_SIGNOF(lo) == DO_SUM_SIGNOF(p[n-1])) { + double t = lo * 2.0; + if (t == (_do_sum_add2(hi, t, NULL) - hi)) { + hi += t; + p[n] = -lo; + } + } + } + return hi; +} + +static int +_do_sum(PyObject *seq, /* iterable sequence */ + PyObject *sum, /* optional start value (0) */ + double (*func)(PyObject *), /* callback returning item as double */ + double *sum_ptr) /* sum result */ +{ + int err = 1; /* meaning error */ + Py_ssize_t i, j, n = 0, m = DO_SUM_PS; + double x, y, t, hi, lo, ps[DO_SUM_PS], *p = ps, r = 0.0; /* result */ + PyObject *iter, *item; + + iter = PyObject_GetIter(seq); + if (iter == NULL) + return err; + + if (sum != NULL) { + t = func(sum); /* e.g. PyFloat_AsDouble() */ + if (PyErr_Occurred()) + goto _do_sum_exit; + p[n++] = t; + } else + p[0] = 0.0; + + PyFPE_START_PROTECT("sum", goto double_sum_exit) + + /* for x in iterable */ + for(;;) { + /* some invariants */ + assert(err != 0); + assert(n >= 0); + assert(m >= n); + assert((m == DO_SUM_PS && p == ps) || + (m > DO_SUM_PS && p != NULL)); + + item = PyIter_Next(iter); + if (item == NULL) { + if (PyErr_Occurred()) + goto _do_sum_error; + else + break; + } + x = func(item); /* e.g. PyFloat_AsDouble() */ + Py_DECREF(item); + if (PyErr_Occurred()) + goto _do_sum_error; + if (Py_IS_INFINITY(x)) + r += x; + else if (Py_IS_NAN(x)) { + r = x; + goto _do_sum_ok; + } else if (x != 0.0) { + /* i = 0; for y in partials */ + for (i = j = 0; j < n; j++) { + y = p[j]; + hi = _do_sum_add2(x, y, &lo); + if (/* errno == ERANGE && */ Py_IS_INFINITY(hi)) { + /* intermediate overflow, retry + after adjusting largest */ + r += fabs(x) < fabs(y) + ? _do_sum_adjust(&y) + : _do_sum_adjust(&x); + hi = _do_sum_add2(x, y, &lo); + } + if (errno) { + r = hi; + goto _do_sum_errno; + } + /* ignore non-finite lo */ + if (lo != 0.0 && Py_IS_FINITE(lo)) + p[i++] = lo; + x = hi; + } + /* ps[i:] = [x] */ + n = i; + if (n >= m) { /* extend partials */ + void *v = NULL; + m += m / 2; /* plus one half */ + 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 _do_sum_error; + } + p = (double*) v; + } + if (x != 0.0) + p[n++] = x; + } + } + + if ((t = r) != 0.0) { + /* adjust for overflow if valid */ + if (n > 0) + t = p[--n]; + if (t < 0.0 && r == 1.0) + x = _do_sum_pow_2; + else if (t >= 0.0 && r == -1.0) + x = -_do_sum_pow_2; + else if (Py_IS_INFINITY(r)) + goto _do_sum_ok; + else { /* rule 2 or overflow */ + errno = Py_IS_NAN(r) ? EDOM : ERANGE; + goto _do_sum_errno; + } + y = t / 2.0; + hi = _do_sum_add2(x, y, &lo); + r = hi * 2.0; + if (Py_IS_INFINITY(r)) { + /* overflow except corner cases */ + if (n > 0 && DO_SUM_SIGNOF(lo) == DO_SUM_SIGNOF(p[n-1])) { + x = lo * 2.0; + y = _do_sum_add2(hi, x, NULL); + if (x == (y - hi)) { /* corner case */ + r = y * 2.0; + goto _do_sum_ok; + } + } + /* overflow */ + errno = ERANGE; + goto _do_sum_errno; + } else { /* replace top */ + if (lo != 0.0) + p[n++] = lo * 2.0; + t = hi * 2.0; + } + } + /* sum(ps, t) */ + r = _do_sum_exact(t, p, &n); +_do_sum_ok: + errno = err = 0; + +_do_sum_errno: + *sum_ptr = r; + if (err && errno) + is_error(1.0); /* *** any non-zero */ + +_do_sum_error: + PyFPE_END_PROTECT(r) + +_do_sum_exit: + Py_DECREF(iter); + if (p != ps) + PyMem_Free(p); + return err; +} + +#undef DO_SUM_PS +#undef DO_SUM_SIGNOF + +static PyObject* +math_sum(PyObject *self, PyObject *args) +{ + double s = 0.0; + PyObject *seq, *sum = NULL; + + if (!PyArg_UnpackTuple(args, "sum", 1, 2, &seq, &sum)) + return NULL; + + if (_do_sum(seq, sum, PyFloat_AsDouble, &s)) + return NULL; + + return PyFloat_FromDouble(s); +} + +PyDoc_STRVAR(math_sum_doc, +"sum(sequence[, start=0[, limit=4e3])\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.\n"); + static PyObject * math_trunc(PyObject *self, PyObject *number) { @@ -718,6 +1011,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},