ラマヌジャンの公式で円周率を計算してみる

QuoraにGoogleの ”Jin Yamanaka” 氏が投稿していた、
ラマヌジャンの公式を使った円周率の計算法を自分でも動かしてみた。

\cdot \frac{1}{\pi}=\frac{2\sqrt{2}}{99^2}\sum_{n=0}^\infty \frac{(26390n+1103)\cdot(4n)!}{(4^n99^n\cdot n!)^4}

普通にmathライブラリをインポートしたpythonを使って計算すると、20桁程度の精度しかない。
しかもラマヌジャンの式は収束が速いようで、ほんの３ループほどで計算値の違いがなくなった。
これなら確かに、NASAが数値演算で"π"として扱っているのが３ループさせた値というのは本当だろう。
1 : 3.141592 7300133056603
2 : 3.141592653589793 878
3 : 3.1415926535897932385
4 : 3.1415926535897932385
5 : 3.1415926535897932385
6 : 3.1415926535897932385
7 : 3.1415926535897932385
8 : 3.1415926535897932385
9 : 3.1415926535897932385

９ループも実装してあるのに３ループ目から全く同じ数字というのもつまらないので、任意精度の結果が得られるという便利なmpmathライブラリを使って書き直してみた。
もちろん、まだ入れてない場合は追加でのインストールが必要。
精度の値(mpath.mp.prec)は、何度か計算させて十分なビット幅をカット＆トライで調整した値。256では大きすぎたので減らし、224=(128+64+32)という感じで適当に設定した。

import mpmath
mpmath.mp.prec = 224

def get_one_over_pie(total_n):
    v = 0
    for n in range(total_n):
        v += mpmath.factorial(4 * n) * (1103 + 26390 * n) / ((mpmath.factorial(n) ** 4) * (396 ** (4 * n) ))
    return mpmath.sqrt(8) / 9801 * v

for n in range(1, 10):
    print(n, ": ", 1 / get_one_over_pie(n))

1 : 3.141592 73001330566031399618902521551859958160711003355965653629013
2 : 3.141592653589793 87799890582630601309421664502932284887917396379151
3 : 3.14159265358979323846264 906570275889815667748046233478116839959564
4 : 3.1415926535897932384626433832795 5527315997421042037991121670389601
5 : 3.141592653589793238462643383279502884197 663818133030623976165591
6 : 3.14159265358979323846264338327950288419716939937 984683274351250728
7 : 3.14159265358979323846264338327950288419716939937510582 102093324247
8 : 3.141592653589793238462643383279502884197169399375105820974944592 76
9 : 3.14159265358979323846264338327950288419716939937510582097494459231

おそろしく収束が速い式だし、なんでこんな式が浮かんだのか想像もつかない。
天才って本当に居るんだな。