自然数2つで円周率が求まる！？

# はじめに
皆さんこんにちは！
つい最近、このような動画を視聴しました。
<iframe width="560" height="315" src="https://www.youtube.com/embed/Yl6Rb-ri158?si=QyhP8WKGd7bLY4k3" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture; web-share" allowfullscreen></iframe>
おすすめに流れてきて、しっかり最後まで見させてもらいました！

てことで、この動画で紹介されている方法を使って円周率を求めていきたいと思います！
ただ、普通に求めるのも面白くないので、いくつかのパラメータを動かして結果を比べていきたいと思います！

# 円周率の求め方
※動画でとても丁寧に解説されているので、ここでは証明はせずに要点のみを抽出します。
1. xを大きな自然数とする
1. ある自然数a,b(1≦a,b≦x)を2つ選ぶ
1. a,bが互いに素であるかを調べる
1. これをn回繰り返し、以下の式に当てはめる
```math
\displaylines{
\frac{互いに素であった回数}{n(繰り返した回数)}\fallingdotseq\frac{6}{\pi^2}
}
```
こちらの内容が動画内[11:11](https://youtu.be/Yl6Rb-ri158?t=671)から解説されています

# プログラムの構想
ここでは、c++を使用します
### ランダム関数
ある自然数a,bを取得するためには、random関数を作成します
```c++
#include <random>
using namespace std;

// random関数.
random_device rd;
mt19937 gen(rd());
unsigned long long random(unsigned long long low, unsigned long long high)
{
    uniform_int_distribution<> dist(low, high);
    return dist(gen);
}
```
生成方法は以下のとおりです
1. デバイスのノイズからseedを取得
1. mt19937(メルセンヌ・ツイスター法)での乱数生成

これで1つ完成です
### GCD(最大公約数)関数
互いに素であることを調べるために、動画と同様にユークリッドの互除法を活用して、
最大公約数を求めて1であるかを調べていこうと思います
```c++
using namespace std;

// 最大公約数.
long long GCD(long long a, long long b)
{
    if (b == 0)
        return a;
    else
        return GCD(b, a % b);
}
```
whileループなどの他のやり方もありますが、ここでは再帰関数で求めていきます。

### 仕上げ
最後にメイン関数を書いていきます！
`forループでn回繰り返し、a,bが互いに素であるかを調べ、その個数/nをする`
では、先程の数式をπについて解くために少し変形させます
先程の数式は日本語が多く含まれていたので、全て文字に置き換えます

c: `a,bが互いに素であった回数`
n: `forループ回数`
```math
\displaylines{
\frac{c}{n}\fallingdotseq\frac{6}{\pi^2}
\\\Rightarrow\frac{c}{6n}\fallingdotseq\frac{1}{\pi^2}
\\\Rightarrow\frac{6n}{c}\fallingdotseq{\pi^2}
\\\Rightarrow\sqrt{\frac{6n}{c}}\fallingdotseq{\pi}
}
```
1,2行目: 両辺を1/6倍
2,3行目: 逆数を取る(-1乗)
3,4行目: 正の平方根を取る

では、コーディングしていきます！
```c++
#include <iostream>
#include <iomanip>
using namespace std;

int main()
{
    // 試行回数.
    long long n = 1000000;

    // 乱数最大値.
    long long x = 100000000;

    // メイン.
    double counter = 0;
    for (int i = 0; i < n; i++)
    {
        if (GCD(random(1, x), random(1, x)) == 1)
            counter++;
    }

    // output.
    cout << fixed << setprecision(16);
    cout << sqrt(6 * n / counter) << endl;
    return 0;
}
```
### 完成
そうして、なんやかんやで完成したコードがこちらです！
```c++
#include <random>
#include <iostream>
#include <iomanip>
using namespace std;

// random関数.
random_device rd;
mt19937 gen(rd());
unsigned long long random(unsigned long long low, unsigned long long high)
{
    uniform_int_distribution<> dist(low, high);
    return dist(gen);
}

// 最大公約数.
long long GCD(long long a, long long b)
{
    if (b == 0)
        return a;
    else
        return GCD(b, a % b);
}

int main()
{
    // 試行回数.
    long long n = 10000000;

    // 乱数最大値.
    long long x = 100000000;

    // メイン.
    double counter = 0;
    for (int i = 0; i < n; i++)
    {
        if (GCD(random(1, x), random(1, x)) == 1)
            counter++;
    }

    // output.
    cout << fixed << setprecision(16);
    cout << sqrt(6 * n / counter) << endl;
    return 0;
}
```
近似とはいえ、こんなに簡単に円周率が求まってしまうんですね～
では、早速実験していきましょう！

# 実験
ここでは、10回の平均を取っていきます！
### n=100万, x=1億
まずは、動画の通りのパラメータでやっていきます
結果はこちら！
|回数|結果|
|---|---|
|1|3.1437112362736928|
|2|3.1404695417881592|
|3|3.1416445951819094|
|4|3.1427304313757491|
|5|3.1426605937551630|
|6|3.1420684544114663|
|7|3.1424381790289990|
|8|3.1427097382620417|
|9|3.1406605595944024|
|10|3.1416704353462390|

平均は、**3.1420763765017821**になりました！
シミュレーションにしては意外と精度良く、3.14まで合いました！

時間に関しては約1.7秒です。
|real|user|sys|
|---|---|---|
|0m1.735s|0m1.716s|0m0.010s|

### n=1000万, x=1億
次に試行回数を増やして誤差が減るかの実験です！
結果はこちら！
|回数|結果|
|---|---|
|1|3.1420245098138810|
|2|3.1417843980774904|
|3|3.1413578121695238|
|4|3.1415024856743785|
|5|3.1416425279963072|
|6|3.1408614240144335|
|7|3.1411695082100399|
|8|3.1416872317949758|
|9|3.1421778069125055|
|10|3.1412196159707375|

平均は、**3.1415427320634275**となりました！
3.1415まで一致していますね！

時間に関しては約17秒です。
|real|user|sys|
|---|---|---|
|0m17.027s|0m17.010s|0m0.010s|

つまり、このままもう1桁増やしても170秒程度ですね
次で、nを増やすのを終了します

### n=1億, x=1億
いきなりですが、結果はこちら！
|回数|結果|
|---|---|
|1|3.1417626640916740|
|2|3.1415038291648765|
|3|3.1418313822422599|
|4|3.1416302283263460|
|5|3.1416480835668654|
|6|3.1414748669920374|
|7|3.1415890151491426|
|8|3.1415277280819738|
|9|3.1417297409342684|
|10|3.1415731504462032|

平均は、**3.1416270688995649**となりました！
試行回数を増やしていくのには限界がありましたね～
最大でも3.1415までしか合わないみたいです

では、次に乱数の数を大きくしてみます！
### n=100万, x=10億
|回数|結果|
|---|---|
|1|3.1423037186130411|
|2|3.1428571856171050|
|3|3.1422390803257998|
|4|3.1408128825513790|
|5|3.1417634652175086|
|6|3.1417841396445696|
|7|3.1404385692118013|
|8|3.1408077186983947|
|9|3.1413009816180719|
|10|3.1412983984844733|

平均は、**3.1415606139982146**です！
やっぱり乱数だと良い時と悪い時がありますね...

### n=100万, x=100億
|回数|結果|
|---|---|
|1|3.1405573023979874|
|2|3.1398889584298484|
|3|3.1420115852831327|
|4|3.1429580826080277|
|5|3.1409419871542825|
|6|3.1429477337514795|
|7|3.1409600630691923|
|8|3.1433695324688835|
|9|3.1422494221837050|
|10|3.1396877646815917|

平均は、**3.1415572432028132**です！
意外と3.1415までは取れるみたい...
まあ、xに関しては対して値の変動はありませんでしたね
# 結論
:::note
シミュレーションで円周率を求めるのには限界がある！(この場合だと3.1415)
:::
# 最後に
久しぶりのqiitaの投稿ですが、面白い動画を見て、実験したいと思いやってみました！
また何かあればメモ程度にqiitaに投稿します。
最後まで見てくださりありがとうございました！