[issue2819] Full precision summation

[Mark Dickinson]

Mark Dickinson <dickinsm at gmail.com> added the comment:

Here's (msum.py) an example in Python of one fairly straightforward way of 
dealing with overflow correctly, without needing more than one pass 
through the data, and without significant slowdown in the normal case.
(The Python code is needlessly inefficient in places, notably in that 
partials[1:] creates a new list;  this is obviously not a problem in C.)

The idea is essentially just to maintain the sum modulo integer multiples 
of 2.**1024.  partials[0] is reserved for keeping track of the current 
multiple of 2.**!024.  So at each stage, the sum so far is 
sum(partials[1:], 0.0) + 2.**1024 * partials[0].

I'm 97.3% convinced that the proof of correctness goes through:  it's 
still true with this modification that partials always consists of 
nonadjacent, nonzero values of increasing magnitude.  One of the keys to proving this is to note that for any value x between 2**1023 and 2**1024, 
both x-2**1023 and x-2**1024 are exactly representable.

---

There's one more 'nice-to-have' that I think should be considered:  it 
would be nice if the result of msum were always correctly rounded.  One 
aspect of correct rounding is that it provides a guarantee that the sum is 
independent of the order of the summands, so

msum(list) == msum(sorted(list))

would hold true.

Added file: http://bugs.python.org/file10345/msum.py

__________________________________
Tracker <report at bugs.python.org>
<http://bugs.python.org/issue2819>
__________________________________