モンテカルロ法で円周率を計算するベンチマーク
モンテカルロ法で円周率を計算するベンチマークで、Juliaと比較してgccが悪い言われているのですが、本当なのか検証してみましょう。
CとJuliaでこの分岐が単純なマイクロベンチマークを実行すると、乱数生成の速度が結果を大きく左右しえます。後方互換性の問題で、Cの標準の乱数生成関数は線形合同法になっています。Juliaの方は標準がdouble precision SIMD oriented Fast Mersenne Twister pseudorandom number generator (dSFMT)で、明示的にMersenneTwisterを使った場合もdSFMTを用いているようです。Cの方もdSFMTを用いてベンチマークを取ることにします。
Cのソースコード
dSFMTがよくチューニングされているので、私には工夫の余地は無く、オーソドックスに書いています。なお、生成した実行ファイルをreadelfで調べたところ、dSFMTがヘッダーファイルでinline static宣言している関数はインライン展開されていました。10回計測して、一番短い経過時間を採用しています。
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <limits.h>
#include "dSFMT.h"
#include <float.h>
#include <math.h>/*
dSFMT-src-2.2.3.tgz を展開したディレクトリで、ソースコードのファイル名をpi_mc_dsfmt.cとして、clang -o pi_mc pi_mc_dsfmt.c dSFMT-src-2.2.3/dSFMT.c -IdSFMT-src-2.2.3/ -DDSFMT_MEXP=521 -msse2 -DHAVE_SSE2 -O3 -lm
とコンパイルしました。gccとclangはオプションが同じなので、gccとの比較は容易に行えます。
*/double pi_mc(long n_iter) {
double x, y;
dsfmt_t rng;dsfmt_init_gen_rand(&rng, time(NULL));
int n_in = 0;
for (long i = 1; i <= n_iter; i++) {x = dsfmt_genrand_open_open(&rng);
y = dsfmt_genrand_open_open(&rng);if (x * x + y * y < 1.0)
n_in++;
}return (double)4.0*n_in/n_iter;
}int main(){
struct timespec t1, t2;
double dt, sum_dt, min_dt;
double r = 0;
size_t n;
int j, nos = 10;printf("log10(n)\telapsed time (seconds)\n");
for(n=1; n<=1000000*1000; n*=10){
sum_dt = 0;
min_dt = DBL_MAX;
for(j=0; j<nos; j++){
clock_gettime(CLOCK_THREAD_CPUTIME_ID, &t1); // clock()を使って計測しても、ほぼ同じ数字になりますr = pi_mc(n);
clock_gettime(CLOCK_THREAD_CPUTIME_ID, &t2);
// fprintf(stderr, "estimated pi: %f\n", r);
dt = (double)(t2.tv_sec - t1.tv_sec) + (double)(t2.tv_nsec - t1.tv_nsec)/1000000000;
sum_dt += dt;
if(min_dt > dt)
min_dt = dt;
}printf("%.0f\t%.09f\n", log10(n), min_dt);
}return 0;
}
Juliaのスクリプト
Cのとだいたい同じですが、グローバル変数へのアクセスが無いように一番外側ループ部分も関数内にしまいました。
using Random
# using Random: default_rngfunction pi_mc(n_iter)
n_in = 0
# rng = default_rng()
rng = MersenneTwister(time_ns())
for i in 1:n_iter
x, y = rand(rng), rand(rng)
n_in += ifelse(x^2 + y^2 ≤ 1, 1, 0)
end
4n_in/n_iter
endfunction loop(s, e)
n = s
while(n <= e)
sum_dt = 0;
min_dt = 3600;
for i in 1:10
t1 = time_ns()
r = pi_mc(n)
t2 = time_ns()
dt = (t2 - t1)/1000000000
sum_dt = sum_dt + dt
if(min_dt > dt)
min_dt = dt
end
# println("estimated pi: ", r)
end
println(log10(n), "\t", min_dt)
n = n*10
end
endprintln("log10(n)", "\t", "elapsed time (seconds)")
loop(1, 10^9)