円周率を計算した

円周率を求めるプログラムを作りました。
目標としてはとりあえず100万桁、最終的には1億桁くらい計算したいです。
【環境】
PC:MacBook Air(2018年モデル)
CPU:i5-8210Y
RAM:8GB
Python 3.7.0
Anaconda 3 ,Jupyter 5.6.0

円周率の計算方法

とりあえず、Wikipedia先生¹によるとガウス＝ルジャンドルのアルゴリズムというのが速いそうです。
↓計算方法(入力が面倒だったのでWikipediaをスクショコピペ)

何も考えず実装

小数第n位まで求める時はlog₂n回くらい反復を行えばいいらしいので、log₂10⁶≒19.9にマージンつけて21回計算させてみました。

gl1.py

import math

def updater(pilist):
    a_n ,b_n, t_n, p_n = pilist
    a_n1 = (a_n + b_n) /2
    b_n1 = math.sqrt(a_n * b_n)
    t_n1 = t_n - p_n * (a_n - a_n1)**2
    p_n1 = 2 * p_n
    return a_n1, b_n1, t_n1, p_n1

a,b,t,p = 1,1/math.sqrt(2),1/4,1
for i in range (21):
    a,b,t,p = updater((a,b,t,p))

pi = ((a + b)**2)/(4 * t)

print(pi)

結果：3.141592653589794

知ってた。

Decimalを使って桁数を増やす

前項では浮動小数点の桁数を何も考えていなかったので、小数点以下15位までしか計算(表示)できませんでした。当然のことです。
そこで、(十進法の世界において)任意の桁数で正確な計算ができるDecimalを使いました。

Decimalを使った実装は以下の通りです。上のコードにDecimalとかを付け加えたりするだけです。
とりあえず10000桁の計算です。

gl2.py

import math
from decimal import *
dig = 10000 #計算させたい桁数
N = math.ceil(math.log2(dig))+1 #反復する回数(自動計算)
getcontext().prec = dig + 1 #10000桁には整数部分(3)も含まれているため、1を足す
digset = Context(prec = dig + 1)

def updater(pilist):
    a_n ,b_n, t_n, p_n = pilist
    a_n1 = (a_n + b_n) /2
    b_n1 = Decimal(a_n * b_n).sqrt()
    t_n1 = t_n - p_n * digset.power((a_n - a_n1),2)
    p_n1 = 2 * p_n
    return a_n1, b_n1, t_n1, p_n1

a,b,t,p = Decimal(1),Decimal(2).sqrt()/2,Decimal(1)/4,Decimal(1)
for i in range (N):
    a,b,t,p = updater((a,b,t,p))

pi = (digset.power((a + b),2))/(4 * t)
print(pi)

結果：
3.141592653589793238462643383279502884197169399375105820974……
(中略)
……64836999892256959688159205600101655256375676

実際の円周率は3.14……64836999892256959688159205600101655256375678のため、最後の1桁だけズレてますがそれ以外は合っています。(最終桁の丸めの問題？数桁だけ余分に計算したらこの問題は解消できそうです)

問題点

遅い。とにかく遅い。
桁数と処理時間の関係は以下の通り。

100万桁計算するのに278秒かかりました。
主犯は100万桁の桁数を持つ数字の根号を取るところな気がするので、アルゴリズムを変えて計算してみます。
(同じアルゴリズムを使っているスーパーπはもっと早いのでPython自体の速度も含め、別のところにも原因はあると思いますが。。。)

今度は高速なアルゴリズムらしいChudnovskyの公式を使ってやってみようと思います。

Chudnovskyの公式を用いた円周率の計算

このブログ²を参考にしました。
以下の式で円周率を求めることができる、らしいです。

\frac{1}{\pi} = 12 \sum_{n=0}^\infty \frac{(-1)^n (6n)! (A + B n)}{(3n)!(n!)^3 C^{3n + 3/2}} \qquad \left(A = 13591409, B = 545140134 , C = 640320\right)

式変形を繰り返すことにより、以下の漸化式を計算していくことで円周率を算出できるようになる、らしいです。

π 〜 \frac{C^\frac{3}{2} Q_N}{12(T_N+A Q_N)} \\

ここで、

P_N =P_{N-1} × p_N \\
Q_N = Q_{N-1} × q_N \\
T_N = T_{N-1} × q_N + a_N × P_N \\
a_N = (-1)^N(A+B×N)\\
q_N = \frac{N^3C^3}{24}\\
p_N = (2N-1)(6N-1)(6N-5)\\
P_0 = Q_0 = 1, T_0 = 0

この漸化式のNを1つ進めるたびに14桁くらい円周率に近づいていくらしいです。
pythonで実装すると以下のようになります。

chudnovsky.py

from decimal import *
import math
dig = 10000
getcontext().prec = dig + 1 #10000桁には整数部分(3)も含まれているため、1を足す
digset = Context(prec = dig + 1)

A = 13591409
B = 545140134
C = 640320
C3_24 = int(C**3/24)

P = 1
Q = 1
T = 0
N = math.ceil(dig/14)+10

def update(N ,P, Q, T):
    for i in range(1,N):
        a = A+B*i if i%2 == 0 else -A-B*i
        q = i**3 * C3_24
        p = (2*i-1)*(6*i-1)*(6*i-5)
        P_new = P * p
        Q_new = Q * q
        T_new = T * q + a * P_new
        P, Q, T =  P_new, Q_new, T_new
    return P,Q,T

P,Q,T = update(N,P,Q,T)
pi = digset.power(C,Decimal("1.5")) /12 * Q / (T + A * Q)

↓計算に要する時間(ガウス＝ルジャンドル法といっしょに)

遅すぎる・・・
ガウス＝ルジャンドル法との比較をするために両方１００万桁まで計算したかったのですが、MacBook Airが唸り声をあげたので3万桁までで断念しました。

参考にしたブログ²では、Core 2 Duo P8600のパソコンでこの計算を0.0数秒〜数秒でこなしています。
アルゴリズムはほとんど変わらず、計算能力も劣っていないはずなのにここまで遅いのは何故でしょうか。
今のところ、Python自体の計算が遅いことに起因してるのかなと思っています(ブログはC)。

改善案

【Cythonを用いてみる】
①そのままCythonを使って動かす
　⇨ほとんど改善せず(数割早くなっている、らしい)
②変数の型を宣言して使ってみる
　⇨PやQをintとして宣言するとあっという間にオーバーフローして終了(16?32?bitの数字しか使えない)
③Cythonで使える任意精度計算ライブラリを導入する
　⇨やり方がさっぱりわからない

結論

Pythonは円周率計算に向かない

参考文献

[1] Wikipedia 円周率
https://ja.wikipedia.org/wiki/%E5%86%86%E5%91%A8%E7%8E%87
[2]円周率を1億桁計算しました！ ― その試行錯誤の詳しい経緯と結果 ー
https://itchyny.hatenablog.com/entry/20120304/1330870932