6
6

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

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

Posted at

概要

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

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

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 --- 十進固定及び浮動小数点数の算術演算

6
6
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
6
6

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?