概要
PythonのdecimalモジュールでChudnovskyのアルゴリズム
\frac{1}{\pi} = 12 \sum_{n=0}^\infty \frac{(-1)^n (6n)! (13591409 + 545140134 n)}{(3n)!(n!)^3 640320^{3n + 3/2}}
を実装して円周率を1億桁計算しました(Qiita初投稿です。間違い等がありましたらコメントでご指摘お願いします)。
アルゴリズムの詳細についてはperiaさんの記事(Chudnovsky の公式を用いた円周率の計算用メモや円周率計算の高速化)で解説されているため、本記事では主にdecimalモジュール固有の内容を扱います。
プログラムは以下の環境で実行しています。
- OS: Arch Linux
- CPU: Intel(R) Core(TM) i5-8400
- Memory: 16GB
- Python: 3.9.6
目次
- decimalモジュールで巨大な桁数を高速に計算するために
- pydecimal vs cdecimal
- 指数部の設定
- 平方根関数の実装
- Chudnovskyのアルゴリズムの実装
- 計算時間
- TODO
- 参考
decimalモジュールで巨大な桁数を高速に計算するために
pydecimal vs cdecimal
pythonの多倍長浮動小数点数ライブラリであるdecimalには2つの実装があります。
- pydecimal (Pythonで書かれた実装)
- cdecimal (C言語で書かれた実装で、Python 3.3から標準ライブラリとして含まれる。大きな桁数の乗算を数論変換で行うため高速)
Python (CPython) 3.3以降では、decimalモジュールをimportすると自動的にcdecimalを使うことができますが、古いバージョンやPyPyにはcdecimalが含まれていません。数万桁以上の計算を現実的な時間で行うにはcdecimalが不可欠なため、本記事のプログラムはcdecimalを含まないPythonインタープリタでの実行をサポートしていません。
指数部の指定
decimalモジュールで扱える数の指数部は、デフォルトでは-999999から999999に制限されています。
from decimal import Decimal
Decimal("10") ** 100000000
# Traceback (most recent call last):
# File "<stdin>", line 1, in <module>
# decimal.Overflow: [<class 'decimal.Overflow'>]
Decimal("10") ** -100000000
# Decimal('0E-1000026')
これを回避して巨大な桁数を計算するには、getcontext().Emaxとgetcontext().Eminを設定します。
from decimal import Decimal, getcontext, MAX_EMAX, MIN_EMIN
getcontext().Emax = MAX_EMAX # MAX_EMAX = 999999999999999999
getcontext().Emin = MIN_EMIN # MIN_EMIN = -999999999999999999
Decimal("10") ** 100000000
# Decimal('1.000000000000000000000000000E+100000000')
Decimal("10") ** -100000000
# Decimal('1E-100000000')
平方根関数の実装
decimalモジュールでは、大きな桁数の乗算に数論変換を使用しているため乗算は高速です。一方、平方根を求めるDecimal.sqrt()メソッドは低速なため、Chudnovskyのアルゴリズムで大きな桁数の円周率を計算するには、自分で平方根関数を実装しなければいけません。平方根関数はNewton-Raphson法で実装します。
def my_sqrt(x):
x2 = Decimal(x) # 入力としてintを受け取るため
if getcontext().prec <= 128: # 桁数が小さいときはDecimal.sqrt()を使用
return x2.sqrt()
else:
prec_orig = getcontext().prec
t = (prec_orig * 53) // 100 # 理論上は50でよいが丸め誤差対策のため少し多めに
prec_list = []
while t > 128:
prec_list.append(t)
t = t // 2
prec_list.reverse()
getcontext().prec = 128
a = 1 / x2.sqrt() # 初期値はDecimal.sqrt()で生成
for i in prec_list:
getcontext().prec = i + 10
a = a + (a * (Decimal(1) - x2 * (a * a)) / Decimal(2))
getcontext().prec = prec_orig
a = a + (a * (Decimal(1) - x2 * (a * a)) / Decimal(2))
return x2 * a
Decimal.sqrt()と自作の平方根関数で$\sqrt{2}$の計算時間を計測すると、
桁数 | Decimal.sqrt()の計算時間(s) | 自作関数の計算時間(s) |
---|---|---|
1万 | 0.0390 | 0.00387 |
10万 | 0.729 | 0.0138 |
100万 | 10.2 | 0.139 |
1000万 | - (> 120) | 2.04 |
1億 | - | 18.7 |
となり、大幅に高速化できました。
Chudnovskyのアルゴリズムの実装
完成版のプログラムは以下のものです。
# 1億桁計算するには、"python3 chudnovsky_pi.py 100000000" を実行
from decimal import Decimal, getcontext, MAX_EMAX, MIN_EMIN
import math, time
from sys import argv
getcontext().Emax = MAX_EMAX # 指数部の指定
getcontext().Emin = MIN_EMIN # 指数部の指定
# cdecimalが含まれるか確認
def confirm_cdecimal():
a = True
try:
import _decimal
except ModuleNotFoundError:
a = False
return a
# 平方根関数
def my_sqrt(x):
x2 = Decimal(x)
if getcontext().prec <= 128:
return x2.sqrt()
else:
prec_orig = getcontext().prec
t = (prec_orig * 53) // 100
prec_list = []
while t > 128:
prec_list.append(t)
t = t // 2
prec_list.reverse()
getcontext().prec = 128
a = 1 / x2.sqrt()
for i in prec_list:
getcontext().prec = i + 10
a = a + (a * (Decimal(1) - x2 * (a * a)) / Decimal(2))
getcontext().prec = prec_orig
a = a + (a * (Decimal(1) - x2 * (a * a)) / Decimal(2))
return x2 * a
# Binary Splitting
def computePQT(a, b):
if b == a + 1:
p_int = (1 - 2*b) * (6*b - 5) * (6*b - 1)
q_int = (b**3) * 10939058860032000
t_int = p_int * (13591409 + 545140134*b)
return [Decimal(p_int), Decimal(q_int), Decimal(t_int)]
else:
m = (a + b) // 2
[pam, qam, tam] = computePQT(a, m)
[pmb, qmb, tmb] = computePQT(m, b)
p = pam * pmb
q = qam * qmb
t = tam * qmb + pam * tmb
return [p, q, t]
# 円周率計算
def chudnovsky_pi(n):
getcontext().prec = n
terms = math.ceil(n / 14.181647462) + 10
[_, q, t] = computePQT(0, terms)
num = q * Decimal(426880) * my_sqrt(10005)
den = t + q * Decimal(13591409)
return num / den
def main():
if not confirm_cdecimal(): # cdecimalが含まれない場合はプログラムを停止
print("Decimal (cdecimal) module is not found.")
else:
prec = int(argv[1]) + 10
print("This program might take much time for large digits. Press Ctrl + C to quit.")
t1 = time.time()
pi = chudnovsky_pi(prec)
file_name = "pi_" + argv[1] + "digits.txt"
f = open(file_name, "w")
f.write(str(pi)[0:int(argv[1])+2])
f.close()
t2 = time.time()
print("The result is written in " + file_name + ".")
print("Elapsed time: " + str(t2 - t1) + " seconds.")
if __name__ == "__main__":
main()
計算時間
1万桁から1億桁の計算時間は以下の表のようになりました(ファイル出力の時間も含みます)。
桁数 | 計算時間(s) |
---|---|
1万 | 0.00840 |
10万 | 0.181 |
100万 | 2.84 |
1000万 | 48.6 |
1億 | 644 |
円周率1億桁を11分弱で計算できました。
TODO
他の公式でも計算して、計算結果を検証する。