Rcppを使う：その２

C言語による乱数生成
ここにメルセンヌツイスターのソースコードがあるので利用してみる。

あらかじめMT.hをincludeフォルダにおいておく必要がある。

以下のファイルをpi4.cppとして保存

#include <Rcpp.h>
#include <stdio.h>
#include <MT.h>
using namespace Rcpp;
// [[Rcpp::export]]
double pi4Cpp(){
long circle_in = 0;
init_genrand(10);
for (long i = 0; i < 100000000; ++i) {
double  l = pow(genrand_real3(),2) + pow(genrand_real3(),2);
if (l < 1.0){
      circle_in++;   
}
}
//return(circle_in);
return(4.0 * circle_in / 100000000);
}

そして、以下のRコードを実行する。

library(Rcpp)
library(rbenchmark)

sourceCpp("pi4.cpp")

cppFunction('double pi3Cpp(){
      long circle_in = 0;
      Rcpp::NumericVector x = runif(200000000);
      for (long i = 0; i < 100000000; ++i) {
            double  l = x[i]*x[i] + x[i+100000000]*x[i+100000000];
            if (l < 1.0){
                  circle_in++;   
            }
      }
      return(4.0 * circle_in / 100000000);
}')

f0 <- function(){
      n = 100000000
      x <- runif(n = n)
      y <- runif(n = n)
      print(4*sum(x^2 + y^2 <= 1)/n)
}

benchmark(f0(),pi3Cpp(),pi4Cpp(),replications = 1)

> benchmark(f0(),pi3Cpp(),pi4Cpp(),replications = 1)
[1] 3.141453
[1] 3.141787
      test replications elapsed relative user.self sys.self user.child sys.child
1     f0()            1   11.09    8.338     10.17     0.83         NA        NA
2 pi3Cpp()            1    2.87    2.158      2.59     0.26         NA        NA
3 pi4Cpp()            1    1.33    1.000      1.32     0.00         NA        NA
> 

１１秒だったのが、Rcppを使用して毎ループでメルセンヌツイスター乱数を
呼び出しても、１．３３秒まで短くなった。