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: [email protected]
[[email protected]] on behalf ofDirk Eddelbuettel
[[email protected]]
Sent: 09 July 2013 07:57mbe
To: Peng Yu
Cc: [email protected]
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 | [email protected] | http://dirk.eddelbuettel.com
_______________________________________________
Rcpp-devel mailing list
[email protected]
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
_______________________________________________
Rcpp-devel mailing list
[email protected]
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
