並列計算の効果検証用プログラム

並列化の試行やJITの効果テストなどのために，Forをたくさん回す計算例が必要になる。 そのような例としてモンテカルロ法を置いておく。

ToDo 並列計算のテスト用としては，モンテカルロ法のような単純に並列化される(データ共有が無い)ものだけでなく，適時データ共有が必要となる例を別途置いておく。

モンテカルロ法による円周率計算

• モンテカルロ法の例題でよく見る計算。[0,1]の範囲の乱数を生成して，その点が半径1の円の中に入るか否か？(点が円の中か否かは，点の距離dがsqrt(x²+y²)で計算できるから，円周率piを知らなくても判定できる)をランダムに試行して，N回サンプリングすると，円の中に入る比率が1:pi/4になることから推定する(円の面積がpirrなのは知っている)

# monte_carlo.py
import random

def calc_pi(iter):
    cnt = 0
    for _ in range(iter):
        x, y = random.random(), random.random()
        if x*x + y*y < 1:
            cnt += 1
    return (cnt*4.0)/iter

• 全く同じ書き方でC言語でも。

#include <stdio.h>
#include <stdlib.h>

double calc_pi(int iter)
{
  int cnt=0;
  double x, y;
  
  for (int i=0; i<iter; i++) {
    x = (double)rand()/RAND_MAX;
    y = (double)rand()/RAND_MAX;
    if ((x*x + y*y) < 1)
      cnt++;
  }
  return cnt*4.0/iter;
}

int main(int argc, char** argv) {
  int iter;
  double pi;
  
  iter = atoi(argv[1]);
  pi = calc_pi(iter);
  printf("pi:%f", pi);

  return 0;
}

• ちなみに，生のPythonは20倍くらい遅い。

# Python
$ time python3 monte_carlo.py 10000000 # 10M
real 0m3.838s                          # 0.4us/iter

# C(Clang, -O3)
$ time ./a.out 10000000                # 10M
real 0m0.214s                          # 0.02us/iter