[Rd] Discrepancy: R sum() VS C or Fortran sum

[Pierre Chausse]

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.

-- 
Pierre Chaussé
Assistant Professor
Department of Economics
University of Waterloo