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

はじめに

モンテカルロ法による各言語の速度比較という記事を読みました。
D言語のことも紹介したいと思ったのが、この記事を書くきっかけとなります。
※こちらの記事では、各言語との速度比較は行いません。

コンパイル環境

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

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

ソースコード１

元記事を読んで、乱数生成方法についての明確なルールはなさそうですが、D言語の標準ライブラリを使った場合のソースコードです。
標準ライブラリにあるRandom関数は、メルセンヌ・ツイスタを使っています。

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を使うことで最適化が効いています。

コマンドプロンプト

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使った場合のソースコードです。
ソースコード１との差異は、乱数生成部分だけです。元記事にもあるように、乱数生成処理を変えると処理速度が大きく改善しました。

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 );
}

dub.json

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

コマンドプロンプト

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

ソースコード３

モンテカルロ法では、乱数を使うのが前提となっています。
処理全体の中で、乱数生成処理が占めるコスト（処理時間）が高いと思われるため、乱数を使わない場合をベンチマークとして考えてみました。
乱数を使わない場合の方式として思いついたのが、

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

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

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) );
}

コマンドプロンプト

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%