Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
C vs Rcpp for monte carlo pi approximation
Rcpp_code <- "
#include <Rcpp.h>
// [[Rcpp::export]]
double mcsim_rcpp(const int n){
int i, r = 0;
double u, v;
Rcpp::RNGScope scope;
for (i=0; i<n; i++){
u = R::runif(0, 1);
v = R::runif(0, 1);
if (u*u + v*v <= 1)
r++;
}
return (double) 4.*r/n;
}
"
C_code <- "
#include <R.h>
#include <Rinternals.h>
// [[Rcpp::export]]
SEXP mcsim_c(SEXP n_){
SEXP ret;
int i, r = 0;
const int n = INTEGER(n_)[0];
double u, v;
GetRNGstate();
for (i=0; i<n; i++){
u = unif_rand();
v = unif_rand();
if (u*u + v*v <= 1)
r++;
}
PutRNGstate();
PROTECT(ret = allocVector(REALSXP, 1));
REAL(ret)[0] = (double) 4.*r/n;
UNPROTECT(1);
return ret;
}
"
library(Rcpp)
sourceCpp(code=Rcpp_code)
sourceCpp(code=C_code)
mcsim_c(100L)
mcsim_rcpp(100L)
mcsim_rcpp(100)
mcsim_c(100) ### uh-oh
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment