今朝はとても夜のような朝でしたね。おはようございます。まがりかどです。
今回は円周率をC++でモンテカルロ法という計算方法で求めていきたいと思います。
実はこのブログで前一回円周率を同じ方法で求めたんですけど、その時はPythonで可視化しながら求めてだいぶ計算速度が遅くなってしまったので今回はC++でもっと精度を上げていきたいなと思います。
なので今回はモンテカルロ法の手法については触れませんので、もし見ておきたかったらこれをどうぞ。
ということでもうコードを貼っていきますよ。
円周率を求めるC++コード
#include<bits/stdc++.h>
using namespace std;
using lint = long long;
int main(){
cout << "How many times?" << endl;
lint loop;
cin >> loop;
lint cnt = 0;
lint range = 10000;
random_device rnd;
mt19937_64 mt(rnd());
uniform_int_distribution<> randrange(0, range);
for(lint i = 0; i < loop; i++){
lint x = randrange(mt);
lint y = randrange(mt);
if(x * x + y * y <= range * range) cnt++;
}
double ans = ((double)cnt / (double)loop) * 4;
cout << setprecision(5) << ans << endl;
}わたしはlong longをlintにする派なので気を付けてくださいね。
たぶん ll にするほうがマジョリティーなんでしょうけど僕は嫌だ。
なんか整数感がないじゃないですよね。
というわけでlintを使っているのでご了承ください。
まあそんなことは置いておいて上から解説です。
シミュレーションだ
random_device rnd;
mt19937_64 mt(rnd());
uniform_int_distribution<> randrange(0, range);
for(lint i = 0; i < loop; i++){
lint x = randrange(mt);
lint y = randrange(mt);
if(x * x + y * y <= range * range) cnt++;
}まずrandom_device という非決定的なハードウェアの状態から乱数を発生してくれるやつがあるんですけど、それをループに使うと計算速度が遅くなってしまうので別の決定的な乱数生成アルゴリズムのシードにしてあげよう。
その別の決定的な乱数生成アルゴリズムというのはメルセンヌツイスタというものです。
最初この名前見たときメルセンヌツイスタっていう人が考えたのかなとか思ってたんですけど、これ日本人が考えたらしいですよ。
つまり二行目でメルセンヌツイスタのシードを設定しています。
そして三行目では一様分布器を宣言してこれの引数にさっきのメルセンヌツイスタを入れたら一様分布乱数ができるらしいです。
一様分布とか正直よくわからないんですが要はどれもこれも等しい確率で出るらしいです。
で、それを使って打つ点のxとyを求めて
$$x^2 + y^2 <= range^2$$
かどうかで点が円の中に入っているかどうかを判定します。
ここでわざわざ平方根を取って距離を求めていると計算量がちょっと増えるのでいやですね。
なので距離の二乗で比較します。
そして最後
double ans = ((double)cnt / (double)loop) * 4;
cout << setprecision(5) << ans << endl;ちゃんと(double)でキャストしておかないと多分０になるので注意ですね。
setprecision(5)っていうのは小数型の変数を５桁まで表示しますよっていうやつですね。
C++で円周率を求めた結果
これで実行すると
How many times?
100000000
3.1416結構近いですね。いい感じです。
Pythonで同じように1億回シミュレーションしようとすると全然終わらないのでやっぱりC++ありがたいですね。
ありがとうございました
