【数学】円周率の計算_モンテカルロ法
円周率をモンテカルロ法を使い、円の面積と正方形の面積の比にて求める方法は有名ですが、
円じゃなくとも、球、超球の体積からも求められるのではと思いついて、実際に計算してみました。
下記、「円周率の計算（2次元の場合）」、「円周率の計算（3次元の場合）」、「円周率の計算（4次元の場合）」の内容は、ほぼほぼ同じで、
5次元以上の場合も同様の考えで計算できると思われます。
モンテカルロ法（モンテカルロ・シミュレーション）は、疑似乱数を用いて、数学的問題を解決する手法。
サンプル数を大きくするほど、確率的に一定値に収束する性質が利用されている。
■円周率の計算（２次元の場合）
下図のような半径$r$の円（内接円）と正方形を考えます。
正方形、円の面積はそれぞれ、
\begin{aligned}
正方形&・・・ \left( 2r\right) ^{2}\\
円&・・・ \pi r^{2}
\end{aligned}
となる。
正方形内にランダムで点を打っていくことを考えて、
正方形に入る数を$N$、円に入る数を$n$とする。
\begin{aligned}
正方形に入る数&・・・ N\\
円に入る数&・・・ n
\end{aligned}
※当然、正方形に入る数$N$は、サンプル数と一致する。
正方形、円の面積と、正方形に入る数、円に入る数で比を取ると、
\begin{aligned}
\dfrac{円の面積}{正方形の面積}=\dfrac{円に入る数}{正方形に入る数}
\end{aligned}
となる。
式を展開していくと、
\begin{aligned}
\dfrac{円の面積}{正方形の面積}&=\dfrac{円に入る数}{正方形に入る数}\\
\dfrac{\pi r^{2}}{\left( 2r\right) ^{2}}&=\dfrac{n}{N}\\
\dfrac{\pi }{4}&=\dfrac{n}{N}
\end{aligned}
よって、
\begin{aligned}
\pi =4\cdot\dfrac{n}{N}
\end{aligned}
となり、円周率$\pi$を近似にて推定することができます。
■円周率の計算（３次元の場合）
下図のような半径$r$の球と立方体を考えます。
立方体、球の体積はそれぞれ、
\begin{aligned}
立方体&・・・ \left( 2r\right) ^{3}\\
球&・・・ \dfrac{4}{3}\pi r^{3}
\end{aligned}
となる。
立方体内にランダムで点を打っていくことを考えて、
立方体に入る数を$N$、球に入る数を$n$とする。
\begin{aligned}
立方体に入る数&・・・ N\\
球に入る数&・・・ n
\end{aligned}
※当然、立方体に入る数$N$は、サンプル数と一致する。
立方体、球の体積と、立方体に入る数、球に入る数で比を取ると、
\begin{aligned}
\dfrac{球の体積}{立方体の体積}=\dfrac{球に入る数}{立方体に入る数}
\end{aligned}
となる。
式を展開していくと、
\begin{aligned}
\dfrac{球の体積}{立方体の体積}&=\dfrac{球に入る数}{立方体に入る数}\\
\dfrac{\dfrac{4}{3}\pi r^{3}}{\left( 2r\right) ^{3}}&=\dfrac{n}{N}\\
\dfrac{\pi }{6}&=\dfrac{n}{N}
\end{aligned}
よって、
\begin{aligned}
\pi =6\cdot\dfrac{n}{N}
\end{aligned}
となり、円周率$\pi$を近似にて推定することができます。
■円周率の計算（４次元の場合）
半径$r$の４次元の超球と４次元の超立方体を考えます。
超立方体、超球の体積はそれぞれ、
\begin{aligned}
超立方体&・・・ \left( 2r\right) ^{4}\\
超球&・・・ \dfrac{\pi ^{2}}{2}r^{4}
\end{aligned}
となる。
超立方体内にランダムで点を打っていくことを考えて、
超立方体に入る数を$N$、超球に入る数を$n$とする。
\begin{aligned}
超立方体に入る数&・・・ N\\
超球に入る数&・・・ n
\end{aligned}
※当然、超立方体に入る数$N$は、サンプル数と一致する。
超立方体、超球の体積と、超立方体に入る数、超球に入る数で比を取ると、
\begin{aligned}
\dfrac{超球の体積}{超立方体の体積}=\dfrac{超球に入る数}{超立方体に入る数}
\end{aligned}
となる。
式を展開していくと、
\begin{aligned}
\dfrac{超球の体積}{超立方体の体積}&=\dfrac{超球に入る数}{超立方体に入る数}\\
\dfrac{\dfrac{\pi ^{2}}{2}r^{4}}{\left( 2r\right) ^{4}}&=\dfrac{n}{N}\\
\dfrac{\pi ^{2}}{32}&=\dfrac{n}{N}
\end{aligned}
よって、
\begin{aligned}
\pi ^{2}&=32\cdot\dfrac{n}{N}\\
\pi &=4\sqrt{2\cdot\dfrac{n}{N}}
\end{aligned}
となり、円周率$\pi$を近似にて推定することができます。
<プログラム>
#include
#include
#include
#include
#define MAX 10 // 1行の最大文字数(バイト数)
double calcPi2d(int num) {
int i;
int n;
double x1, x2;
double r;
double pi;
n = 0;
for (i = 0; i < num; i++) {
x1 = 2 * ((double)rand() / RAND_MAX) - 1;
x2 = 2 * ((double)rand() / RAND_MAX) - 1;
r = sqrt(x1 * x1 + x2 * x2);
if (r <= 1) {
n++;
}
}
pi = 4.0 * n / num;
return pi;
}
double calcPi3d(int num) {
int i;
int n;
double x1, x2, x3;
double r;
double pi;
n = 0;
for (i = 0; i < num; i++) {
x1 = 2 * ((double)rand() / RAND_MAX) - 1;
x2 = 2 * ((double)rand() / RAND_MAX) - 1;
x3 = 2 * ((double)rand() / RAND_MAX) - 1;
r = sqrt(x1 * x1 + x2 * x2 + x3 * x3);
if (r <= 1) {
n++;
}
}
pi = 6.0 * n / num;
return pi;
}
double calcPi4d(int num) {
int i;
int n;
double x1, x2, x3, x4;
double r;
double pi;
n = 0;
for (i = 0; i < num; i++) {
x1 = 2 * ((double)rand() / RAND_MAX) - 1;
x2 = 2 * ((double)rand() / RAND_MAX) - 1;
x3 = 2 * ((double)rand() / RAND_MAX) - 1;
x4 = 2 * ((double)rand() / RAND_MAX) - 1;
r = sqrt(x1 * x1 + x2 * x2 + x3 * x3 + x4 * x4);
if (r <= 1) {
n++;
}
}
pi = 4.0 * sqrt(2.0 * n / num);
return pi;
}
int main(void) {
char buf[MAX];
int num;
double pi1, pi2, pi3;
srand((unsigned)time(NULL));
fgets(buf, sizeof(buf), stdin);
num = atoi(buf);
pi1 = calcPi2d(num);
pi2 = calcPi3d(num);
pi3 = calcPi4d(num);
printf("円周率1：%20.18lf\n", pi1);
printf("円周率2：%20.18lf\n", pi2);
printf("円周率3：%20.18lf\n", pi3);
return 0;
}
