Install the gmp package, run your code, and then try this:
bu <- gmp::as.bigq(u)
bs4 <- bu[1] + bu[2] + bu[3] + bu[4] + bu[5]
s4 <- as.double(bs4)
s1 - s4
## [1] 0
s2[[2]] - s4
## [1] 7.105427e-15
s3 - s4
## [1] 7.105427e-15
identical(s1, s4)
## [1] TRUE
`bs4` is the exact sum of the binary rationals in your `u` vector;
`s4` is the closest double precision to this exact sum.
Looking at the C source code for sum() will show you that it makes
some extra efforts to get a more accurate sum than your simple
version.
Best,
luke
On Fri, 16 Mar 2018, Pierre Chausse wrote:
Hi all,
I found a discrepancy between the sum() in R and either a sum done in C or
Fortran for vector of just 5 elements. The difference is very small, but this
is a very small part of a much larger numerical problem in which first and
second derivatives are computed numerically. This is part of a numerical
method course I am teaching in which I want to compare speeds of R versus
Fortran (We solve a general equilibrium problem all numerically, if you want
to know). Because of this discrepancy, the Jacobian and Hessian in R versus
in Fortran are quite different, which results in the Newton method producing
a different solution (for a given stopping rule). Since the solution computed
in Fortran is almost identical to the analytical solution, I suspect that the
sum in Fortran may be more accurate (That's just a guess). Most of the time
the sum produces identical results, but for some numbers, it is different.
The following example, shows what happens:
set.seed(12233)
n <- 5
a <- runif(n,1,5)
e <- runif(n, 5*(1:n),10*(1:n))
s <- runif(1, 1.2, 4)
p <- runif(5, 3, 10)
x <- c(e[-5], (sum(e*p)-sum(e[-5]*p[-5]))/p[5])
u <- a^(1/s)*x^((s-1)/s)
dyn.load("sumF.so")
u[1] <- u[1]+.0001 ### If we do not add .0001, all differences are 0
s1 <- sum(u)
s2 <- .Fortran("sumf", as.double(u), as.integer(n), sf1=double(1),
sf2=double(1))[3:4]
s3 <- .C("sumc", as.double(u), as.integer(n), sC=double(1))[[3]]
s1-s2[[1]] ## R versus compiler sum() (Fortran)
[1] -7.105427e-15
s1-s2[[2]] ## R versus manual sum (Fortran
[1] -7.105427e-15
s1-s3 ## R Versus manual sum in C
[1] -7.105427e-15
s2[[2]]-s2[[1]] ## manual sum versus compiler sum() (Fortran)
[1] 0
s3-s2[[2]] ## Fortran versus C
[1] 0
My sumf and sumc are
subroutine sumf(x, n, sx1, sx2)
integer i, n
double precision x(n), sx1, sx2
sx1 = sum(x)
sx2 = 0.0d0
do i=1,n
sx2 = sx2+x(i)
end do
end
void sumc(double *x, int *n, double *sum)
{
int i;
double sum1 = 0.0;
for (i=0; i< *n; i++) {
sum1 += x[i];
}
*sum = sum1;
}
Can that be a bug? Thanks.
--
Luke Tierney
Ralph E. Wareham Professor of Mathematical Sciences
University of Iowa Phone: 319-335-3386
Department of Statistics and Fax: 319-335-3017
Actuarial Science
241 Schaeffer Hall email: luke-tier...@uiowa.edu
Iowa City, IA 52242 WWW: http://www.stat.uiowa.edu
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
