I noticed the output of sum() is different from custom loops

julia> function mysum(A)
out = 0.0
for ii in eachindex(A)
out += A[ii]
end
return out
end
mysum (generic function with 1 method)
julia> A=rand(1000);
julia> sum(A)-mysum(A)
-1.1368683772161603e-13

In FORTRAN(gfortran) this only happens when -ffastmath is enabled. Is this by design, or there is something else I missed? If this is expected, what other functions have this behavior?

Good to know. Thanks.
And why does sum() sometimes not agree with itself either?

julia> function rowsum(A)
out = Array{Float64,2}(undef, size(A, 1),1)
for ii in axes(A, 1)
out[ii,1] = sum(A[ii, :])
end
return out
end
rowsum (generic function with 1 method)
julia> A=rand(10,1000);
julia> sum(A,dims=2)-rowsum(A)
10×1 Matrix{Float64}:
2.2737367544323206e-13
2.2737367544323206e-13
4.547473508864641e-13
-2.2737367544323206e-13
7.389644451905042e-13
1.7053025658242404e-13
-2.8421709430404007e-13
-2.2737367544323206e-13
-5.115907697472721e-13
5.684341886080802e-14

In the case where you have a large number of similar-magnitude integers (so their sum is BigInt territory) I think pairwise summation is asymptotically faster than sequential

That is interesting to me. I would have thought that base sum always used Kahan summation, at least on items with known size. I wonder why that is not the case.

Ah sorry, I tend to always equate the term pairwise sum with Kahan sum. Anyways, what I meant is I’m wondering why summing along other dimensions would be any different than a standard sum.

Not entirely sure. My guess would be that it’s just that it’s easier to write it as a maprecuce than writing the code to do pairwise summation over arbitrary dimensions.

In general, the reduction operations in Julia (mapreduce, sum, …) are documented to treat the operator as associative, and to use an implementation-dependent associativity. If you want a particular associativity you have to use foldl/foldr or similar (or write your own loop).