Full precision summation

[MrJean1]

BPO	2819
Nosy	@tim-one, @rhettinger, @terryjreedy, @jcea, @mdickinson, @vstinner
Files	cmathmodule.c.2.6a3.diff mathmodule.c.2.6a3.diff test_math_sum1.py mathmodule.c.2.6a3.take2.diff floatsum.c msum.py: Version of msum that takes care of overflow msum3.py: Modified msum.py with overflow file crsum.py: Correctly rounded sum of nonadjacent floats msum4.py: msum, with overflow, correct rounding and special-value handling mathmodule4.c.2.6a3.diff: msum4 patch for mathmodule.c of Python 2.6a3 cmathmodule4.c.2.6a3.diff: msum4 patch for cmathmodule.c of Python 2.6a3 test_math_sum4.py: msum4.py tests for math- and cmathmodule.c mathmodule5.c.2.6a3.diff: update Convert READMEs to rst #5 of mathmodule.c for Python 2.6a3 cmathmodule5.c.2.6a3.diff: update Convert READMEs to rst #5 of cmathmodule.c for Python 2.6a3 test_math_sum5.py: update Convert READMEs to rst #5 of the test script mathmodule11.c.2.6a3.diff: rev 11 patch for mathmodule.c for Python 2.6a3 test_math_sum11.py: rev 11 of test_math_sum.py for Python 2.6a3 mathmodule12.c.2.6a3.diff: rev 12 of test_math_sum.py for Python 2.6a3 test_math_sum12.py: rev 12 of test_math_sum.py for Python 2.6a3 mathmodule19.c.2.6a3.diff: rev 19 patch for mathmodule.c for Python 2.6a3 test_math_sum19.py: rev 19 of test_math_sum.py for Python 2.6a3 mathsum_IA32.patch fsum11.patch

^{Note: these values reflect the state of the issue at the time it was migrated and might not reflect the current state.}

Show more details

GitHub fields:

assignee = 'https://github.com/mdickinson'
closed_at = <Date 2008-12-11.19:54:45.020>
created_at = <Date 2008-05-11.16:16:51.847>
labels = ['extension-modules', 'performance']
title = 'Full precision summation'
updated_at = <Date 2008-12-11.19:54:45.008>
user = 'https://bugs.python.org/MrJean1'

bugs.python.org fields:

activity = <Date 2008-12-11.19:54:45.008>
actor = 'mark.dickinson'
assignee = 'mark.dickinson'
closed = True
closed_date = <Date 2008-12-11.19:54:45.020>
closer = 'mark.dickinson'
components = ['Extension Modules']
creation = <Date 2008-05-11.16:16:51.847>
creator = 'MrJean1'
dependencies = []
files = ['10285', '10286', '10287', '10307', '10308', '10345', '10347', '10351', '10357', '10371', '10372', '10373', '10379', '10380', '10381', '10405', '10406', '10407', '10408', '10410', '10411', '10505', '11108']
hgrepos = []
issue_num = 2819
keywords = ['patch']
message_count = 47.0
messages = ['66641', '66655', '66663', '66666', '66734', '66745', '66749', '66750', '66753', '66755', '66764', '66844', '66947', '66953', '66954', '66965', '67006', '67021', '67034', '67045', '67053', '67067', '67080', '67197', '67198', '67210', '67211', '67212', '67213', '67647', '70295', '70317', '70409', '70424', '70430', '70438', '70439', '70440', '70508', '70522', '70545', '70585', '70611', '70623', '70641', '71118', '77628']
nosy_count = 7.0
nosy_names = ['tim.peters', 'rhettinger', 'terry.reedy', 'jcea', 'mark.dickinson', 'vstinner', 'MrJean1']
pr_nums = []
priority = 'normal'
resolution = 'accepted'
stage = None
status = 'closed'
superseder = None
type = 'performance'
url = 'https://bugs.python.org/issue2819'
versions = ['Python 2.6']

[MrJean1]

Attached are 2 patches and a test script adding a function sum to the
math and cmath modules of Python 2.6a3. The sum is calculated using a
full precision summation method.

The test script compares the result of the functions with the original
implementation in Python. All tests pass with 4 different builds of
Python 2.6a3:

• GNU gcc 4.0.1 on MacOS X 10.4.11 (Intel Core Duo), 32-bit
• GNU gcc 4.1.2 on RHEL 3 update 7 (Opteron), 64-bit
• Sun C 5.8 on Solaris 10 (Opteron), both 32- and 64-bit

[rhettinger]

This looks pretty good at first glance. Will review more throughly
later this week. It does need docs and unittests.

[mdickinson]

Some comments/questions:

(1) It seems wasteful to wrap every addition in PyFPE_START/END_PROTECT,
and to check for NaNs and infinities after every addition. I'd wrap the
whole thing in a single PyFPE_START/END_PROTECT, replace _math_sum_add
with an inline addition, and just let NaNs and Infs sort themselves out.

If the result comes out finite (as it surely will in almost all
applications), then all the summands were necessarily finite and there's
nothing more to do. If the result comes out as an infinity or NaN, you
need to decide whether it's appropriate to return a NaN, an infinity, or
to raise OverflowError or ValueError. I'm not sure it's worth trying to
do the right thing for all special value cases, but if you do want to
follow 'spirit of IEEE 754' rules for special values, they should look
something like this:

(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.)

Note that some sums involving overflow won't be computed correctly:
e.g. [1e308, 1e308, -1e308] will likely sum to infinity instead of
returning 1e308. I don't see any easy way around this, and it's
probably not worth worrying about.

(2) The algorithm is only guaranteed to work correctly assuming IEEE 754
semantics. Python currently doesn't insist on IEEE 754 floating point,
so what should happen on non IEEE-754 machines?

(3) Rather than duplicating the math module code in cmathmodule.c, why
not have the complex version simply sum real parts and imaginary parts
separately, using a version of the code that's already in mathmodule.c?

[mdickinson]

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.

[MrJean1]

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 <report@bugs.python.org>
wrote:

Mark Dickinson <dickinsm@gmail.com> 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 <report@bugs.python.org>
<http://bugs.python.org/issue2819\>

---

[MrJean1]

My apologies for the messy post. I replied to the email instead of
posting my response.

/Jean Brouwers

PS) Attached is *an* example of the math_sum() and cmath_sum() functions
using the same, shared function float_sum(). Perhaps, that resides in
Objects/floatobject.c?

[mdickinson]

One related issue is testing these, how can a NaN and +/-Infinity
float object be created in Python?

In 2.6 and 3.0 (but not 2.5 and older), float('nan'), float('inf') and
float('-inf') should all work reliably across platforms (or at least
those platforms that support infs and nans). If they don't it's a bug.

(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?

Actually, I think you could probably just leave the algorithm exactly as
it is, but put a warning in the documentation that the exactness only
applies in the presence of IEEE 754 semantics. Practically everybody's
on an IEEE 754 platform anyway.

There are two symbols in pyconfig.h HAVE_IEEEFP_H and
HAVE_LIBIEEE. Are those the proper ones to determine IEEE floating
point support?

I'm not sure that either of these is the right thing. Neither is
defined on my MacBook, for example.

(3) On the 3rd comment, Raymond and I did discus using a single
function to be called by the math and cmath modules.

I think you're right that it's easier to just duplicate the code.
It's a nice feature that this function only has to pass once through the
data, and it doesn't seem worth losing that to save a little bit of code
duplication.

I still wonder whether there's a way to avoid incorrectly signaling
overflow in the case where the result is finite. (e.g. sum([1e308,
1e308, -1e308])).

[MrJean1]

It turns out, float('nan') creates a Nan and float('[+|-]inf') creates a
[signed] Infinity object. That answers my earlier question.

/Jean Brouwers

[MrJean1]

The current results are quite "interesting"

>>> math.sum([1e308, 1e308, -1e308])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
OverflowError: math range error

>>> math.sum([1e308, -1e308, 1e308])
1e+308

Handling this case would require holding the finite, prior sum at an
intermediate overflow in a separate array, say overflows. Then, add the
partials which may create additional overflows. Finally, keep adding
the overflows (accurately?) until none remain or until none can be added
without overflow.

/Jean Brouwers

[MrJean1]

There may be another reason to use a single summation function. The
summation functions does need (a copy of) the is_error() function from
the math module.

The cmath module has a similar function called math_error() which
slightly different from is_error().

It would be better, more consistent if both modules used the same
function, say is_error() moved** to file Object/floatobject.c

Then, the math and cmath module can use ist_error() instead of each
their own. Also, the summation function can use it and function
float_pow() in floatobject.c could.

/Jean Brouwers

**) and renamed to e.g. float_error().

[rhettinger]

Rather keep it in mathmodule.c, not floatobject.c.
We should keep extension code out of the core types.

[rhettinger]

When you need full precision, the Kahan approach helps but doesn't make
guarantees and can sometimes hurt (it makes some assumptions about the
data). One use case in is computing stats like a mean where many of
the larger magnitude entries tend to cancel out but only after clipping
bits off of the lower magnitude components.