> > On 2019-02-19 2:08 p.m., William Dunlap via R-devel wrote: >> The algorithm does make a differece. You can use Kahan's summation >> algorithm (https://en.wikipedia.org/wiki/Kahan_summation_algorithm) to >> reduce the error compared to the naive summation algorithm. E.g., in R >> code: >> >> naiveSum <- >> function(x) { >> s <- 0.0 >> for(xi in x) s <- s + xi >> s >> } >> kahanSum <- function (x) >> { >> s <- 0.0 >> c <- 0.0 # running compensation for lost low-order bits >> for(xi in x) { >> y <- xi - c >> t <- s + y # low-order bits of y may be lost here >> c <- (t - s) - y >> s <- t >> } >> s >> } >> >>> rSum <- vapply(c(1:20,10^(2:7)), function(n) sum(rep(1/7,n)), 0) >>> rNaiveSum <- vapply(c(1:20,10^(2:7)), function(n) naiveSum(rep(1/7,n)), 0) >>> rKahanSum <- vapply(c(1:20,10^(2:7)), function(n) kahanSum(rep(1/7,n)), 0) >>> >>> table(rSum == rNaiveSum) >> >> FALSE TRUE >> 21 5 >>> table(rSum == rKahanSum) >> >> FALSE TRUE >> 3 23
If you use the vector c(1,10^100,1,-10^100) as input then sum, naiveSum or kahanSum will all give an incorrect answer. All return 0 instead of 2. >From the wikipedia page we can try the pseudocode given of the modification by >Neumaier. My R version (with a small correction to avoid cancellation?) is neumaierSum <- function (x) { s <- 0.0 z <- 0.0 # running compensation for lost low-order bits for(xi in x) { t <- s + xi if( abs(s) >= abs(xi) ){ b <- (s-t)+xi # intermediate step needed in R otherwise cancellation z <- z+b # If sum is bigger, low-order digits of xi are lost. } else { b <- (xi-t)+s # intermediate step needed in R otherwise cancellation z <- z+b # else low-order digits of sum are lost } s <- t } s+z # correction only applied once in the very end } testx <- c(1,10^100,1,-10^100) neumaierSum(testx) gives 2 as answer. Berend Hasselman ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel