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
See do_summary() in summary.c, rsum() for doubles. R uses long double type as accumulator on systems where available. Best, Tomas 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 at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel
(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 at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel
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 <pgilbert902 at 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 at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-devel
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
[[alternative HTML version deleted]]
Dear Paul, Thank you for thinking with me. I will respond to your options:> > 0/ Your code is wrong, but that seems unlikely on such a simple > calculations. >It's really just a comparison of the sum function in Fortran with that of R. If instead I use the naive summation with a for loop in both languages I get the same answer.> > 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. >Yes, this is what's happening and why I need to know what algorithm R uses to overcome or reduce these precision issues.> > 2/ You are using single precision in fortran rather than double. Double > is needed for all floating point numbers you use! >I use doubles.> > 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. >It doesn't matter if I double them in this way or not (they are declared as doubles and I think the compiler treats then as doubles). So my question remains what algorithm R uses. Cheers, Rampal> > 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 at r-project.org mailing list > > https://stat.ethz.ch/mailman/listinfo/r-devel >[[alternative HTML version deleted]]
Dear Tomas, Where do I find these files? Do they contain the code for the sum function? What do you mean exactly with your point on long doubles? Where can I find documentation on this? Cheers, Rampal On Mon, Feb 18, 2019, 15:38 Tomas Kalibera <tomas.kalibera at gmail.com wrote:> See do_summary() in summary.c, rsum() for doubles. R uses long double > type as accumulator on systems where available. > > Best, > Tomas > > 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 at r-project.org mailing list > > https://stat.ethz.ch/mailman/listinfo/r-devel > > >[[alternative HTML version deleted]]