円周率を求めてみる【Python】

Python のプログラミングの一例として、自力で円周率を求めるスクリプトについてまとめておきます。
以下の２つの環境で動作確認をしています。

環境：
・　Windows パソコン、Python 3.x
・　Raspberry Pi (Raspbian Linux)

背景　～　目指せ、π計算 100兆桁！

円周率計算の世界記録は100兆桁を超えるそうです（Google Cloud）。
計算は、Debian Linux で実行されているとのことで、Raspberry Pi と同等ということになります。

ということで、円周率を計算してみます。
円周率を、どれくらいの時間で何桁まで計算できるか実測してみると、お使いの Linux / 計算機環境の素の実力を把握することが可能です。

なお、円周率を求める（π計算）のではなく、Python に保存されている円周率の値をたんに呼び出して使うだけの場合は、以下のスクリプトで、表示・利用が可能です。

Python で円周率 π を使う方法　～　mathパッケージを使う場合

import math
pi1 = math.pi 
print( pi1 )                                          # >> 3.141592653589793 

radius1 = 10.0 
S1 = pi1 * radius1 ** 2 
print( S1 )                                          # >> 314.1592653589793

円周率（π）を求める計算式

さて、円周率を自力で求める場合、いくつかの方法が考えられます。
どのやり方でもよいのですが、ここでは以下の公式（グレゴリー級数）を使うことにします。

グレゴリー級数の式

（式１）
$$ \frac{\pi}{4} = 1 – \frac{1}{3} + \frac{1}{5} – \frac{1}{7} + \frac{1}{9} – …$$
この式の証明は、末尾に概要を記載しておくことにします。
プログラミングのため、両辺に４を掛けて以下のように変形します。

（式２）
$$ \pi = \frac{4}{1} – \frac{4}{3} + \frac{4}{5} – \frac{4}{7} + \frac{4}{9} – …$$

サンプルコードと実行結果

上記の（式２）をそのまま、プログラムにします。

円周率 π を求める Python プログラム

pi1 = 0 
for i1 in range( 100000 ): 
    odd1 = 2.0*i1 + 1.0 
    v1 = ((-1.0)**i1) * 4.0/odd1
    pi1 = pi1 + v1 
    print( pi1 )

上記のプログラムを実行すると、3.14… の数字がずらずらと出てきて、だんだん１つの値（真値）に漸近する動きとなります。
プログラムを６行ほど書くだけで、円周率を計算できることがわかります。

サンプルコードの説明

・　pi1 が円周率です。
・　（式２）より、各項を無限に足していけば、円周率が求められます。
・　range() のカッコ内の数値が級数の項の数です。（式２）からは、無限個の数値を入れることになりますが、所定の有効桁で求められればよいため、まずは 10000 としています。
・　各項の分母は奇数となっていますので、最初に奇数 odd1 を作ります。
・　つぎに、（式２）の符号と分数の記載に基づき、各項の値を作ります。
この各項を pi1 に足しこんでいけば、円周率を求めることができます。
・　range() 内の数値は、最初は少ない数値を入れてスクリプトを実行してみてください。パソコンの計算能力に応じて大きな値に修正し、計算精度を上げてください。
・　（式２）の右辺は、絶対値が次第に小さくなる符号の異なる数を足し続ける形となっています。
したがって、円周率の真の値は、プログラムの実行中に最後に示される２つの数の間に存在することになります。加えて、表示される数値の変動幅はかならず小さくなっていくことになります。
この結果、計算中に表示される数値は、かならず、真の値に近づいていき、確定した桁数が次第に増えていく、といった動きになります。
・　グレゴリー級数は収束が遅いといわれます。その反面、グレゴリー級数を使うメリットとしては、計算式がシンプルであることに加え、加えていく項の絶対値が次第に小さくなっていくため、計算結果が収束していくという点にあるといえます。
例えば、円周率を計算するのに、乱数を使う方法も考えられます。しかし乱数を使うと、乱数の出方によって都度、計算結果が変わります。また、得られた値の何桁目までが信頼性のある値なのか、わからなくなると思います。グレゴリー級数を使うと、浮動小数点演算の精度の範囲内で、最後に加える値の大きさ程度のばらつきを考慮しておけばよいことがわかるという点で、メリットがあるといえます。

式の証明（概要）

まず、$tan y$ を $ y $ で微分します。
$$ \frac{d}{dy}tan y = \frac{1}{cos^2 y }
= \frac{cos^2 y + sin^2 y}{cos^2 y } = 1+ tan^2 y $$
つぎに、上式で、$tan y = x $ とおき、逆数を取ります。
$$ \frac{dy}{dx} = \frac{1}{1+ x^2} $$
ここで、$ y = tan^{-1} x $ であり、$ | x | < 1$ と仮定すると、右辺は等比数列の和の形となっているので、以下となります。
$$ \frac{d}{dx} tan^{-1} x= 1-x^2+x^4-x^6+ … $$
両辺が$x $ で積分できるものと仮定（※）し、0 から $x$ まで積分します。
$$ tan^{-1} x= x-\frac{1}{3} x^3+\frac{1}{5} x^5-\frac{1}{7} x^7+ … $$
$x → 1$ とすると、 $tan \frac{\pi}{4} = 1$ ですので、（式１）が得られます。
※　途中で仮定を入れていますので、厳密な証明とはいえません。しかしながら、とくに問題なく円周率の計算はできます。