update: based on what others have found in the source, I was wrong about this - sum() does not sort. The patterns of consistency I found below stem from the fact that sorting (as done in some cases below) and using extended-precision intermediate values (as done in sum()) can have similar effects on precision ...
@user2357112 comments below:
src/main/summary.c ... doesn't do any sorting. (That'd be a lot of expense to add to a summation operation.) It's not even using pairwise or compensated summation; it just naively adds everything up left to right in an LDOUBLE (either long double or double, depending on HAVE_LONG_DOUBLE).
I have exhausted myself looking for this in the R source code (without success - sum is hard to search for), but I can show by experiment that when executing sum(), R sorts the input vector from smallest to largest in order to maximize accuracy; the difference between sum() and Reduce() results below is due to use of extended precision. I don't know what accu does ...
set.seed(101)
vec <- runif(100, 0, 0.00001)
options(digits=20)
(s1 <- sum(vec))
## [1] 0.00052502325481269514554
Using Reduce("+",...) just adds the elements in order.
(s2 <- Reduce("+",sort(vec)))
## [1] 0.00052502325481269514554
(s3 <- Reduce("+",vec))
## [1] 0.00052502325481269503712
identical(s1,s2) ## TRUE
?sum() also says
Where possible extended-precision accumulators are used, but this is platform-dependent.
Doing this in RcppArmadillo on the sorted vector gives the same answer as in R; doing it on the vector in the original order gives yet a different answer (I don't know why; my guess would be the aforementioned extended-precision accumulators, which would affect the numerical outcome more when the data are unsorted).
suppressMessages(require(inline))
code <- '
arma::vec ax = Rcpp::as<arma::vec>(x);
return Rcpp::wrap(arma::accu(ax));
'
## create the compiled function
armasum <- cxxfunction(signature(x="numeric"),
code,plugin="RcppArmadillo")
(s4 <- armasum(vec))
## [1] 0.00052502325481269525396
(s5 <- armasum(sort(vec)))
## [1] 0.00052502325481269514554
identical(s1,s5) ## TRUE
But as pointed out in comments this doesn't work for all seeds: in this case the Reduce() result is closer to the results of sum()
set.seed(123)
vec2 <- runif(50000,0,0.000001)
s4 <- sum(vec2); s5 <- Reduce("+",sort(vec2))
s6 <- Reduce("+",vec2); s7 <- armasum(sort(vec2))
rbind(s4,s5,s6,s7)
## [,1]
## s4 0.024869900535651481843
## s5 0.024869900535651658785
## s6 0.024869900535651523477
## s7 0.024869900535651343065
I'm stumped here. I would have expected at least s6 and s7 to be identical ...
I will point out that in general when your algorithm depends on these kinds of tiny numeric differences you're likely to be getting very frustrated, as the results are likely to differ on the basis of many small and possibly-out-of-your-control factors like particular operating system, compiler, etc. you work with.
