モンテカルロ法で円周率を求めてみた

# はじめに
[モンテカルロ法による各言語の速度比較](https://qiita.com/kurooolong/items/2c95ac40a0c051c5a6ae)という記事を読みました。
D言語のことも紹介したいと思ったのが、この記事を書くきっかけとなります。
※こちらの記事では、各言語との速度比較は行いません。

# コンパイル環境
* OS : Windows 10 64ビット版 (バージョン 1909)
* D言語コンパイラ : LDC - the LLVM D compiler (1.20.1)

速度比較という観点で、コンパイラは`LLVM`系の`LDC`を使いました。

# ソースコード１
元記事を読んで、乱数生成方法についての明確なルールはなさそうですが、D言語の標準ライブラリを使った場合のソースコードです。
標準ライブラリにある[Random](https://dlang.org/phobos/std_random.html#Random)関数は、[メルセンヌ・ツイスタ](https://ja.wikipedia.org/wiki/%E3%83%A1%E3%83%AB%E3%82%BB%E3%83%B3%E3%83%8C%E3%83%BB%E3%83%84%E3%82%A4%E3%82%B9%E3%82%BF)を使っています。

```d:monte1.d
// ldc2 -m64 -release -O monte1.d
import std;

void main()
{
	double    ave   = 0.0;
	immutable times = 100;
	for ( int i = 1; i <= times; i++ ){
		writefln("times : %03d", i);
		ave += montecarlo();
	}
	writefln("AVE   : %.8f\n", ave / times);
}

double montecarlo()
{
	immutable times = 10_000_000;
	long time_start = Clock.currStdTime();					// 戻り値の単位hnsecs
	auto rnd = Random(cast(uint)time_start);
	long cnt =
		generate!(() => [uniform01(rnd), uniform01(rnd)])	// 乱数を2つ生成
		.takeExactly(times)									// times回繰り返す
		.map!(a => a[0] * a[0] + a[1] * a[1])				// それぞれ2乗して足して
		.count!("a <= b")(1.0);								// 1.0以下となる場合をカウント
	long time_end = Clock.currStdTime();
	double time_split = cast(double)(time_end - time_start) / 10_000_000;	// 時間を秒に変換
	writefln("pi    : %.8f", 4.0 * cnt / times);
	writefln("split : %.8f", time_split);
	
	return( time_split );
}
```
私の環境での実行例となります。`LDC`を使うことで最適化が効いています。

```text:コマンドプロンプト
D:\Dev> ldc2 -m64 -release -O monte1.d

D:\Dev> monte1
times : 001
pi    : 3.14075760
split : 0.09999950
times : 002
pi    : 3.14132680
split : 0.10001290

・・・省略・・・

times : 099
pi    : 3.14214240
split : 0.09903040
times : 100
pi    : 3.14100760
split : 0.10102790
AVE   : 0.09983340
```

# ソースコード２
D言語のパッケージマネージャ`DUB`で提供されている`Xoroshiro128Plus`使った場合のソースコードです。
ソースコード１との差異は、乱数生成部分だけです。元記事にもあるように、乱数生成処理を変えると処理速度が大きく改善しました。

```d:src\app.d
// dub --build=release --arch=x86_64 --compiler=ldc2
import std;
import mir.random.engine.xoshiro;

void main()
{
	double    ave   = 0.0;
	immutable times = 100;
	for ( int i = 1; i <= times; i++ ){
		writefln("times : %03d", i);
		ave += montecarlo();
	}
	writefln("AVE   : %.8f\n", ave / times);
}

double montecarlo()
{
	immutable times = 10_000_000;
	long time_start = Clock.currStdTime();					// 戻り値の単位hnsecs
	auto rnd = Xoroshiro128Plus(time_start);
	long cnt =
		generate!(() => [uniform01(rnd), uniform01(rnd)])	// 乱数を2つ生成
		.takeExactly(times)									// times回繰り返す
		.map!(a => a[0] * a[0] + a[1] * a[1])				// それぞれ2乗して足して
		.count!("a <= b")(1.0);								// 1.0以下となる場合をカウント
	long time_end = Clock.currStdTime();
	double time_split = cast(double)(time_end - time_start) / 10_000_000;	// 時間を秒に変換
	writefln("pi    : %.8f", 4.0 * cnt / times);
	writefln("split : %.8f", time_split);
	
	return( time_split );
}
```

```json:dub.json
{
	"name": "montecarlo",
	"dependencies": {
		"mir-random": "~>2.2.14",
	},
	"license": "BSL-1.0",
}
```

```text:コマンドプロンプト
D:\Dev\montecarlo> dub --build=release --arch=x86_64 --compiler=ldc2
times : 001
pi    : 3.14199840
split : 0.04703700
times : 002
pi    : 3.14141240
split : 0.04804310

・・・省略・・・

times : 099
pi    : 3.14194520
split : 0.04702730
times : 100
pi    : 3.14133000
split : 0.04799930
AVE   : 0.04728063
```

# ソースコード３
[モンテカルロ法](https://ja.wikipedia.org/wiki/%E3%83%A2%E3%83%B3%E3%83%86%E3%82%AB%E3%83%AB%E3%83%AD%E6%B3%95)では、乱数を使うのが前提となっています。
処理全体の中で、乱数生成処理が占めるコスト（処理時間）が高いと思われるため、乱数を使わない場合をベンチマークとして考えてみました。
乱数を使わない場合の方式として思いついたのが、

* `x`方向に`n`、`y`方向に`n`に分割する。
* 分割された各点を乱数で求める点の代わりに使う。

ソースコードにすると、以下の通りです。
ソースコード１と２では乱数で1000万の点を取得しているので、同程度の点（3200*3200）としました。

```d:monte2.d
// ldc2 -m64 -release -O monte2.d
// 円周率PI算出
// x方向にn、y方向にnに分割して、各点が円の内側か外側かを調べる
import std;

void main()
{
	int    n   = 3200;
	double pi  = montecarlo(n);
	double err = fabs(pi - PI) / PI;
	writefln("n, pi, err = %5d, %.8f, %.8f", n, pi, err);
}

double montecarlo(int n)
{
	long time_start = Clock.currStdTime();
	long cnt = 0;
	for ( int x = 0; x < n; x++ ){
		double dx = cast(double)x / n;
		for ( int y = 0; y < n; y++ ){
			double dy = cast(double)y / n;
			if ( dx * dx + dy * dy <= 1.0 ){
				cnt++;
			}
		}
	}
	long time_end = Clock.currStdTime();
	immutable times = 10_000_000;
	double time_split = cast(double)(time_end - time_start) / times;
	writefln("split : %.8f", time_split);
	return( 4.0 * cnt / (n * n) );
}
```

```text:コマンドプロンプト
D:\Dev> ldc2 -m64 -release -O monte2.d

D:\Dev> monte2
split : 0.01201970
n, pi, err =  3200, 3.14282539, 0.00039239
```

# まとめ
No.1とNo.2では円周率の平均を算出していないため、ソースコード１と２の実行結果100番目から引用しました。
処理時間はNo.3が圧倒的に速いですが、No.1やNo.2の方が誤差が少ないことは興味深いです。
規則正しく点を選ぶより、乱数を使って点を選ぶことに効果があるんだなと思いました。

|NO.|乱数生成|処理時間|円周率の計算結果|円周率との誤差|
|:--:|:--:|:--:|:--:|:--:|
|1|Random(メルセンヌ・ツイスタ)|0.09983340秒|3.14100760|0.018613%|
|2|Xoroshiro128Plus|0.04728063秒|3.14133000|0.008351%|
|3|乱数でない分割方式|0.01201970秒|3.14282539|0.039239%|