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

Python

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

環境:
・ Windows 10、Python 3.x
・ Raspberry Pi (Raspbian Linux)

背景 ~ 目指せ、100兆桁!

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

さて以前に、このサイト上で LaTex を使った数式が自由に扱えるようになりましたので、何か計算をしたくなってきます。

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

計算式

円周率の計算は、いくつか方法が考えられます。
どのやり方でもよいのですが、ここでは以下の公式(グレゴリー級数)を使うことにします。

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

(式2)
$$ \pi = \frac{4}{1} – \frac{4}{3} + \frac{4}{5} – \frac{4}{7} + \frac{4}{9} – …$$

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

上記の(式2)をそのまま、プログラムにします。

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… の数字がずらずらと出てきて、だんだん1つの値(真値)に漸近する動きとなります(詳細は省略します)。
円周率を計算するプログラムが6行ほどで書けることがわかります。

サンプルコードの説明

・ pi1 が円周率です。
・ (式2)より、各項を無限に足していけば、円周率が求められます。
・ range() のカッコ内の数値が級数の項の数です。(式2)からは、無限個の数値を入れることになりますが、所定の有効桁で求められればよいため、まずは 10000 としています。
・ 各項の分母は奇数となっていますので、最初に奇数 odd1 を作ります。
・ つぎに、(式2)の符号と分数の記載に基づき、各項の値を作ります。
この各項を pi1 に足しこんでいけば、円周率を求めることができます。
・ range() 内の数値は、最初は少ない数値を入れてスクリプトを実行してみてください。パソコンの計算能力に応じて大きな値に修正し、計算精度を上げてください。
・ (式2)の右辺は、絶対値が次第に小さくなる符号の異なる数を足し続ける形となっています。
したがって、円周率の真の値は、プログラムの実行中に最後に示される2つの数の間に存在することになります。加えて、表示される数値の変動幅はかならず小さくなっていくことになります。
この結果、計算中に表示される数値は、かならず、真の値に近づいていき、確定した桁数が次第に増えていく、といった動きになります。

式の証明(概要)

まず、\(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\) ですので、(式1)が得られます。
※ 途中で仮定を入れていますので、厳密な証明とはいえません。しかしながら、とくに問題なく円周率の計算はできます。

まとめ

Python での数値計算の事例として、円周率を計算してみました。
円周率を求める程度であれば、6行ほどのプログラムで書けることがわかります。
パソコンの性能評価などに使えそうです。

また、ざっと見たところ、LaTeX 形式で記載した数式もとくに違和感はなく、サイト上で数式が自由に扱えるようになりました。

とはいえ、グレゴリー級数を使った場合/いまのハードウェア環境の場合、100兆桁の計算ができるころには人生が終わっているような感じです。Debian ベースの Raspberry Pi は持っているのですけれども。。

関連リンク
・ ネイピア数を求めてみる 【Python】
・ Python で自由曲線を描く 【HTML & SVG】
・ 起動用SDカードを設定する 【Raspberry Pi】
・ WordPress で数式を扱う 【LaTeX】

タイトルとURLをコピーしました