ホーム / 技術開発(r&dディスクロージャ) / 基本物理定数と単位換算 / 円周率(pi) 100万(1,000,000)桁 円周率(PI) 100万(1,000,000)桁 10億桁はこちら(ダウンロード[ZIP,447MB]) GMP で浮動小数点数を扱う場合は MPFR(The GNU Multiple Precision Floating-Point Reliably) ライブラリが必要であるので、インストール済みであること。（参照：「MPFR - ソースビルドでインストール (on Linux Mint)！ …
ガウス＝ルジャンドルのアルゴリズムでπ計算 ガウスルジャンドルのアルゴリズムで、円周率を1000桁まで求める 前回、円周率を求めるプログラムをPythonで書いたが、途中までしか求めることが出来なかった。 idoushiki.hatenablog.com 64ビット浮動小数点で計算しているので、小数点以下14桁分まで求… 演算には GMP(The GNU Multi Precision Arithmetic Library) 任意精度算術ライブラリを Python 用にラップした gmp という PyPI ライブラリを使用するので、インストール済みであること。 4.
1. とりあえず円周率1000桁！ 3.1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164 0628620899 8628034825 3421170679 8214808651 3282306647 0938446095 5058223172 5359408128 4811174502 8410270193 8521105559 6446229489 5493038196 4428810975 6659334461 2847564823 3786783165 2712019091 いかがでしたか？今回は Python で円周率を計算する方法を13通り扱いました。かなりのボリュームで、執筆にも大変な時間がかかりました(汗)以下が実行結果です。多角形近似に比べると僅かに劣りますが、長方形近似に比べると収束が速くなっていることが分かります。円周率が出てくる級数は他にもあります。以下のライプニッツ級数もその１つです。以下が結果です。１回目は小数第２桁まで、２回目は小数第７桁まで、３回目は表示されている最後の桁まで円周率に一致しています。収束がかなり速いですね！長さ $l$ の針を間隔 $t$ でたくさん引かれた平行線の上にたくさん落とします。$N$ 個の針のうち線と交わった針の本数を $n$ とします。針の数 $N$ が十分大きいとき、位置 $x$ が時刻 $t$ の関数として $\sin$ で表されるので、時間経過に従って $x$ が振動することがわかります。$A$ は振動の振幅、$\phi$ は初期位相と呼ばれる定数で、物体の $t=0$ での位置に関係する量です。そして、この振動の周期は$1\times 1$ の正方形と、その頂点の1つを中心とする四分円を考えます。正方形と四分円の面積比は $1 : \pi/4$ です。この正方形内に、完全にランダムに点を打っていきます。点が多くなると、「点の総数：四分円の内部に含まれる点の数 $\simeq1 : \pi/4$」となります。例えば、$N$個ののうち $n$ 個の点が四分円に含まれた場合、では、実装していきましょう！(18)式のルンゲ=クッタ法はデフォルトでは１階微分方程式に用いられますが、(15)式の微分方程式は２階の微分方程式なので、以下のように２つの１階微分方程式に分けてしまいます。以下の図を見れば分かるように、長方形近似では上部に余計な、あるいは足りない部分がかなりありますが、台形にすることでかなり減少します。以下が実行結果です。本当に円周率がN桁目まで求まっていますね！しかし、精度を１桁増やすのに計算量が１００倍になるので、使い勝手は悪そうですねwガウス=ルジャンドルのアルゴリズムとは以下のように「連立漸化式に従って数列の項を繰り返し計算していく」という手順に従って円周率を計算する方法です。となりますが、プラスの方が上半分、マイナスの方が下半分を表します。今回は上半分のうちさらに $x>=0$ の部分(四分円)の面積を求めます。求める面積の区間は $[a,\ b]=[0,\ 1]$ となるので、これを区分求積の公式(4)式に当てはめると最後に、今まで実装してきた手法を全てクラスにまとめてみました！こうして全部集めてみると、かなりのボリュームですね。。。！数学の世界だけでなく、自然界にも円周率が隠れています。以下では、物理を利用して円周率を計算する方法を紹介します！以下が結果です。和を取る項数を増やすに従って円周率に収束しているのがわかります。ラマヌジャン級数の収束が非常に速いこと、そして式に登場する階乗($n!$とか)は発散が急激で大きい値を入れるとオーバーフローする恐れがあることより、和を取る項はバーゼル級数やライプニッツ級数よりも遥かに少なくしています。この方法を円の面積に応用しましょう。$x-y$ 平面における単位円(原点中心、半径1)の方程式はそれでは、実装してみましょう！今回は気分で(11)式を実装します。theta = np.radians(rand.rand(N) * 90) で $N$ 本の針の角度を生成します。rand.rand() は区間 $[0, 1)$ の一様分布なので、90をかけることで $[0, 90)$ の一様分布になります。数値積分の時と同様に、コンピューターでは無限に続く演算を行うことは不可能なので、ここでは途中で和を打ち切って近似値を求めます。では、実装してみましょう！以下が結果です。第１項だけでも少数第４桁まで一致し、第３項以降はパソコンに表示される桁数の範囲で同じ値になってしまいました。圧倒的な収束の速さですね！PythonにはNumpyと呼ばれるライブラリがあり(Anacondaの場合は初めから用意されています)、ベクトルや行列その他数学の計算において非常に便利な機能が備わっています。以下では、Numpyの基本的な機能で円周率を計算する方法を紹介します！では、実装してみましょう！針の中心と線との距離は $0\sim t/2$ で一様、また針と平行線との角度は $0^\circ \sim 90^\circ$ で一様とします。また、今回は針の長さ $l=2$、平行線の間隔 $t=4$ とします。となります。$\omega = 1$とすれば、振動の半周期が円周率となります。以下が結果です。$h$の幅を小さくするほど精度が上がっていることがわかります。しかし、それに伴い必要なループ数も増加するので、時間がかかってきます。無限に続く数列の総和である無限級数を用いて円周率を表す方法が存在します。それを利用して円周率を求めてみましょう！（尚、証明はかなり面倒或いは筆者には無理なので割愛します。。。）円に内接する正 $n$ 角形の面積を増やしていくことで、値がどんどん円の面積に近づきます。半径 $r$ の円に内接する正 $n$ 角形の面積は質量比が $1 : 100^N$ の２球を用意し、軽い球を静止した状態で重い球を速さ1で軽い球にぶつけます。衝突後の軽い球が進む先に壁を用意し、反発係数1で跳ね返るようにします。衝突を繰り返すと、重い球が壁と反対方向に進むようになり、軽い球が重い球に追いつけなくなったところで衝突しなくなります。ここまでの、軽い球が重い球及び壁と衝突した回数を数えると、円周率の $N$ 桁分の数字になるということです。出力結果は以下のようになります。徐々に円周率に近づいていますね！今回、衝突は全て弾性衝突なので、$e=1$ となります。連立方程式を $v', V'$ について解くと、王道の方法から非自明な方法、そしてクセのある方法がたくさんありましたね！どの方法たちも見ていて飽きることがなく、個人的には甲乙付け難かったです。。。以下が結果です。一般に乱数による値の誤差はサンプル数 $N$ に対して $1 / \sqrt{N}$ のオーダーで収束する(割と遅い)ので、なかなか精度は低いです。で求められます。半径 $r=1$ の場合を想定して、実装してみましょう！乱数を使うことで $\pi$ を求めることもできます。乱数とは、その名の通りランダムな数で、たとえばサイコロの出る目などが該当します。しかし、乱数を用いて $\pi$ を求めるには大量のサンプル($\sim$数百万)が必要で、人間が１００万回もサイコロを降るのはあまりにも退屈です。めんどくさいルーティンワークはコンピュータを導入して効率化しましょうということで、以下では乱数を用いて $\pi$ を求める方法を紹介していきます！次に、y = 2 * rand.rand(N) + np.sin(theta) で針の先端の $y$ 座標を計算します。2 * rand.rand(N) で針の中心の座標を、np.sin(theta) で針の中心と先端の高さの差を出します。針の先端が平行線と交わっているかどうかは [y>2] で判断します($y>2$ なら交わっている、そうでなければ交わっていない)。あとは y.size で交わっている針の本数を計算し、ごちゃごちゃいじれば終了です。多角形近似に比べると収束は遅いですが、等分数 $n$ の増加にしたがって値が円周率に収束しています。 LMDE 2 (Linux Mint Debian Edition 2; 64bit) での作業を想定。 2. そもそも円周率とは何かわかりますか？小学生くらいのときに教わるそうなのですが、私は忘れていました（笑）。 円周率は「円の円周の直径に対する比率」です。求めた方としては、円の円周を直径で割るだけです。 しかし、円の円周はどのように求めるのでしょうか？円の円周を正確に求めるのは難しいので、円周率を求めるアルゴリズムを使って円周率を求めます。 私はガウス＝ルジャンドルのアルゴリズムを使って求めてみました。この方法をpythonを使って解説していきたいと思います。 ちなみに …