　宇宙の深淵からこんにちは！

\int d r {\hat{ψ}}^{†} (r) ぽ っ ぴ ん フ レ ン ズ \hat{ψ} (r)

です！今回は、Pythonで円周率 $π$ を計算する方法を13通り紹介したいと思います！

　以前、「もう円周率で悩まない！πの求め方10選」という(円周率界隈では有名?な)記事を拝読し、「π求めてぇ」という欲望が増幅したので、今回の記事を執筆するに至りました。そして今回、記事にある $π$ を求める方法のほとんど全てをPythonで実装してみました！

　いざやってみると、Pythonの基礎知識(関数、Numpy、for文、while文、if文、クラスなど)がふんだんに盛り込まれていたので、(筆者のような)プログラミング初心者が Python の練習としてやってみるのはすごくオススメです！

環境は以下の通りです。

バソコン：MacBook Pro(15-inch, 2018)
CPU：2.2 GHz 6コアIntel Core i7
メモリ：16 GB 2400 MHz DDR4
Pythonの環境：Anaconda, Jupyter Notebook(6.0.1), python 3.7.4

　それでは、円周率を計算する方法たちが織りなすハーレムをとくとご堪能ください！

Numpyの基本的な機能を使う

　PythonにはNumpyと呼ばれるライブラリがあり(Anacondaの場合は初めから用意されています)、ベクトルや行列その他数学の計算において非常に便利な機能が備わっています。以下では、Numpyの基本的な機能で円周率を計算する方法を紹介します！

1. numpy.pi

　1つ目は、numpy.pi を使う方法です。Numpyには、numpy.piという円周率の値を出してくれるやつが存在し、これを使うのが最も手っ取り早いです。もはや計算すらしていないですねw

import numpy as np #numpyをインポート

pi = np.pi #変数piにnp.piを代入
print(pi) #出力

出力は以下のようになります。

3.141592653589793

2. tan の逆関数

　あなたはある日突然「Numpy の機能のうち numpy.pi だけが使えなくなる呪い」にかけられてしまったとします。

　さぁ大変！1つ目の方法で円周率を出すことはできません！しかし、Numpyの他の機能はまだ使えるので、他の方法を探ってみましょう。ところで、有史以前から伝わる伝説の方法があります。それは、tan の逆関数を使う方法です。

\begin{matrix} (1) & \tan (\frac{π}{4}) = 1 ⟺ \arctan (1) = \frac{π}{4} \end{matrix}

であることを利用すれば、

\begin{matrix} (2) & π = 4 \arctan (1) \end{matrix}

と求めることができるのです。実際にやってみましょう！

import numpy as np

pi = 4 * np.arctan(1) #変数piに4*tanh(1)を代入
print(pi) #出力

出力は以下のようになります。

3.141592653589793

1つ目の方法とまったく同じ結果がでました！ちなみにこの方法は、人類初の「言語」であるFortranが生まれた先史時代から伝わるものです。まさに先人の知恵ですね！

面積から求める

「半径 $r$ の円の面積は $S = π r^{2}$ 」といえば、人類全員の深層心理に深く刻まれた普遍の常識ということで皆様ご存知かと思います！以下では、この誰もが知る円の面積を利用して円周率を求めてみましょう！

3. 正多角形近似

　円に内接する正 $n$ 角形の面積を増やしていくことで、値がどんどん円の面積に近づきます。半径 $r$ の円に内接する正 $n$ 角形の面積は

\begin{matrix} (3) & S_{n} = n \times \frac{r^{2}}{2} \sin (\frac{360^{\circ}}{n}) \to π r^{2} (n \to \infty) \end{matrix}

で求められます。半径 $r = 1$ の場合を想定して、実装してみましょう！

import numpy as np

#多角形近似
def RegularPolygon(N):
    theta = np.radians(360 / N) #３６０°をN等分し、度数法から弧度法に変換
    pi = N * np.sin(theta) / 2 #正N角形の面積を計算
    return pi

#出力
print("RegularPolygon(100)     = " + str(RegularPolygon(100)))
print("RegularPolygon(1000)    = " + str(RegularPolygon(1000)))
print("RegularPolygon(10000)   = " + str(RegularPolygon(10000)))
print("RegularPolygon(100000)  = " + str(RegularPolygon(100000)))
print("RegularPolygon(1000000) = " + str(RegularPolygon(1000000)))

　コードの中で用いているnp.radius()という関数は円周率の値に依存している関数で、円周率を求めるのに円周率に依存する関数を用いるのは個人的に気持ち悪いのですが、他にいい方法が思いつかなかったので、今回はこのままにしています(いい方法を知っていたら教えてください)。

　出力結果は以下のようになります。徐々に円周率に近づいていますね！

RegularPolygon(100)     = 3.1395259764656687
RegularPolygon(1000)    = 3.141571982779475
RegularPolygon(10000)   = 3.141592446881286
RegularPolygon(100000)  = 3.141592651522708
RegularPolygon(1000000) = 3.1415926535691225

4. 数値積分(長方形近似)

　円の面積を求める方法は多角形近似以外にも存在します。今回は積分(区分求積)の考えを使ってみましょう。

　関数 $y = f (x)$ 、直線 $x = a, x = b (a < b)$ 、および $x$ 軸で囲まれた部分の面積 $S$ を求めることを考えます。 $x$ 軸の区間 $[a, b]$ を $n$ 等分し、等分点をそれぞれ $x_{0}, x_{1}, x_{2}, \dots, x_{n} (ただし x_{0} = a, x_{n} = b)$ とします。また、隣り合う分点の間隔を $Δ x = (b - a) / n$ とします。このとき、区間 $[x_{i}, x_{i + 1}]$ の部分の面積 $S_{i}$ を $S_{i} = f (x_{i}) Δ x$ 、つまり底辺が $Δ x$ で高さが $f (x_{i})$ の長方形として近似します。あとは $S_{i}$ を全ての区間 $[x_{i}, x_{i + 1}] (i = 0, 2, \dots, n - 1)$ について求め、それらを合計することで近似的な面積 $S_{n}$ が求まります。そして、区間 $[a, b]$ を等分する数 $n$ を大きくしていけば、やがて真の面積 $S$ に近づいていきます。

\begin{matrix} (4) & S_{n} = \sum_{i = 0}^{n - 1} f (x_{i}) \frac{b - a}{n} \to S (n \to \infty) \end{matrix}

この方法を円の面積に応用しましょう。 $x - y$ 平面における単位円(原点中心、半径1)の方程式は

\begin{matrix} (5) & x^{2} + y^{2} = 1 ⟺ y = \pm \sqrt{1 - x^{2}} \end{matrix}

となりますが、プラスの方が上半分、マイナスの方が下半分を表します。今回は上半分のうちさらに $x >= 0$ の部分(四分円)の面積を求めます。求める面積の区間は $[a, b] = [0, 1]$ となるので、これを区分求積の公式(4)式に当てはめると

\begin{matrix} (6) & S_{n} = \sum_{i = 0}^{n - 1} \frac{\sqrt{1 - {(\frac{i}{n})}^{2}}}{n} \to \frac{π}{4} (n \to \infty) \end{matrix}

となります。

それでは実装していきましょう！

import numpy as np

#長方形近似の関数。区間[0,1]をN等分する。
def Rectangle(N):
    x = np.arange(N) / N #0~(N-1)/NまでのN要素配列
    y = np.sqrt(1 - x**2) #y=root(1-x^2)を計算
    pi = 4 * np.sum(y) / N #面積を計算
    return pi

#出力
print("Rectangle(100)     = " + str(Rectangle(100)))
print("Rectangle(1000)    = " + str(Rectangle(1000)))
print("Rectangle(10000)   = " + str(Rectangle(10000)))
print("Rectangle(100000)  = " + str(Rectangle(100000)))
print("Rectangle(1000000) = " + str(Rectangle(1000000)))

結果は以下のようになります

Rectangle(100)     = 3.160417031779046
Rectangle(1000)    = 3.1435554669110277
Rectangle(10000)   = 3.141791477611322
Rectangle(100000)  = 3.1416126164019866
Rectangle(1000000) = 3.141594652413813

多角形近似に比べると収束は遅いですが、等分数 $n$ の増加にしたがって値が円周率に収束しています。

5. 数値積分(台形近似)

　先ほどの区分求積にひと工夫を加えるだけで、長方形近似より精度がかなり高くなります。先ほどは微小区間の面積を長方形で近似していましたが、それを上底 $f (x_{i})$ 、下底 $f (x_{i + 1})$ 、高さ $Δ x$ の台形として近似します。つまりは以下の公式を用います。

\begin{matrix} (7) & S_{n} = \sum_{i = 0}^{n - 1} \frac{f (x_{i}) + f (x_{i + 1})}{2} \cdot \frac{b - a}{n} \to S (n \to \infty) \end{matrix}

以下の図を見れば分かるように、長方形近似では上部に余計な、あるいは足りない部分がかなりありますが、台形にすることでかなり減少します。

では、実装してみましょう！

import numpy as np

#台形近似
def Trapezoid(N):
    x = np.arange(N + 1) / N #0~1までのN＋1要素配列
    y = np.sqrt(1 - x**2) #y=root(1-x^2)を計算
    z = (y[0:N] + y[1:N+1]) / 2 #台形の(上底＋下底)/2を計算
    pi = 4 * np.sum(z) / N #面積を計算
    return pi

#出力
print("Trapezoid(100)     = " + str(Trapezoid(100)))
print("Trapezoid(1000)    = " + str(Trapezoid(1000)))
print("Trapezoid(10000)   = " + str(Trapezoid(10000)))
print("Trapezoid(100000)  = " + str(Trapezoid(100000)))
print("Trapezoid(1000000) = " + str(Trapezoid(1000000)))

以下が実行結果です。多角形近似に比べると僅かに劣りますが、長方形近似に比べると収束が速くなっていることが分かります。

Trapezoid(100)     = 3.1404170317790454
Trapezoid(1000)    = 3.141555466911028
Trapezoid(10000)   = 3.141591477611322
Trapezoid(100000)  = 3.1415926164019865
Trapezoid(1000000) = 3.1415926524138094

級数から求める

　無限に続く数列の総和である無限級数を用いて円周率を表す方法が存在します。それを利用して円周率を求めてみましょう！（尚、証明はかなり面倒或いは筆者には無理なので割愛します。。。）

6. バーゼル級数

　以下の級数で円周率が出現することがわかっています。

\begin{matrix} (8) & 1 + \frac{1}{2^{2}} + \frac{1}{3^{2}} + \dots + \frac{1}{n^{2}} + \dots = \sum_{n = 1}^{\infty} \frac{1}{n^{2}} = \frac{π^{2}}{6} \end{matrix}

自然数の２乗の逆数を全て足し合わせた式ですが、意外や意外、なんと結果に円周率が含まれています。数学をやっていると、このような「そこで $π$ 出てくるんかい！」といった楽しみが味わえて面白いですね！
バーゼル級数の詳細はこちら

　数値積分の時と同様に、コンピューターでは無限に続く演算を行うことは不可能なので、ここでは途中で和を打ち切って近似値を求めます。では、実装してみましょう！

import numpy as np

#バーゼル級数。第N項まで和を計算する。
def Basel(N):
    x = np.arange(1, N + 1) #1~Nまでの自然数の配列
    pi = np.sqrt(6 * np.sum(1 / x**2)) #配列の各要素の逆２乗和を計算し、６をかけて平方する。
    return pi

#出力
print("Basel(100)     = " + str(Basel(100)))
print("Basel(1000)    = " + str(Basel(1000)))
print("Basel(10000)   = " + str(Basel(10000)))
print("Basel(100000)  = " + str(Basel(100000)))
print("Basel(1000000) = " + str(Basel(1000000)))

以下が結果です。和を取る項数を増やすに従って円周率に収束しているのがわかります。

Basel(100)     = 3.132076531809106
Basel(1000)    = 3.1406380562059932
Basel(10000)   = 3.14149716394721
Basel(100000)  = 3.1415831043264415
Basel(1000000) = 3.1415916986604673

7. ライプニッツ級数

　円周率が出てくる級数は他にもあります。以下のライプニッツ級数もその１つです。

\begin{matrix} (9) & 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \dots + \frac{(- 1)^{n - 1}}{2 n - 1} + \dots = \sum_{n = 1}^{\infty} \frac{(- 1)^{n - 1}}{2 n - 1} = \frac{π}{4} \end{matrix}

今度は奇数の逆数をプラスマイナス交互に足して行く計算ですが、なんとここにも円周率が出現します！
ライプニッツ級数の詳細はこちら

　バーゼル級数と同様に実装していきましょう。

import numpy as np

#ライプニッツ級数。第N項まで和を計算する。
def Leibniz(N):
    x = np.arange(1, N + 1) #1~Nまでの自然数の配列
    y = 1 / (2 * x - 1) #奇数の逆数の配列を計算
    pm = (-1) ** (x - 1) #プラマイ１が交互に並ぶ数列の計算
    pi = 4 * np.dot(y, pm) #yとpmの内積によりプラマイ交互の和を計算し、最後に４をかける
    return pi

#出力
print("Leibniz(100)     = " + str(Leibniz(100)))
print("Leibniz(1000)    = " + str(Leibniz(1000)))
print("Leibniz(10000)   = " + str(Leibniz(10000)))
print("Leibniz(100000)  = " + str(Leibniz(100000)))
print("Leibniz(1000000) = " + str(Leibniz(1000000)))

以下が実行結果です。順調に収束していますね！

Leibniz(100)     = 3.131592903558553
Leibniz(1000)    = 3.140592653839795
Leibniz(10000)   = 3.141492653590045
Leibniz(100000)  = 3.141582653589825
Leibniz(1000000) = 3.141591653589752

8. ラマヌジャン級数

　最後に、化け物じみた級数「ラマヌジャン級数」を紹介します。以下がその式です。

\begin{array}{r} (10) & \frac{1}{π} = \frac{2 \sqrt{2}}{99^{2}} \sum_{n = 0}^{\infty} \frac{(4 n)! (1103 + 26390 n)}{(396^{n} n!)^{4}} \\ (11) & \frac{4}{π} = \sum_{n = 0}^{\infty} \frac{(- 1)^{n} (4 n)! (1123 + 21460 n)}{882^{2 n + 1} (4^{n} n!)^{4}} \end{array}

もはや意味がわからないですね。。。ラマヌジャンはインドの数学者で、上の式のような数々の意味不明な式を残しています。また、「式をどのように思いついたか」と聞かれると「夢の中で女神が教えてくれる」と答えたそうです。変態です。
ラマヌジャン級数の詳細はこちら
 ラマヌジャンについてはこちら

　それでは、実装してみましょう！今回は気分で(11)式を実装します。

import numpy as np

#ラマヌジャン級数。第N項まで和を取る。
def Ramanujan(N):
    sum = 0.0 #和の値。初期値０
    for i in range(N): #iが０からN-1までの間以下を実行
        numerator = ((-1)**i) * np.math.factorial(4 * i) * (1123 + 21460 * i) #第i項の分子の計算
        denominator = (882**(2 * i + 1)) * ((4**i) * np.math.factorial(i))**4 #第i項の分母の計算
        sum = sum + numerator / denominator #sumに第i項を加算
    pi = 4 / sum #円周率の計算
    return pi

#出力
print("Ramanujan(1) = " + str(Ramanujan(1)))
print("Ramanujan(2) = " + str(Ramanujan(2)))
print("Ramanujan(3) = " + str(Ramanujan(3)))
print("Ramanujan(4) = " + str(Ramanujan(4)))
print("Ramanujan(5) = " + str(Ramanujan(5)))

　ラマヌジャン級数の収束が非常に速いこと、そして式に登場する階乗( $n!$ とか)は発散が急激で大きい値を入れるとオーバーフローする恐れがあることより、和を取る項はバーゼル級数やライプニッツ級数よりも遥かに少なくしています。

　以下が結果です。第１項だけでも少数第４桁まで一致し、第３項以降はパソコンに表示される桁数の範囲で同じ値になってしまいました。圧倒的な収束の速さですね！

Ramanujan(1) = 3.1415850400712375
Ramanujan(2) = 3.141592653597622
Ramanujan(3) = 3.1415926535897936
Ramanujan(4) = 3.1415926535897936
Ramanujan(5) = 3.1415926535897936

乱数を使う

　乱数を使うことで $π$ を求めることもできます。乱数とは、その名の通りランダムな数で、たとえばサイコロの出る目などが該当します。しかし、乱数を用いて $π$ を求めるには大量のサンプル( $\sim$ 数百万)が必要で、人間が１００万回もサイコロを降るのはあまりにも退屈です。めんどくさいルーティンワークはコンピュータを導入して効率化しましょうということで、以下では乱数を用いて $π$ を求める方法を紹介していきます！

9. モンテカルロ法

　 $1 \times 1$ の正方形と、その頂点の1つを中心とする四分円を考えます。正方形と四分円の面積比は $1 : π / 4$ です。この正方形内に、完全にランダムに点を打っていきます。点が多くなると、「点の総数：四分円の内部に含まれる点の数 $≃ 1 : π / 4$ 」となります。例えば、 $N$ 個ののうち $n$ 個の点が四分円に含まれた場合、

\begin{matrix} (12) & \frac{n}{N} ≃ \frac{π}{4} \end{matrix}

となります。これを利用し、円周率を求めます。

では、実装してみましょう！

import numpy as np

#モンテカルロ法。N個のサンプルを取る。
def MonteCarlo(N):
    xy = np.random.rand(2, N) #区間[0,1)の一様分布乱数の(2×N)行列。各列が平面上の点の座標に対応。
    r = np.sum(xy**2, axis = 0) #各サンプル点の原点からの距離を計算。各成分２乗し、同じ列の要素の和を取る。
    r = r[r<=1] #距離が1以下（=四分円内にある）サンプルのみを抽出

    n = r.size #円内にある点の数を計算 #円周率を計算

    pi = 4 * n / N #円周率を計算
    return pi

np.random.seed(seed=0) #乱数の種の設定（これをすると出力が毎回同じになる）

#出力
print("MonteCarlo(100)       = " + str(MonteCarlo(100)))
print("MonteCarlo(1000)      = " + str(MonteCarlo(1000)))
print("MonteCarlo(10000)     = " + str(MonteCarlo(10000)))
print("MonteCarlo(100000)    = " + str(MonteCarlo(100000)))
print("MonteCarlo(1000000)   = " + str(MonteCarlo(1000000)))
print("MonteCarlo(10000000)  = " + str(MonteCarlo(10000000)))

　以下が結果です。一般に乱数による値の誤差はサンプル数 $N$ に対して $1 / \sqrt{N}$ のオーダーで収束する(割と遅い)ので、なかなか精度は低いです。

MonteCarlo(100)       = 3.24
MonteCarlo(1000)      = 3.068
MonteCarlo(10000)     = 3.1508
MonteCarlo(100000)    = 3.13576
MonteCarlo(1000000)   = 3.140048
MonteCarlo(10000000)  = 3.1424264

10. ビュフォンの針

　長さ $l$ の針を間隔 $t$ でたくさん引かれた平行線の上にたくさん落とします。 $N$ 個の針のうち線と交わった針の本数を $n$ とします。針の数 $N$ が十分大きいとき、

\begin{matrix} (13) & \frac{2 l N}{t n} \sim π \end{matrix}

となることが知られています。また来ましたね、「ここで $π$ 出てくるんかい！」のやつですw ただ平行線の上に針を落とすだけなのに、その確率には円周率が出現するという一見不思議なことが起こります。
ビュフォンの針の詳細はこちら

　では、実装してみましょう！針の中心と線との距離は $0 \sim t / 2$ で一様、また針と平行線との角度は $0^{\circ} \sim 90^{\circ}$ で一様とします。また、今回は針の長さ $l = 2$ 、平行線の間隔 $t = 4$ とします。

　針の角度を乱数からサンプルする上で気をつけなければならないことがあります。それは、直接角度をサンプルする場合、度数法を弧度法に変換する必要がありますが、この際に用いるであろうnp.radius()は円周率の値に依存した関数となっています。それを使ってもシミュレーションはできるのですが、個人的に気持ち悪いので、今回は少し工夫をします。単位四分円中の点の座標をサンプルし、そこからsinの値を求めるという方法をとります。サンプル点の範囲を単位四分円の中に限定することで、点と $x$ 軸のなす角度が
$0^{\circ} \sim 90^{\circ}$ の一様な分布となります。

import numpy as np

#ビュフォンの針。N本の針を落とす。
def Buffon(N):
    n = 0; #線に重なる針の本数。初期値０。
    for i in range(N): #iが0からN-1までの間、以下を繰り返す。
        #針の角度のサンプリング。角度を直接サンプルする代わりに単位円内の点をサンプルすることでpiの値依存性を排除。
        r = 2 #サンプル点の原点からの距離。以下のwhile文を行うため初期値を2に。
        while r > 1: #r<=1となるまでサンプルを繰り返す
            dx = np.random.rand() #x座標
            dy = np.random.rand() #y座標
            r = np.sqrt(dx ** 2 + dy ** 2) #原点からの距離の計算

        h = 2 * np.random.rand() + dy / r #房の先端の高さ（高い方）の計算
        if h > 2: #針の先端の高さが平行線の高さを終えた場合
            n += 1 #線に重なる針の本数を加算

    pi = N / n #円周率の計算
    return pi

np.random.seed(seed=0) #乱数の種の設定（これをすると出力が毎回同じになる）

#出力
print("Buffon(100)        = " + str(Buffon(100)))
print("Buffon(1000)       = " + str(Buffon(1000)))
print("Buffon(10000)      = " + str(Buffon(10000)))
print("Buffon(100000)     = " + str(Buffon(100000)))
print("Buffon(1000000)    = " + str(Buffon(1000000)))
print("Buffon(10000000)   = " + str(Buffon(10000000)))

　theta = np.radians(rand.rand(N) * 90) で $N$ 本の針の角度を生成します。rand.rand() は区間 $[0, 1)$ の一様分布なので、90をかけることで $[0, 90)$ の一様分布になります。

　次に、y = 2 * rand.rand(N) + np.sin(theta) で針の先端の $y$ 座標を計算します。2 * rand.rand(N) で針の中心の座標を、np.sin(theta) で針の中心と先端の高さの差を出します。針の先端が平行線と交わっているかどうかは [y>2] で判断します( $y > 2$ なら交わっている、そうでなければ交わっていない)。あとは y.size で交わっている針の本数を計算し、ごちゃごちゃいじれば終了です。

　以下が出力結果です。なかなか結果が安定しませんね。。。(汗)

Buffon(100)        = 2.7777777777777777
Buffon(1000)       = 2.881844380403458
Buffon(10000)      = 3.1259768677711786
Buffon(100000)     = 3.11284046692607
Buffon(1000000)    = 3.1381705093564554
Buffon(10000000)   = 3.1422055392056127

アルゴリズムを使う

　アルゴリズムとは、ざっくり言えばある特定の手順をまとめたものです。以下では、特定の手順を繰り返すことで円周率を求める方法を紹介します！

11.ガウス=ルジャンドルのアルゴリズム

　ガウス=ルジャンドルのアルゴリズムとは以下のように「連立漸化式に従って数列の項を繰り返し計算していく」という手順に従って円周率を計算する方法です。

\begin{matrix} (14) & \begin{aligned} a_{0} = 1, b_{0} = \frac{1}{\sqrt{2}}, t_{0} = \frac{1}{4}, p_{0} = 1 \\ a_{n + 1} = \frac{a_{n} + b_{n}}{2} \\ b_{n + 1} = \sqrt{a_{n} b_{n}} \\ t_{n + 1} = t_{n} - p_{n} (a_{n} - a_{n + 1})^{2} \\ p_{n + 1} = 2 p_{n} \\ \frac{a_{n} + b_{n}}{4 t_{n}} \to π (n \to \infty) \end{aligned} \end{matrix}

　このアルゴリズムは、計算を進めていくたびに厳密値と一致する桁数がおよそ２倍に増えることが知られており、簡単かつ高精度で円周率を求めることができます。
　このアルゴリズムをそれぞれ独立に考え出したガウスとルジャンドルは、二人とも数学・物理界のスーパースターです。
ガウス=ルジャンドルのアルゴリズムの詳細はこちら
 ガウスの詳細はこちら
 ルジャンドルの詳細はこちら

　では、実装していきましょう！

import numpy as np

#ガウス=ルジャンドルのアルゴリズム。N回繰り返す
def GaussRegendre(N):
    #初項
    a = 1.0
    b = np.sqrt(2) / 2
    t = 1 / 4
    p = 1

    for i in range(N): #iが0からN-1までの間以下を繰り返す
        #漸化式に従い、次の項を計算
        a_new = (a + b) / 2
        b_new = np.sqrt(a * b)
        t_new = t - p * (a - a_new)**2
        p_new = 2 * p

        #古い値を新しい値に置き換え
        a = a_new
        b = b_new
        t = t_new
        p = p_new

    pi = (a + b)**2 / (4 * t) #円周率の計算
    return pi

#出力
print("GaussRegendre(1) = " + str(GaussRegendre(1)))
print("GaussRegendre(2) = " + str(GaussRegendre(2)))
print("GaussRegendre(3) = " + str(GaussRegendre(3)))

以下が結果です。１回目は小数第２桁まで、２回目は小数第７桁まで、３回目は表示されている最後の桁まで円周率に一致しています。収束がかなり速いですね！

GaussRegendre(1) = 3.1405792505221686
GaussRegendre(2) = 3.141592646213543
GaussRegendre(3) = 3.141592653589794

物理を使う

　数学の世界だけでなく、自然界にも円周率が隠れています。以下では、物理を利用して円周率を計算する方法を紹介します！

12. 単振動の周期

　物理学において、振動をはじめとする盛んに研究されてきた運動ですが、全ての周期運動の基本となるのが単振動です。単振動は、位置 $x$ が以下の方程式に従う場合の運動です。

\begin{matrix} (15) & \frac{d^{2} x}{d t^{2}} = - ω^{2} x \end{matrix}

(15)式のように、関数の微分を含む方程式を微分方程式と呼びます。世にある多くの微分方程式は厳密に解くことができませんが、(15)式は厳密に解くことができ、解は以下のようになります。

\begin{matrix} (16) & x = A \sin (ω t + ϕ) \end{matrix}

位置 $x$ が時刻 $t$ の関数として $\sin$ で表されるので、時間経過に従って $x$ が振動することがわかります。 $A$ は振動の振幅、 $ϕ$ は初期位相と呼ばれる定数で、物体の $t = 0$ での位置に関係する量です。そして、この振動の周期は

\begin{matrix} (17) & T = \frac{2 π}{ω} \end{matrix}

となります。 $ω = 1$ とすれば、振動の半周期が円周率となります。

　今回は、微分方程式を数値的にシミュレートすることで振動の周期を計算していきます。微分方程式のシミュレートには、ルンゲ≡クッタ法を用います。ルンゲ≡クッタ法の詳しい説明は長いのでここでは割愛しますが、簡単に説明すると、「ある一定の手順に従って、ある時刻での位置や速度の値から少し未来での位置や速度の値を求めるということを繰り返す方法」になります。

\begin{matrix} (18) & \begin{aligned} 時 間 間 隔 : h, 微 分 方 程 式 : \frac{d x}{d t} = f (x, t) \\ k_{1} = h \times f (x_{i}, t) \\ k_{2} = h \times f (x_{i} + \frac{k_{1}}{2}, t + \frac{h}{2}) \\ k_{3} = h \times f (x_{i} + \frac{k_{2}}{2}, t + \frac{h}{2}) \\ k_{4} = h \times f (x_{i} + k_{3}, t + h) \\ k = \frac{k_{1} + 2 k_{2} + 2 k_{3} + k_{4}}{6} \\ x_{i + 1} = x_{i} + k \end{aligned} \end{matrix}

ルンゲ=クッタ法の詳しい説明はこちら

ルンゲ=クッタ法を用いた別の記事も書いているのでよかったらお読みください！
ルンゲクッタ法を用いたUnityでのカオスのシミュレーション

　では、実装していきましょう！(18)式のルンゲ=クッタ法はデフォルトでは１階微分方程式に用いられますが、(15)式の微分方程式は２階の微分方程式なので、以下のように２つの１階微分方程式に分けてしまいます。

\begin{matrix} (19) & \begin{aligned} \frac{d x}{d t} = v \\ \frac{d v}{d t} = - ω^{2} x \end{aligned} \end{matrix}

初期位置及び初速度は $(x_{0}, v_{0}) = (1, 0)$ とします。こうすると、物体の運動は振動の最高点から始まります。その後しばらく位置が減少し、丁度半周期経った時に振動の最下点に達し、その後は位置が上昇します。そうすると、物体の位置が減少している間だけ時刻を進め、位置が増加したところで打ち切ることで、半周期の時刻(≒円周率)が得られます。

import numpy as np

#微分方程式。r:位置と速度をまとめたベクトル
def f(r):
    x = r[0] #位置の値を取り出す
    v = r[1] #速度の値を取り出す

    #微分方程式の右辺を計算
    xdot = v
    vdot = -x

    return np.array([xdot, vdot])

#ルンゲ=クッタ法
def RungeKutta(r, h):
    k1 = h * f(r)
    k2 = h * f(r + k1 / 2)
    k3 = h * f(r + k2 / 2)
    k4 = h * f(r + k3)

    k = (k1 + 2 * k2 + 2 * k3 + k4) / 6

    return r + k

#単振動。時間間隔h=10^(-N)でルンゲ・クッタ法を行う。
def Vibration(N):
    h = 10**(-N) #hの計算
    t = 0 #時刻に対応する整数の初期値
    isDecrease = True #位置が減少しているかどうかのフラグ
    r = np.array([1.0, 0.0]) #位置・速度の初期値
    x_old = r[0] #古い位置の初期値代入。位置の増減判定に使う。

    while isDecrease == True and t * h < 4.0: #位置が減少かつ時刻が4以下の間、以下を繰り返す
    #※t*h<4.0を課すことで万が一無限ループに入ってしまうのを防止

        r = RungeKutta(r, h) #ルンゲ=クッタ法で位置・速度を更新
        x_new = r[0] #新しい位置の値

        if x_old > x_new: #古い位置が新しい位置より大きい（=位置が減少している）場合
            x_old = x_new #古い位置の値を置き換える
            t += 1 #時刻を進める
        else : #それ以外の場合
            isDecrease = False #減少しているというフラグをFalseに→ループが終了
    pi = t / 10 ** N
    return pi #半周期(≒π)を返す

#出力
print("Vibration(1) = " + str(Vibration(1)))
print("Vibration(2) = " + str(Vibration(2)))
print("Vibration(3) = " + str(Vibration(3)))
print("Vibration(4) = " + str(Vibration(4)))
print("Vibration(5) = " + str(Vibration(5)))

以下が結果です。 $h$ の幅を小さくするほど精度が上がっていることがわかります。しかし、それに伴い必要なループ数も増加するので、時間がかかってきます。

Vibration(1) = 3.1
Vibration(2) = 3.14
Vibration(3) = 3.142
Vibration(4) = 3.1416
Vibration(5) = 3.14159

13. 2つの球の衝突回数

　１３番目、すなわち最後の方法として、円周率を求める衝撃の方法を紹介します。なんと、２つの球と壁との衝突回数が円周率になるとのことです。何を言っているのかさっぱりわからないと思うので、動画と原論文を貼っておきます。

動画：

原論文：PLAYING POOL WITH π (THE NUMBER π FROM A BILLIARD POINT OF VIEW)

質量比が $1 : 100^{N}$ の２球を用意し、軽い球を静止した状態で重い球を速さ1で軽い球にぶつけます。衝突後の軽い球が進む先に壁を用意し、反発係数1で跳ね返るようにします。衝突を繰り返すと、重い球が壁と反対方向に進むようになり、軽い球が重い球に追いつけなくなったところで衝突しなくなります。ここまでの、軽い球が重い球及び壁と衝突した回数を数えると、円周率の $N$ 桁分の数字になるということです。

　では、実装してみましょう！まず、衝突による速度の変化を事前に計算しておきます。計算方法は、高校物理でもお馴染みの運動量保存則と反発係数の定義式の連立です。衝突前の軽い球、重い球の速度をそれぞれ $v, V$ 、衝突後の速度をそれぞれ $v^{'}, V^{'}$ とします。軽い球と重い球の質量比を $1 : r$ とします。この時、衝突前後で成り立つ方程式は以下の２つです。

\begin{matrix} (20) & \begin{aligned} v + r V = v^{'} + r V^{'} \\ e = - \frac{v^{'} - V^{'}}{v - V} \end{aligned} \end{matrix}

今回、衝突は全て弾性衝突なので、 $e = 1$ となります。連立方程式を $v^{'}, V^{'}$ について解くと、

\begin{matrix} (21) & v^{'} = \frac{(1 - r) v + 2 r V}{1 + r}, V^{'} = \frac{2 v + (r - 1) V}{1 + r} \end{matrix}

となります。これを踏まえて、実装していきましょう！

#２球の衝突。第N桁まで求める。
def BallCollision(N):
    v = 0.0  #軽球の初速
    V = -1.0 #重球の初速
    r = 100 ** N #質量比
    count = 0 #衝突回数の初期値

    while v > V: #軽球の速度が重球の速度より大きい間、以下を繰り返す
        v_new = ((1 - r) * v + 2 * r * V) / (r + 1) #衝突後の軽球の速度
        V_new = (2 * v + (r - 1) * V) / (r + 1) #衝突後の重球の速度
        count += 1 #衝突数追加

        if(v_new < 0.0): #軽球の速度が負の場合
            v_new = -v_new #壁との衝突により速度が反転
            count += 1 #衝突数追加

        v = v_new #軽球の速度の置き換え
        V = V_new #重球の速度の置き換え

    pi = count / (10 ** N) #円周率の計算
    return pi

#出力
for i in range(8):
    print("BallCollision(" + str(i) + ") = " + str(BallCollision(i)))

以下が実行結果です。本当に円周率がN桁目まで求まっていますね！しかし、精度を１桁増やすのに計算量が１００倍になるので、使い勝手は悪そうですねw

BallCollision(0) = 3.0
BallCollision(1) = 3.1
BallCollision(2) = 3.14
BallCollision(3) = 3.141
BallCollision(4) = 3.1415
BallCollision(5) = 3.14159
BallCollision(6) = 3.141592
BallCollision(7) = 3.1415926

クラスに全部まとめてみた

　最後に、今まで実装してきた手法を全てクラスにまとめてみました！こうして全部集めてみると、かなりのボリュームですね。。。！

#クラスPi_Haremの説明
#    円周率を求める１３通りの方法が織りなすハーレムclass。
#    好みの方法で円周率を求めよう！
#    ※使用には import numpy as np が必要。
#
#関数の説明：
#    関数は全て返り値を持つ。
#    1.  pi()　: 
#        np.piの値を返す。引数なし。
#
#     2.  Arctan()　: 
#        4*np.arctan(1)の値を返す。引数なし。
#
#    3.  RegularPokygon(N) : 
#       正多角形の面積から円周率を求め、値を返す。引数N＝正多角形の頂点の数。
#
#    4.  Rectangle(N) : 
#        数値積分（長方形近似）から円周率を求め、値を返す。引数N＝積分区間の等分数。
#
#    5.  Trapezoid(N) : 
#        数値積分（台形近似）から円周率を求め、値を返す。引数N＝積分区間の等分数。
#
#    6.  Basel(N) : 
#        バーゼル級数から円周率を求め、値を返す。引数N＝和を取る項数。
#
#    7.  Leibniz(N) : 
#        ライプニッツ級数から円周率を求め、値を返す。引数N＝和を取る項数。
#
#    8.  Ramanujan(N) : 
#        ラマヌジャン級数から円周率を求め、値を返す。引数N＝和を取る項数。
#
#    9.  MonteCarlo(N) : 
#        四分円のモンテカルロ求積から円周率を求め、値を返す。引数N＝サンプル数。
#
#    10. Buffon(N) : 
#        ビュフォンの針のシミュレーションから円周率を求め、値を返す。引数N＝サンプル数。
#
#    11. GaussRegendre(N) : 
#        ガウス=ルジャンドルのアルゴリズムから円周率を求め、値を返す。引数N＝アルゴリズムの反復回数。
#
#    12. Vibration(N) : 
#        単振動のシミュレーションから円周率を求め、値を返す。引数N＝時間の刻み幅hを決める変数。h=10^(-N)
#
#    13. BallCollision(N) : 
#        ２球の衝突のシミュレーションから円周率を求め、値を返す。引数N＝求めたい桁数。

import numpy as np

#クラスの定義
class Pi_Harem():
    #Numpyの基本機能
    # 1.πの厳密値(numpy.piを使用)
    def pi(self):
        return np.pi

    # 2.tanの逆関数
    def Arctan(self):
        return 4 * np.arctan(1)

    #面積
    # 3.多角形近似
    def RegularPolygon(self, N):
        theta = np.radians(360 / N)
        pi = N * np.sin(theta) / 2

        return pi

    # 4.長方形近似
    def Rectangle(self, N):
        x = np.arange(N) / N
        y = np.sqrt(1 - x**2)
        pi = 4 * np.sum(y) / N

        return pi

    # 5.台形近似
    def Trapezoid(self, N):
        x = np.arange(N + 1) / N
        y = np.sqrt(1 - x**2)
        z = (y[0:N] + y[1:N+1]) / 2    
        pi = 4 * np.sum(z) / N

        return pi

    #級数
    # 6.バーゼル級数
    def Basel(self, N):
        x = np.arange(1, N + 1)
        pi = np.sqrt(6 * np.sum(1 / x**2))

        return pi

    # 7.ライプニッツ級数
    def Leibniz(self, N):
        x = np.arange(1, N + 1)
        y = 2 * x - 1
        pi = 4 * np.dot(1 / y, (-1)**(x - 1))

        return pi

    # 8.ラマヌジャン級数    
    def Ramanujan(self, N):
        sum = 0.0

        for i in range(N):
            numerator = ((-1)**i) * np.math.factorial(4 * i) * (1123 + 21460 * i)
            denominator = (882**(2 * i + 1)) * ((4**i) * np.math.factorial(i))**4
            sum = sum + np.sum(numerator / denominator)

        pi = 4 / sum

        return pi

    #乱数
    # 9.モンテカルロ法
    def MonteCarlo(self, N):
        xy = np.random.rand(2, N)
        r = np.sum(xy**2, axis = 0)
        r = r[r<=1]
        n = r.size

        pi = 4 * n / N

        return pi

    # 10.ビュフォンの針
    def Buffon(self, N):
        n = 0;
        for i in range(N):
            r = 2
            while r > 1:
                dx = np.random.rand()
                dy = np.random.rand()
                r = np.sqrt(dx ** 2 + dy ** 2)

            h = 2 * np.random.rand() + dy / r
            if h > 2:
                n += 1

        pi = N / n
        return pi

    #アルゴリズム
    # 11.ガウス=ルジャンドルのアルゴリズム
    def GaussRegendre(self, N):
        a = 1.0
        b = np.sqrt(2) / 2
        t = 1 / 4
        p = 1

        for i in range(N):
            a_new = (a + b) / 2
            b_new = np.sqrt(a * b)
            t_new = t - p * (a - a_new)**2
            p_new = 2 * p

            a = a_new
            b = b_new
            t = t_new
            p = p_new

        pi = (a + b)**2 / (4 * t)
        return(pi)

    #物理
    # 12.単振動
    #　　微分方程式(Private関数)
    def __f(self, r):
        x = r[0]
        v = r[1]
        xdot = v
        vdot = -x
        return np.array([xdot, vdot])

    #　　ルンゲクッタ法(Private関数)
    def __RungeKutta(self, r, h):
        k1 = h * self.__f(r)
        k2 = h * self.__f(r + k1 / 2)
        k3 = h * self.__f(r + k2 / 2)
        k4 = h * self.__f(r + k3)

        k = (k1 + 2 * k2 + 2 * k3 + k4) / 6

        return r + k

    #　　　周期の計算
    def Vibration(self, N):
        h = 10 ** (-N)
        t = 0
        isDecrease = True
        r = np.array([1.0, 0.0])
        x_old = r[0]

        while isDecrease == True and t * h < 4.0:
            r = RungeKutta(r, h)
            x_new = r[0]

            if x_old > x_new:
                x_old = x_new
                t += 1
            else :
                isDecrease = False

        pi = t / 10 ** (N)
        return pi

    # 13.2球の衝突
    def BallCollision(self, N):
        v = 0.0
        V = -1.0
        r = 100 ** N
        count = 0
        while v > V:
            v_new = ((1 - r) * v + 2 * r * V) / (r + 1)
            V_new = (2 * v + (r - 1) * V) / (r + 1)
            count += 1

            if(v_new < 0.0):
                v_new = -v_new
                count += 1

            v = v_new
            V = V_new

        pi = count / (10 ** N)
        return pi

    #以上

jupyter notebook の場合、上のコードをあるセルで実行してから別のセルで以下のようなコードを実行して使用できます(クラス中のBasel()関数を使用する場合)。

pi_harem = Pi_Harem() #Pi_Harem()を変数pi_haremに代入（一度だけ行えば良い）
pi = pi_harem.Basel(1000) #変数piにクラスの中にある関数Barsel(1000)の値を代入
print(pi) #出力

まとめ

　いかがでしたか？今回は Python で円周率を計算する方法を13通り扱いました。かなりのボリュームで、執筆にも大変な時間がかかりました(汗)

　王道の方法から非自明な方法、そしてクセのある方法がたくさんありましたね！どの方法たちも見ていて飽きることがなく、個人的には甲乙付け難かったです。。。

　今後も物理やプログラミングを中心に記事を投稿していく予定です。もし何か扱って欲しいテーマがあればコメント欄にお願いします！

Qiitaにログインしてダークテーマを使ってみませんか？🌙

【ハーレム】多すぎて選べない！Pythonで円周率πを計算する13の方法