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

# 概要

Pythonの[decimalモジュール](https://docs.python.org/ja/3/library/decimal.html)でChudnovskyのアルゴリズム

```math
\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 の公式を用いた円周率の計算用メモ](https://qiita.com/peria/items/97af2df9845decc108ae)や[円周率計算の高速化](https://qiita.com/peria/items/97af2df9845decc108ae))で解説されているため、本記事では主に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に制限されています。

```python
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を設定します。

```python
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法](https://qiita.com/peria/items/97af2df9845decc108ae)で実装します。

```python
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 の公式を用いた円周率の計算用メモ](https://qiita.com/peria/items/97af2df9845decc108ae)

[円周率計算の高速化](https://qiita.com/peria/items/97af2df9845decc108ae)

[decimal --- 十進固定及び浮動小数点数の算術演算](https://docs.python.org/ja/3/library/decimal.html)