Different Float64 sum on different architectures

Here is a fun example of the lack of associativity (thanks to Stefan Karpinski):

function sumsto(x::Float64)
    0 <= x < exp2(970) || throw(ArgumentError("sum must be in [0,2^970)"))
    n, p₀ = Base.decompose(x) # integers such that `n*exp2(p₀) == x`
    [floatmax(); [exp2(p) for p in -1074:969 if iseven(n >> (p-p₀))]
    -floatmax(); [exp2(p) for p in -1074:969 if isodd(n >> (p-p₀))]]
end

Example:

julia> x = sumsto(2.3);

julia> y = sumsto(1e18);

julia> sort(x) == sort(y)
true

julia> foldl(+,x)
2.3

julia> foldl(+,y)
1.0e18
4 Likes