モンテカルロ法で円周率を求める

モンテカルロ法は乱数を使って計算する手法です。この記事ではモンテカルロ法で円周率を計算します。動作環境はVisual Studio Community 2022のバージョン17.7.4です。使用言語はC++20です。

理論

以下の図のような辺の長さが1の正方形と半径が1、中心角が90°の扇形について考えます。

出典: https://note.com/shimakaze_soft/n/n9547f5c0bae0

この正方形の中のランダムな位置に大量の点を打ちます。点を打っていくと扇形の中に含まれる点と含まれない点が出てきます。モンテカルロ法では、扇形の中の点の総数と正方形の中に打った点の総数の比が、扇形の面積と正方形の面積の比に近似できると考えます。この考えによって以下の式が成り立ちます。

\frac{扇形の中に打った点の合計}{正方形の中に打った点の合計} = \frac{\pi ×1×1×\frac{90}{360}}{1×1}

この式から円周率は以下の式で求められます。

\pi = \frac{扇形の中に打った点の合計}{正方形の中に打った点の合計}×4

実装

以下のコードでモンテカルロ法を用いた円周率の計算が出来ます。

MonteCarlo.cpp

#include <iostream>
#include <random>
#include <numbers>

int main()
{
	int points = 1000000; //点の合計
	int circlePoints = 0; //扇形の中に入った点の数
	std::random_device rng; //乱数インスタンスの生成
	std::uniform_real_distribution<> dist(0, 1); //0以上1未満の実数の乱数を発生させる

	for (int i = 0; i < points; i++)
	{
		double x = dist(rng); //x座標をランダムで指定
		double y = dist(rng); //y座標をランダムで指定

		if (std::pow(x, 2) + std::pow(y, 2) <= 1)
		{
			circlePoints++; //扇形の中に点があれば1を足す
		}
	}

	double pi = (double)circlePoints / points * 4; //点の数の比から円周率を求める

	std::cout << pi << std::endl; //モンテカルロ法で求めた円周率を表示
    //3.14446
	std::cout << std::numbers::pi << std::endl; //標準ライブラリで定義された円周率を表示
    //3.14159
}

扇形の中に点が入っているかどうかは、原点(0, 0)と点の間の距離が扇形の半径1以下であるかどうかで判定しています。
モンテカルロ法による円周率の計算結果は3.14446になっていますが、乱数計算なので実行する度に結果は変わります。しかし複数回実行してもほとんどの場合で約3.14という結果が出ます。このコードでは打った点の合計を1000000にしましたが、点の数を増やせば精度はさらに良くなると考えられます。
C++では標準ライブラリのnumbersに円周率が定数として定義されているので、モンテカルロ法の計算結果と比較してみてください。

参考

https://note.com/shimakaze_soft/n/n9547f5c0bae0
https://anko.education/joho/programming/monte_carlo_pi