概要

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

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のアルゴリズムの実装

完成版のプログラムは以下のものです。

chudnovsky_pi.py

# 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

他の公式でも計算して、計算結果を検証する。

参考

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

円周率計算の高速化

decimal --- 十進固定及び浮動小数点数の算術演算

Pythonのdecimalモジュールを使用した円周率1億桁の計算

概要

目次