円周率の日

円周率の分数近似（4年前）からの続き

3月14日は円周率の日だ。

PC-9801を導入して最初のころにBASICで円周率の計算の例題を入力して試してみた。参考書にあったのがマチンの公式だった。$4 \tan^{-1 }\frac{1}{5}- \tan^{-1}\frac{1}{239}=\frac{\pi}{4}$を級数展開して求めるもので，$\tan^{-1} 1 = \frac{\pi}{4}$より収束が良く計算時間も稼げる。それでも当時のBASICインタープリタでは1000桁の計算に多少時間がかかっていた。今では，MathematicaでN[Pi,1000000]とするだけで，0.28秒でMacbook Air M1 の画面に100万桁の円周率が表示される。

1980年代，まだ日本のベクトル型スーパーコンピュータに勢いがあった頃，東大大型計算機センターを舞台として，金田康正さんが円周率計算桁数の世界記録を更新し続けていた。かれらのグループの最後の輝きは，2002年にHITACHI SR8000により1兆桁を計算したというものだった。

その後，時代はPCによる計算競争に移ったが，現在の記録はやはり日本人が持っている。Google Cloudの技術者である岩尾エマはるかが，2019年に31兆4159億2653万5897桁を記録し，2022年には100兆桁に達した。チュドノフスキーの公式を用いた計算はGoogle Cloudで157日かけて実行された。

円周率の計算は，arctan型の何種類かの公式やガウス=ルジャンドルのアルゴリズムやラマヌジャン型公式が用いられることが多い。チュドヌフスキーの公式は次のようなものだ。$n$が一つ進むごとに14桁近く精度が増えていく。

$ \dfrac{1}{\pi}=12 \Sigma_{n=0}^{\infty} \dfrac{(-)^n (6n)!(545140134n+13591409)}{(3n)!(n!)^3 640320^{n+3/2}}$

これをさらにすすめたものに，ボールウェィンの公式があったので確かめたところ，Juliaの計算値で$n>=1$の場合がおかしい。公式が間違って印刷されているのかと思ったが，Mathematicaではうまくいく。結局，Juliaの変数の任意精度の取り扱いが間違っていた。こちらのほうは，$n$が一つ進むごとに24桁近く精度が増えていく。

$A=1657145277365+212175710912\sqrt{61}$

$B=107578229802750+13773980892672\sqrt{61}$

$C=(5280\,(236674 + 30303\sqrt{61}))^3$

$ \dfrac{1}{\pi}=12 \Sigma_{n=0}^{\infty} \dfrac{(-)^n (6n)!(A+B n)}{(3n)!(n!)^3 C^{n+1/2}}$

以下が問題のJuliaのコードである。

function pit(m)

s61 = sqrt(big(61))

A = 1657145277365 + 212175710912*s61

B = 107578229802750 + 13773980892672*s61

C = (5280*(236674 + 30303*s61))^3

#A = BigFloat(13591409)

#B = BigFloat(545140134)

#C = BigFloat((640320)^3)

sum = big(0)

    for n in 0:m

        n6=factorial(big(6*n))

        n2=factorial(big(n))^3

        n3=factorial(big(3*n))

        sum=sum+(-1)^n*n6/(n2*n3)*(A+B*n)/C^(n+1/2)

    end

return BigFloat(1/(12*sum))

end

bpi = big(pi)

println(bpi-pit(0))

println(bpi-pit(1))

println(bpi-pit(2))

println(bpi-pit(3))