New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Full precision summation #47068
Comments
Attached are 2 patches and a test script adding a function sum to the The test script compares the result of the functions with the original
|
This looks pretty good at first glance. Will review more throughly |
Some comments/questions: (1) It seems wasteful to wrap every addition in PyFPE_START/END_PROTECT, If the result comes out finite (as it surely will in almost all (1) if the summands include a NaN, return a NaN Note that some sums involving overflow won't be computed correctly: (2) The algorithm is only guaranteed to work correctly assuming IEEE 754 (3) Rather than duplicating the math module code in cmathmodule.c, why |
One more question: What are the use cases for an exact summation algorithm? That is, in what |
Mark, (1) Attached is an attempt to address the 1st issue (just the One related issue is testing these, how can a NaN and +/-Infinity float (2) On your 2nd comment, supporting non-IEEE floating point, perhaps the (3) On the 3rd comment, Raymond and I did discus using a single function to Also, depending on the implementation of that function, it may require /Jean Brouwers On Sun, May 11, 2008 at 1:49 PM, Mark Dickinson <report@bugs.python.org>
|
My apologies for the messy post. I replied to the email instead of /Jean Brouwers PS) Attached is *an* example of the math_sum() and cmath_sum() functions |
In 2.6 and 3.0 (but not 2.5 and older), float('nan'), float('inf') and
Actually, I think you could probably just leave the algorithm exactly as
I'm not sure that either of these is the right thing. Neither is
I think you're right that it's easier to just duplicate the code. I still wonder whether there's a way to avoid incorrectly signaling |
It turns out, float('nan') creates a Nan and float('[+|-]inf') creates a /Jean Brouwers |
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 /Jean Brouwers |
There may be another reason to use a single summation function. The The cmath module has a similar function called math_error() which It would be better, more consistent if both modules used the same Then, the math and cmath module can use ist_error() instead of each /Jean Brouwers **) and renamed to e.g. float_error(). |
Rather keep it in mathmodule.c, not floatobject.c. |
When you need full precision, the Kahan approach helps but doesn't make |
Here's (msum.py) an example in Python of one fairly straightforward way of The idea is essentially just to maintain the sum modulo integer multiples I'm 97.3% convinced that the proof of correctness goes through: it's --- There's one more 'nice-to-have' that I think should be considered: it msum(list) == msum(sorted(list)) would hold true. |
Two tests failed with Python 2.6a3 on MacOS X Intel. Test 11 failed: 9007199254740992.0 vs 9007199254740991.0 expected for Test 12 failed: inf vs 1.7976931348623157e+308 expected for /Jeab Brouwers |
Yes: that's the lack of correct rounding rearing its ugly head...
I'm still trying to work out how to get around this one; again, the |
Ensuring correct rounding isn't as onerous as I expected it to be. |
Okay, just to show it's possible: Here (msum4.py) is a modified version of Raymond's recipe that deals (1) intermediate overflows The file contains more tests, and a proof of correctness. The algorithm |
I intend to submit a C version of msum4 shortly. |
Attached is the patch for the the mathmodule.c file of Python 2.6a3 Please review the C code, in particular the setting of _do_sum_pow_2 for The results of the C and Python versions match (on 32-bit MacOS X More testing and the cmath version are still to be done. /Jean Brouwers |
Attached are updated patches for both the mathmodule.c and cmathmodule.c |
The tests of the test_math_sum4.py script all pass with the latest math-
Also interesting might be that the tests run 25+ times faster with the C /Jean Brouwers |
The latest mathmodule patch (file 10371) looks pretty good to me. I
|
Here is another patch for the mathmodule.c and cmathmodule.c files updated An update of the test script is also attached. The tests all pass on the /Jean Brouwers |
Attached is revision 11 of the mathmodule.c patch for Python 2.6a3. This one includes Raymond's full precision summation, Mark's rounding However, intermediate overflow will raise an OverflowError and a An updated test_math_sum11.py script is also attached. All test cases
/Jean Brouwers |
Here is rev 12 of the mathmodule.c patch. It is the same as rev 11 but
The test_math.sum12.py script and results are the same as rev 11. /Jean Brouwers |
Here is another, cleaner revision 19 of the same mathmodule.c patch and /Jean Brouwers |
Nice work Jean. Marking the patch accepted. Mark please go ahead with Once the commit has settled for a couple of days, go ahead with a After that one settles, close this issue and open a new one for a |
math module patch committed, r63542. I'm still working on converting the tests to the unittest framework. |
Tests committed in r63543 |
Here's a patch that fixes the problems with math.sum on x86 machines that On these machines, the patched version of math.sum does almost the entire |
Here's a patch giving an alternative implementation of math.fsum; it's The patch is experimental: it doesn't have sufficient tests, has no On my MacBook, math.fsum is a factor of 2-3 times slower than math.sum. We *really* need to sort math.sum out, one way or another, before the Some options: (1) leave math.sum as it is, skip all tests on x86/Linux, and document (2) investigate a version of math.sum that plays with the FPU control (3) replace math.sum with the slower but (presumably) less erratic Jean, Raymond: what do you think? |
Tests related to overflow, special-values, intermediate overflow, and |
See r65292 for more updates to the test-suite: I've replaced the Python |
Minor code cleanups, and fixes to special-value handling in r65299 |
Renamed math.sum to math.fsum (as previously discussed) in r65308. I think all that's left now is to add a note to the docs about the |
Added documentation note about x86 problems in r65315. Jean, Raymond: is it okay to close this issue now? |
Here (fsum8.patch) is a clean version of the alternative fsum algorithm. I'd like to push for using this in place of the existing algorithm. |
Mark, I don't currently have a machine with SVN and a compiler Something that's surprised me for decades is how slow platform ldexp() |
[Tim]
Indeed it would! Here's a revised patch that avoids use of fmod. Speed is comparable with the current
------------------------------------------------+----------+--------- Note that lsum doesn't suffer from the 'fragmentation of partials' problem that slows down msum for sorted I'll do some timings on Linux/x86 as well. |
Timings on x86/Linux are similar: the lsum-based version is around There's probably still some snot left to optimize out, though. Some (1) to try using doubles instead of longs for the accumulator digits Anything else? |
Toned down note in docs (at Raymond's request) in r65366. While I'd really like to see an lsum-based version go in instead, I There was also originally some talk of a complex math version. What do |
Mark, changing API is something to avoid after beta cycles begin (so, The big attraction to the lsum-like approach for users is that it "just |
Comment: as a user (on x86) I agree with Tim. I could almost see the Question: I see "math module patch committed, r63542" in May 22. But in
3.0b2, there is no math.fsum and math.sum seems to be a wrapper for
builtin sum. Has this not been forward ported (merged) yet?
>>> math.sum
<built-in function sum>
>>> sum
<built-in function sum>
>>> math.sum is sum
False |
I'm pretty sure it *was* merged: math.sum should be the full-precision The name change to fsum only happened in the last few days. |
They do. I realize now that two different built-in funcs in two
Which will prevent the confusion I had ;-). Good idea. |
Here's a patch, in final form, that replaces fsum with an lsum-based huge_integer * 2**(smallest_possible_exponent) and the huge_integer is stored in base 2**30, with a signed-digit representation (digits What are the chances of getting this in before next week's beta? I did toy with a base 2**52 version, with digits stored as doubles. It's attractive for |
Closing this. Let's stick with what we have. |
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:
bugs.python.org fields:
The text was updated successfully, but these errors were encountered: