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

luke-tierney at uiowa.edu luke-tierney at uiowa.edu
Fri Mar 16 16:37:03 CET 2018

```Install the gmp package, run your code, and then try this:

bu <- gmp::as.bigq(u)
bs4 <- bu + bu + bu + bu + bu
s4 <- as.double(bs4)
s1 - s4
##   0
s2[] - s4
##   7.105427e-15
s3 - s4
##   7.105427e-15
identical(s1, s4)
##   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)
> u <- a^(1/s)*x^((s-1)/s)
>
> u <- u+.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))[]
>
> s1-s2[] ## R versus compiler sum() (Fortran)
>
>  -7.105427e-15
>
> s1-s2[] ## R versus manual sum (Fortran
>
>  -7.105427e-15
>
> s1-s3 ## R Versus manual sum in C
>
>  -7.105427e-15
>
> s2[]-s2[] ## manual sum versus compiler sum() (Fortran)
>
>  0
>
> s3-s2[] ## Fortran versus C
>
>  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-tierney at uiowa.edu
Iowa City, IA 52242                 WWW:  http://www.stat.uiowa.edu

```