Jean Brouwers <[EMAIL PROTECTED]> added the comment:
Mark,
Thank you very much for your comments. Here is my initial response to the
first 3.
(1) Attached is an attempt to address the 1st issue (just the
mathmodule). The macros PyFPE_START_PROTECT/_END_PROTECT have been moved
outside the main loop and the errno is set following the IEEE 754 rules as
you suggested.
One related issue is testing these, how can a NaN and +/-Infinity float
object be created in Python?
(2) On your 2nd comment, supporting non-IEEE floating point, perhaps the
Kahan method should be used in that case. If so, the next question is how
to detect that? There are two symbols in pyconfig.h HAVE_IEEEFP_H and
HAVE_LIBIEEE. Are those the proper ones to determine IEEE floating point
support?
(3) On the 3rd comment, Raymond and I did discus using a single function to
be called by the math and cmath modules. The question is where should that
function reside? The math and cmath modules are not the right place since
both are loadable modules. It will have to be somewhere inside the Python
main/core.
Also, depending on the implementation of that function, it may require
iterating the complex sequence twice. And that will force the C complex
numbers to be created twice by the PyComplex_AsCComplex() call. Would that
be a concern?
/Jean Brouwers
On Sun, May 11, 2008 at 1:49 PM, Mark Dickinson <[EMAIL PROTECTED]>
wrote:
>
> Mark Dickinson <[EMAIL PROTECTED]> added the comment:
>
> One more question:
>
> What are the use cases for an exact summation algorithm? That is, in what
> situations does one care about exactness rather than simply accuracy? I
> know that loss of accuracy is a problem in things like numeric integration
> routines, but something like Kahan summation (faster and simpler, but not
> exact) usually takes care of that.
>
> __________________________________
> Tracker <[EMAIL PROTECTED]>
> <http://bugs.python.org/issue2819>
> __________________________________
>
Added file: http://bugs.python.org/file10306/unnamed
Added file: http://bugs.python.org/file10307/mathmodule.c.2.6a3.take2.diff
__________________________________
Tracker <[EMAIL PROTECTED]>
<http://bugs.python.org/issue2819>
__________________________________
Mark,<div><br></div><div>Thank you very much for your comments. Here is
my initial response to the first 3.</div><div><br></div><div><br></div><div>(1)
Attached is an attempt to address the 1st issue (just the mathmodule). The
macros <span class="Apple-style-span" style="border-collapse: collapse;
">PyFPE_START_PROTECT/_END_PROTECT<span class="Apple-style-span"
style="border-collapse: separate; "> have been moved outside the main loop
and the errno is set following the IEEE 754 rules as you
suggested.</span></span></div>
<div><br></div><div>One related issue is testing these, how can a NaN and
+/-Infinity float object be created in
Python?</div><div><br></div><div><br></div><div><div>(2) On your 2nd comment,
supporting non-IEEE floating point, perhaps the Kahan method should be used in
that case. If so, the next question is how to detect that? There
are two symbols in pyconfig.h HAVE_IEEEFP_H and HAVE_LIBIEEE. Are those
the proper ones to determine IEEE floating point support?</div>
<div><br></div><div><br></div><div>(3) On the 3rd comment, Raymond and I did
discus using a single function to be called by the math and cmath modules.
The question is where should that function reside? The math and
cmath modules are not the right place since both are loadable modules. It
will have to be somewhere inside the Python main/core.</div>
<div><br></div><div>Also, depending on the implementation of that function, it
may require iterating the complex sequence twice. And that will force the
C complex numbers to be created twice by the PyComplex_AsCComplex() call.
Would that be a concern?</div>
<div><br></div><div><br></div><div>/Jean
Brouwers</div><div><br></div><div> </div><div><br></div><div
class="gmail_quote">On Sun, May 11, 2008 at 1:49 PM, Mark Dickinson <<a
href="mailto:[EMAIL PROTECTED]">[EMAIL PROTECTED]</a>> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc
solid;padding-left:1ex;"><div class="Ih2E3d"><br>
Mark Dickinson <<a href="mailto:[EMAIL PROTECTED]">[EMAIL PROTECTED]</a>>
added the comment:<br>
<br>
</div>One more question:<br>
<br>
What are the use cases for an exact summation algorithm? That is, in
what<br>
situations does one care about exactness rather than simply accuracy?
I<br>
know that loss of accuracy is a problem in things like numeric integration<br>
routines, but something like Kahan summation (faster and simpler, but not<br>
exact) usually takes care of that.<br>
<div><div></div><div class="Wj3C7c"><br>
__________________________________<br>
Tracker <<a href="mailto:[EMAIL PROTECTED]">[EMAIL PROTECTED]</a>><br>
<<a href="http://bugs.python.org/issue2819";
target="_blank">http://bugs.python.org/issue2819</a>><br>
__________________________________<br>
</div></div></blockquote></div><br></div>
--- Python-2.6a3/Modules/mathmoduleORIG.c 2008-04-20 18:55:50.000000000
-0700
+++ Python-2.6a3/Modules/mathmodule.c 2008-05-12 08:51:15.000000000 -0700
@@ -307,6 +307,186 @@
FUNC1(tanh, tanh, 0,
"tanh(x)\n\nReturn the hyperbolic tangent of x.")
+/* Summation function as function msum() in this Python recipe
+ <http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>
+
+ 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 double
+_math_sum_add(double x, double y)
+{
+ double s;
+
+ errno = 0;
+ s = x + y;
+ /* following the IEEE 754 rules for summation:
+
+ 1. if the summands include a NaN, return a NaN
+
+ 2. else if the summands include infinities of
+ both signs, raise ValueError,
+
+ 3. else if the summands include infinities of
+ only one sign, return infinity with that sign,
+
+ 4. else (all summands are finite) if the result
+ is infinite, raise OverflowError. The result
+ can never be a NaN if all summands are finite.
+ */
+ if (Py_IS_NAN(s)) {
+ if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
+ errno = EDOM;
+ else
+ errno = 0; /* rule 1 */
+ }
+ else if (Py_IS_INFINITY(s)) {
+ if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
+ errno = ERANGE; /* rule 4 */
+ else if (Py_IS_INFINITY(x) && Py_IS_INFINITY(y) && x != y)
+ errno = EDOM; /* rule 2 */
+ else
+ errno = 0; /* rule 3 */
+ }
+ return s;
+}
+
+static PyObject*
+math_sum(PyObject *self, PyObject *args)
+{
+# define MATH_SUM_PS 128 /* stack 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;
+
+ PyFPE_START_PROTECT("sum", goto math_sum_exit)
+
+ /* for x in iterable */
+ for(;;) {
+ /* some invariants */
+ assert(sum == NULL);
+ assert(n >= 0);
+ assert(m >= n);
+ assert((m == MATH_SUM_PS && p == ps) ||
+ (m > MATH_SUM_PS && p != NULL));
+
+ item = PyIter_Next(iter);
+ if (item == NULL) {
+ if (PyErr_Occurred())
+ goto math_sum_error;
+ else
+ break;
+ }
+ x = PyFloat_AsDouble(item);
+ Py_DECREF(item);
+ if (x == -1.0 && PyErr_Occurred())
+ goto math_sum_error;
+ /* i = 0; for y in partials */
+ for (i = j = 0; j < n; j++) {
+ y = p[j];
+ hi = _math_sum_add(x, y);
+ if (errno && is_error(hi))
+ goto math_sum_error;
+ /* lo = min_x_y - (hi - max_x_y) */
+ 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_error;
+ }
+ p = (double*) v;
+ }
+ p[n++] = x;
+ }
+
+ /* sum(ps, 0.0) */
+ s = p[0]; /* always valid */
+ for (i = 1; i < n; i++) {
+ s = _math_sum_add(p[i], s);
+ if (errno && is_error(s))
+ goto math_sum_error;
+ }
+ sum = PyFloat_FromDouble(s);
+
+math_sum_error:
+ PyFPE_END_PROTECT(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 +898,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},
_______________________________________________
Python-bugs-list mailing list
Unsubscribe:
http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
