[Rcpp-devel] Getting the exact arithmetic as in R when Rcpp is used?

[Wray, Christopher]

This has already been answered by Dirk and Bill. Although the OPs example can, of course, be reproduced in R, sans Rcpp. I use rep(0.1,9) in place of rep(0.1,10).

We can use Rcpp to see various differences explicitly:


> options(digits=22)
> vx=rep(0.1,9)
> vx[1]+vx[2]+vx[3]+vx[4]+vx[5]+vx[6]+vx[7]+vx[8]+vx[9]
[1] 0.8999999999999999111822
> sum(vx[1],vx[2],vx[3],vx[4],vx[5],vx[6],vx[7],vx[8],vx[9])
[1] 0.8999999999999999111822
> sum(vx)
[1] 0.9000000000000000222045
> Sys.setenv("PKG_CXXFLAGS"="-std=c++0x")

#include <Rcpp.h>

// [[Rcpp::export]]
Rcpp::List sumK(Rcpp::NumericVector K, Rcpp::Function f) {
  
  long double sumv=0.0;
    for(int k=0; k<K.size(); ++k)  sumv+=K(k);

  double sumk=0.0;
    for(int k=0; k<K.size(); ++k)  sumk+=K(k);

  auto suma=0.0;
    for(int k=0; k<K.size(); ++k)  suma+=0.1;
  
  long double sumr=Rcpp::sum(K);
  
  Rcpp::NumericVector sumf=f(K);
  
  return Rcpp::List::create(Rcpp::Named( "Ldouble")=Rcpp::wrap(sumv), Rcpp::Named("Double")=Rcpp::wrap(sumk), Rcpp::Named( "auto")=Rcpp::wrap(suma), Rcpp::Named("Rcpp::sum")=Rcpp::wrap(sumr), Rcpp::Named("R sum")=sumf);
}

> sumK(rep(0.1,9),sum)
$Ldouble
[1] 0.9000000000000000222045

$Double
[1] 0.8999999999999999111822

$auto
[1] 0.8999999999999999111822

$`Rcpp::sum`
[1] 0.8999999999999999111822

$`R sum`
[1] 0.9000000000000000222045

chris
________________________________________ 
From: rcpp-devel-bounces at lists.r-forge.r-project.org [rcpp-devel-bounces at lists.r-forge.r-project.org] on behalf ofDirk Eddelbuettel [edd at debian.org]
Sent: 09 July 2013 07:57mbe
To: Peng Yu
Cc: rcpp-devel at lists.r-forge.r-project.org
Subject: Re: [Rcpp-devel] Getting the exact arithmetic as in R when Rcpp is     used?

On 8 July 2013 at 06:51, Peng Yu wrote:
| > vx=rep(.1, 10)
| > options(digits=22)
| > fun(vx)
| [1] 0.9999999999999998889777

You are printing a data type that has around 16 decimals precision with 22.
That is bound to show random stuff at the end.

Otherwise, Bill is quite right that one can (and R Core does) use 'long
double' in the internal parts of loops which could otherwise accumulate
numerical error.  There are many texts on numerical computing which cover
this. A short and classic paper by David Goldbergs is available in many
places under 'what ever computer scientist should know about floating-point
arithmetic'.

Dirk

--
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel at lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel