`sum` float precision issue

[anton3s]

Bug report

Bug description:

I am encountering some strange results in Python 3.12.3 when summing up numbers using the sum function.

It seems that there are some precision issues, which I wouldn't expect with these numbers:

>>> 2.5392 + 0.4608 - 3.0
0.0
>>> sum([2.5392,0.4608,-3.0], 0.0)
1.1102230246251565e-16
>>> import numpy
>>> sum([2.5392,0.4608,-3.0], numpy.float64(0.0))
0.0

I tested on the official docker images of Python 3.11 - 3.13 and version 3.13 also seem affected (3.11 works as expected).

$ docker run -it  python:3.11 python -c "print(sum([2.5392,0.4608,-3.0], 0.0))"
0.0
$ docker run -it  python:3.12 python -c "print(sum([2.5392,0.4608,-3.0], 0.0))"
1.1102230246251565e-16
$ docker run -it  python:3.13 python -c "print(sum([2.5392,0.4608,-3.0], 0.0))"
1.1102230246251565e-16

CPython versions tested on:

3.12, 3.13

Operating systems tested on:

Linux, macOS

[Eclips4]

Duplicate of #111933

[skirpichev]

So, in short:

1. sum(lst) != reduce(lambda x, y: x+y, lst), i.e. (2.5392 + 0.4608) + (-3.0) in your case
2. since CPython 3.12, sum() uses improved variant of compensated summation; more accurate (but slow!) version is in the math.fsum()
3. numpy uses less accurate approach (partial pairwise summation)

Answer 0.0 is not accurate in fact, you can easily do exact calculations for given floats:

>>> from fractions import Fraction as F
>>> lst = [2.5392, 0.4608, -3.0]
>>> lst_F = list(map(F.from_float, lst))
>>> sum(lst_F)  # nonzero!
Fraction(1, 9007199254740992)
>>> sum(lst_F) - F.from_float(sum(lst))  # sum(lst) produces exact answer here
Fraction(0, 1)

[anton3s]

Yes, I understand that I shouldn't do equality comparison with floats.

The problem I have is that those three numbers in my example come from another piece of my program that runs an optimization (using scipy.optimize) with the constraint that the sum of those three numbers is <= 0. But when I then check that constraint in Python, it doesn't hold.

The issue with this is that sum([a,b,c]) != a + b + c which breaks any intuition of what sum means mathematically.
I found it unimaginable that this was intended, but given the outcome of the lengthy discussion in issue #126181, I see that it is, so I won't open up that discussion again.

[skirpichev]

I shouldn't do equality comparison with floats.

That's a bad idea for any computed floating-point values, not just sums.

with the constraint that the sum of those three numbers is <= 0

I would guess, that this constraint is about the sum of three real numbers. Floating-point numbers aren't real numbers from analysis course.

The issue with this is that sum([a,b,c]) != a + b + c which breaks any intuition of what sum means mathematically.

How you intuition live with the fact that (a + b) + c != a + (b + c) in general? ;)

>>> (0.1 + 0.2) + 0.3
0.6000000000000001
>>> 0.1 + (0.2 + 0.3)
0.6
>>> 0.2 + (0.1 + 0.3)
0.6000000000000001

This is not a worst case. It's easy to find floating-point numbers, where all three possible sums (as commutativity holds for addition, only associativity is broken) - are different:

>>> while True:
...     a, b, c = [random.random() for _ in range(3)]
...     ss = [(a + b) + c, a + (b + c), b + (a + c)]
...     if ss[0] != ss[1] and ss[1] != ss[2] and ss[0] != ss[2]:
...         print([a, b, c])
...         print(ss)
...         break
...         
[0.457290196955373, 0.7565955400028612, 0.38914139486258204]
[1.6030271318208165, 1.603027131820816, 1.6030271318208162]

[anton3s]

Ok, point taken. I need to revise my intuitions. ;)

[rhettinger]

On the plus side (pun intended), your intuitions about commutatively are now likely to be correct.

>>> Counter([len(set(map(sum, permutations([random(), random(), random()]))))
...     for i in range(1000)])
...
Counter({1: 1000})

Formerly, this would have returned something like Counter({1: 751, 2: 225, 3: 24}) showing that 25% of random triplets had different sums depending on the order of their arguments.

[skirpichev]

showing that 25% of random triplets had different sums depending on the order of their arguments.

And ~60% already for four arguments.

[Tevemadar]

The '90s would have been a good time to implement this algorithm in Python, early 2000s too. The article was 20 years old by then. Python did not implement it for the 30th and 40th anniversaries either.
As Python gained some momentum as being a "scientific" language during this time, I'd like to point out that the scientific community values both consistency and reproducibility, and +, sum() creating different results now, compared to what they did for circa 30 years, goes against that. While it may sound strange, there also can be people who don't often re-read the documentation of sum() after using it for decades. This 50 years old novelty could have very well been an addition to the math module.

For completeness this example includes Numpy too, as it seems to be affected somehow - or it's just evolving on its own too.
Results on 3.12:

[0.0, 0.06108963827935576] <- this is the 'original'
[0.0, 0.06108963827935582] <- sum()
[0.0, 0.06108963827935587] <- numpy
0.06108963827935576 <- numpy 1D

(IdeOne uses 3.9 at the time of posting, and all 4 displayed results are the 'original', ending with 76)

[skirpichev]

The article was 20 years old by then.

Euclidean algorithm is ~2300 years old. Sometimes it's the best we can have.

and +, sum() creating different results now

a + b and sum([a, b]) is same. (a + b) + c and sum([a, b, c]) - isn't (ditto for more than 3 operands).

But if you want (a + b) + c-like "summation" - it's always possible to do this explicitly, with a loop or using reduce(). I don't see any problems with reproducibility: if you want addition with some given order of operands - you can get it.

creating different results now, compared to what they did for circa 30 years

sum() docs since v3.0 says "Sums start and the items of an iterable from left to right and returns the total. start defaults to 0. [...]". This is still correct. Docs doesn't says how exactly summation is performed.

This 50 years old novelty could have very well been an addition to the math module.

The math module has fsum(). It uses a different algorithm and has better error bounds.

For completeness this example includes Numpy too

Numpy uses partial pairwise summation, a bit faster than fsum(), but less precise.

[tim-one]

I'd like to point out that the scientific community values both consistency and reproducibility

It does. It also values accuracy and speed and portability. And different segments of "the scientific community" have different priorities among these 5. For example, short- and medium-range weather forecasting models favor speed most. They have an enormous amount of low-precision data to crunch, and it's no use to "predict" weather after it's already happened 😉

I vividly recall that in the very first Cray Research vector "supercomputers", multiplication wasn't even commutative: x*y could compute a different result than y*x. Purely so for speed. But that was such a huge hit to consistency and portability than even Seymour Cray was shamed 😉 into admitting that was a design error. By the time the 3rd machine was built, float multiplication was commutative (but still not "accurate" by modern standards: Mr. Cray never thought it was worth the time to compute the full multiplication tree in hardware - and Cray customers were willing to live with that tradeoff).

So it goes. All platforms have tradeoffs to make, and they vary across time too.

Python isn't actually "a scientific language". The numpy extension does the heavy float lifting for most serious float applications, and favors speed more highly than the core Python language. Core Python in turn has tended to favor accuracy and reproducibility more, but things like math.fsum() carry that to such an extreme that they're just impractically slow for many applications. So core Python changes too over time.

In their current states, CPython and numpy both try to improve sum() accuracy over the naïve "strictly left to right` approach at a cost in speed, but typically enough numpy's approach is faster (but also less accurate) than CPython's.

For the vast bulk of users, I think CPython made the right tradeoffs in its current sum(). But if, e.g., another way was found tomorrow that materially improved accuracy without hurting speed, or materially improved speed without hurting accuracy, CPython should move to that instead.

[Tevemadar]

Numpy uses partial pairwise summation, a bit faster than fsum(), but less precise.

Off-topic here, but numpy and numpy 1D lines are different now, so multiplying and summing the same 3 pairs of numbers produces different results for vector-matrix and vector-vector multiplication, à la the Cray-anecdote. Some may say that such behaviour is unexpected.

For us, these changes result in visibly different outputs in an intermediate step of an image processing pipeline that tells apart various brain regions and counts cell bodies. Then of course some test cases suddenly fail. Code is now back to a for loop as the given step has to be consistent with 3 other languages (C++, Java, JavaScript). So in our specific context it's Python what's got broken.

[skirpichev]

I would guess, your tests check either repr() output or compare computed (by sum()) floating-point values with ==. If so, this looks rather as an error in your code: you should use instead math/numpy.isclose() or format() your output with appropriate precision settings.

Will you assume, that mathematical functions in the math module will always return exactly same output for a given input? Hardly so, even for elementary functions there are many algorithms, having different characteristics. New version of the libm might come with a new one. Or this might vary for different C stdlib implementations. The sum() case is just like that. As math.sin() or math.dist().