pi

円周率って不思議だな…と思ったので計算してみた

ex_1.py

• 色々調べる前に何も調べず思いつきでやってみる
• 正2^(n+1)角形の1辺の長さを計算して円周率の近似を求めています
• これはアルキメデスの手法という名前がついているらしい
• 1 - sqrt(a) の部分はnが大きくなると誤差が大きくなるので式変形した
• pythonのfloatでは十数桁の精度が限界

ex_2.py

• floatでは限界が来たのでDecimalを使ってみる
• 10^4回のループでやっと1000桁求まる（とても遅い）

ex_3.py

• アルキメデスの方法では限界がきたので別の方法を調べてみた
• マチンの公式を使っている（めちゃくちゃ中途半端な分数なのが不思議）
• 10^4桁の計算に30秒ほどかかる

ex_4.py

• マチンの公式でも限界がきたのでガウス＝ルジャンドルのアルゴリズムを使う
• bの値はnが大きくなるにつれて収束していくのでニュートン法の初期値に使える
• 10^6桁の計算に60秒ほどかかる

ex_5.py

• チュドノフスキー・アルゴリズムを使ってみる（見た目がすごい）
• ほとんど整数同士の計算で最後だけ分数の計算があるのがポイント
• pythonのintに桁数の制限がないのを初めて知った
• 分母分子の計算はループよりも2分法で計算したほうが高速
• 最後の分数部分でDecimalが遅かったのですべてintで計算を行っています
• 10^6桁の計算に40秒ほどかかる（ガウス＝ルジャンドルのアルゴリズムより少し速い）

ex_6.py

• Pythonの標準機能だけではそろそろ限界になってきたのでgmpy2を使う
• 基本的にはintをmpzに置き換えて最後の分数の計算をmpfrで行うようにした
• 10^7桁の計算に6秒ほどかかる
• 円周率のテキストデータと照らし合わせて検算する方法がそろそろ厳しくなってきた

ex_7.py

• ここからさらに並列化して10^8桁を狙ってみる
• 並列処理で求めた分母分子を1つにまとめるところも並列化しました
• jobsの値は8ぐらいがちょうど良さそう
• 求めた分母分子からpiを求める部分がボトルネックになっている
• 10^8桁の計算に40秒ほどかかる（未検算、500万桁まではあってる）
• これ以上速くするならもうCで書いたほうがいい…？

ex_8.py

• Prime Factorization Optimizationに挑戦してみる
• その前に、まずは素数に慣れるために素数を計算してみる
• エラトステネスの篩を使って計算している（もっと速い方法もあるらしい）
• 10^9までの素数をすべて計算するのに30秒ほどかかる

TODO

• BBPアルゴリズムを使って検算できるようにする
• 10^9桁を目指してみたい（ここまでくると円周率のテキストデータだけで1GB近くなる）
• timeitが引数のデフォルト値を指定した場合に対応していないのをどうにかする