Skip to content

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

4 messages · Peng Yu, William Dunlap, Dirk Eddelbuettel +1 more

Original
Peng Yu
# 
Hi,

The following shows the arithmetic in C++ is slightly different from
in R. Is there a way to get the exact arithmetic as in R?

~/dvcs_src/rexample/Rexample/cran/base/sum$ Rscript main_Rcpp.R
+ Rcpp::NumericVector vec(vx);
+ double sum = 0;
+ for (int i=0; i<vec.size(); i++) {
+   sum += vec[i];
+ }
+ return Rcpp::wrap(sum);
+ '
+   signature(vx='numeric')
+   , src
+   , plugin='Rcpp'
+   )

Attaching package: ?Rcpp?

The following object is masked from ?package:inline?:

    registerPlugin
[1] 0.9999999999999998889777
[1] 1
--
Regards,
Peng
William Dunlap
# 
Make your sum quad precision ('long double' in g++).

% cat s.cpp
#include <iostream>
int main(int argc, char *argv[])
{
    double dSum(0.0);
    long double ldSum(0.0);
    for(int i=0 ; i < 10 ; i++) {
        dSum += 0.1;
        ldSum += 0.1;
    }
    double dLdSum(ldSum);
    std::cout << "diff=" << dSum - dLdSum << "\n";
    std::cout << "  sizeof (double)=" << sizeof (double) << "\n";
    std::cout << "  sizeof (long double)=" << sizeof (long double) << "\n";
    return 0;
}
% g++ s.cpp
% ./a.out
diff=-1.11022e-16
  sizeof (double)=8
  sizeof (long double)=16


Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com
Dirk Eddelbuettel
#
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
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:
[1] 0.8999999999999999111822
[1] 0.8999999999999999111822
[1] 0.9000000000000000222045
#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);
}
$Ldouble
[1] 0.9000000000000000222045

$Double
[1] 0.8999999999999999111822

$auto
[1] 0.8999999999999999111822

$`Rcpp::sum`
[1] 0.8999999999999999111822

$`R sum`
[1] 0.9000000000000000222045

chris