Rcppでループをまわす話

inspired by http://bicycle1885.hatenablog.com/entry/2015/01/04/175916

速度がどうというより、愚直なループをみたらRcppをさくっと適用できるようになりたい。

見本(Julia)はこちら。

circle_in = 0.0
for i in 1:100000000
  l = (rand()^2 + rand()^2) ^ 0.5
  if l <= 1
    circle_in = circle_in + 1
  end
end
println(4 * circle_in / 100000000)

Rで愚直に書くとこんな感じ。

calcPiR <- function(){
  circle_in = 0
  for (i in 1:100000000){
    l = sqrt(runif(1)^2 + (runif(1))^2)
    if (l <= 1){
      circle_in = circle_in + 1      
    }
  }
  4 * circle_in / 100000000 
}

しかしこれだといつまで経っても終わらないので少しはましなコードにしてそちらを実験に使う。

ましなコードはこちらからもってきた。

http://gallery.rcpp.org/articles/simulating-pi/

ということでRと愚直なC++のループとRのコードっぽいC++で比較した。

library(Rcpp)
library(microbenchmark)

calcPiR2 <- function() {
    x <- runif(100000000)
    y <- runif(100000000)
    d <- sqrt(x^2 + y^2)
    return(4 * sum(d < 1.0) / 100000000)
}

cppFunction('double calcPiCpp() {
            double circle_in = 0.0;
              for (int i = 0; i < 100000000; i++) {
                double l = sqrt(pow(rand()/2147483648., 2)+pow(rand()/2147483648.,2));
                if (l <= 1){
                  circle_in = circle_in + 1.0;
                  }
              }
              return 4*circle_in/100000000;
            }')

cppFunction('double calcPiCpp2() {
            NumericVector x = runif(100000000);
            NumericVector y = runif(100000000);
            NumericVector d = sqrt(x*x + y*y);
            return 4.0 * sum(d < 1.0) / 100000000;
            }')

microbenchmark(calcPiR2(), calcPiCpp(), calcPiCpp2(), times=3L)

## Unit: seconds
##          expr       min        lq      mean    median        uq       max
##    calcPiR2() 12.239480 14.391310 15.603079 16.543140 17.284878 18.026615
##   calcPiCpp()  2.522930  2.523942  2.575725  2.524954  2.602123  2.679291
##  calcPiCpp2()  6.314802  6.644208  6.842224  6.973614  7.105935  7.238255
##  neval
##      3
##      3
##      3