--- Python-2.6a3/Modules/mathmoduleORIG.c 2008-04-20 18:55:50.000000000 -0700 +++ Python-2.6a3/Modules/mathmodule.c 2008-05-19 15:04:21.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,347 @@ FUNC1(tanh, tanh, 0, "tanh(x)\n\nReturn the hyperbolic tangent of x.") + +#if FLT_RADIX == 2 /* function sum not implemented otherwise */ + +/* Precision summation function handling intermediate overflow and + rounding errors like function msum() version written by Mark + Dickinson posted at . + + The original 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 is also in ./cmathmodule.c, except for + the 2 lines marked *** below. Make sure to update both. +*/ + +#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 */ + +/* Add two floats returning the sum and optionally the + round-off. Global errno is set to ERANGE on overflow, + to EDOM for a NaN result or to zero otherwise. + */ +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)) { + if (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 */ + } + /* *** see is_error() above */ + if (errno == ERANGE && hi == 0.0) + errno = 0; + if (lo_ptr) + *lo_ptr = fabs(x) < fabs(y) + ? x - (hi - y) + : y - (hi - x); + return hi; +} + +/* Adjust largest summand of an intermediate overflow. + */ +static double +_do_sum_adjust(double *x_ptr) +{ + double x = *x_ptr; + + if (_do_sum_pow_2 == 0.0) { + assert(FLT_RADIX == 2); + /* sys.float_info.radix ** sys.float_info.max_exp / 2.0 */ + _do_sum_pow_2 = ldexp(1.0, DBL_MAX_EXP - 1); + } + if (x > 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; + } +} + +/* Compute the sum of array p[] of nonoverlapping floats. + + On input, p[] is an array of nonzero, nonspecial, non- + overlapping floats, strictly increasing in magnitude, but + possibly not all having the same sign. + + On output, the sum of p[] gives the error in the returned + result, which is correctly rounded (using the round-half- + to-even rule). The elements of p[] remain nonzero, non- + special, nonoverlapping, and increasing in magnitude. + + Assumes IEEE 754 floating point format and semantics. + */ +static double +_do_sum_exact(double top, double *p, Py_ssize_t *n_ptr) +{ + double pn, 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) { + pn = p[--n]; + hi = _do_sum_add2(hi, pn, &lo); + if (lo != 0.0 || errno) { + p[n++] = lo; + *n_ptr = n; + break; + } + } + /* round correctly if necessary */ + if (n > 1) { + pn = p[--n]; + if (DO_SUM_SIGNOF(pn) == DO_SUM_SIGNOF(p[n-1])) { + double t = pn * 2.0; + double x = _do_sum_add2(hi, t, NULL); + if (t == (x - hi) && errno == 0) { + hi = x; + p[n] = -pn; + } + } + } + return hi; +} + +/* Extend the partials array p[] */ +static int +_do_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 / 2; /* extend by one half */ + 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 error"); + return 1; + } + *p_ptr = (double*) v; + *m_ptr = m; + return 0; +} + +/* Full precision summation of an sequence of floats + handling intermediate overflow and rounding errors. + + Callback function func is PyFloat_AsDouble() or any + other function returning the item value as a double. + */ +static int /* non-zero on error */ +_do_sum(PyObject *seq, /* iterable sequence */ + PyObject *sum, /* optional start value (0) */ + double (*func)(PyObject *), /* callback */ + double *sum_ptr) /* result, not set on error */ +{ + 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 _do_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; /* rule 1 */ + 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) + 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 && _do_sum_realloc(&p, n, ps, &m)) + goto _do_sum_error; + if (x != 0.0) + p[n++] = x; + } + } + + assert(err != 0); + if ((t = r) != 0.0) { + 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)) /* rule 3 */ + goto _do_sum_ok; + else { /* rule 2 or 4 */ + errno = Py_IS_NAN(r) ? EDOM : ERANGE; + goto _do_sum_errno; + } + /* adjust for intermediate overflow */ + hi = _do_sum_add2(x, t / 2.0, &lo); + if (errno == 0) { + t = _do_sum_add2(hi, hi, NULL); + if (errno == ERANGE && Py_IS_INFINITY(t)) { + /* overflow except corner cases */ + if (n > 0 && DO_SUM_SIGNOF(lo) == + DO_SUM_SIGNOF(p[n-1])) { + t = lo * 2.0; + x = _do_sum_add2(hi, t, NULL); + if (t == (x - hi) && errno == 0) { + r = x * 2.0; /* corner case */ + goto _do_sum_ok; + } + } + /* overflow */ + errno = ERANGE; + goto _do_sum_errno; + } + /* new top t */ + if (lo != 0.0) + p[n++] = lo * 2.0; + } + } + /* sum(ps, t) */ + r = _do_sum_exact(t, p, &n); +_do_sum_ok: + *sum_ptr = r; + err = 0; + +_do_sum_errno: + 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 +#endif /* FLT_RADIX */ + +static PyObject* +math_sum(PyObject *self, PyObject *args) +{ +#if FLT_RADIX == 2 + 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); + +#else /* FLT_RADIX */ + PyErr_SetString(PyExc_NotImplementedError, + "math sum radix"); + return NULL; +#endif /* FLT_RADIX */ +} + +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.\n\n\ +Note, this function assumes IEEE 754 floating point format and semantics.\n\ +If the floating point radix is not 2, a NotImplementedError is raised."); + static PyObject * math_trunc(PyObject *self, PyObject *number) { @@ -718,6 +1061,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},