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

はじめに

皆さんこんにちは！
つい最近、このような動画を視聴しました。

おすすめに流れてきて、しっかり最後まで見させてもらいました！

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

円周率の求め方

※動画でとても丁寧に解説されているので、ここでは証明はせずに要点のみを抽出します。

1. xを大きな自然数とする
2. ある自然数a,b(1≦a,b≦x)を2つ選ぶ
3. a,bが互いに素であるかを調べる
4. これをn回繰り返し、以下の式に当てはめる

\displaylines{
\frac{互いに素であった回数}{n(繰り返した回数)}\fallingdotseq\frac{6}{\pi^2}
}

こちらの内容が動画内11:11から解説されています

プログラムの構想

ここでは、c++を使用します

ランダム関数

ある自然数a,bを取得するためには、random関数を作成します

#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を取得
2. mt19937(メルセンヌ・ツイスター法)での乱数生成

これで1つ完成です

GCD(最大公約数)関数

互いに素であることを調べるために、動画と同様にユークリッドの互除法を活用して、
最大公約数を求めて1であるかを調べていこうと思います

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ループ回数

\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行目: 正の平方根を取る

では、コーディングしていきます！

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

完成

そうして、なんやかんやで完成したコードがこちらです！

#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に関しては対して値の変動はありませんでしたね

結論

シミュレーションで円周率を求めるのには限界がある！(この場合だと3.1415)

最後に

久しぶりのqiitaの投稿ですが、面白い動画を見て、実験したいと思いやってみました！
また何かあればメモ程度にqiitaに投稿します。
最後まで見てくださりありがとうございました！