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 Bill Dunlap TIBCO Software wdunlap tibco.com On Tue, Feb 19, 2019 at 10:36 AM Paul Gilbert <pgilbert...@gmail.com> wrote: > (I didn't see anyone else answer this, so ...) > > You can probably find the R code in src/main/ but I'm not sure. You are > talking about a very simple calculation, so it seems unlike that the > algorithm is the cause of the difference. I have done much more > complicated things and usually get machine precision comparisons. There > are four possibilities I can think of that could cause (small) differences. > > 0/ Your code is wrong, but that seems unlikely on such a simple > calculations. > > 1/ You are summing a very large number of numbers, in which case the sum > can become very large compared to numbers being added, then things can > get a bit funny. > > 2/ You are using single precision in fortran rather than double. Double > is needed for all floating point numbers you use! > > 3/ You have not zeroed the double precision numbers in fortran. (Some > compilers do not do this automatically and you have to specify it.) Then > if you accidentally put singles, like a constant 0.0 rather than a > constant 0.0D+0, into a double you will have small junk in the lower > precision part. > > (I am assuming you are talking about a sum of reals, not integer or > complex.) > > HTH, > Paul Gilbert > > On 2/14/19 2:08 PM, Rampal Etienne wrote: > > Hello, > > > > I am trying to write FORTRAN code to do the same as some R code I have. > > I get (small) differences when using the sum function in R. I know there > > are numerical routines to improve precision, but I have not been able to > > figure out what algorithm R is using. Does anyone know this? Or where > > can I find the code for the sum function? > > > > Regards, > > > > Rampal Etienne > > > > ______________________________________________ > > R-devel@r-project.org mailing list > > https://stat.ethz.ch/mailman/listinfo/r-devel > > ______________________________________________ > R-devel@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel > [[alternative HTML version deleted]] ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel