mpmathのnsumを使って無限数列の和を高精度で求める

pythonで使えるmpmathモジュールのnsum関数を使ってみたので紹介します。

【例題】次の無限数列の和を小数点以下14桁まで求めよ

1.  \quad \sum_{n=0}^{\infty}\frac{1}{n!} \\
2.  \quad \sum_{n=1}^{\infty}\frac{1}{n^2}

単純ループプログラム

$1.$に関しては以下のように定義通りに$n=25$までループさせると正解（=eの値）が求まります。

from sympy import factorial, E
print(f"{sum([float(1/factorial(n)) for n in range(25)]) :.14f}")
print(E.n())
#2.71828182845905
#2.71828182845905

しかし$2.$に関しては$10^6$までループさせても正解($=\zeta(2) $)にはなりません。

from sympy import zeta
print(f"{sum([float(1/n**2) for n in range(1,10**6)]) :.14f}")
print(f"{zeta(2).n():.14f}")
# 1.64493306684777
# 1.64493406684823

mpmathのnsum

そこでmpmathのnsumの出番にとなります。この関数の詳細は、以下のリンクにありますが、ポイントはmp.dpsで求める値の精度を指定できることです。

mpmath: Sums, products, limits and extrapolation

今回はmp.dps=100として100桁まで求めてみました。例題の$1.$も$2.$も正確に求まっているのが分かると思います。

from mpmath import *
from sympy import N, E, zeta
prec = 100
mp.dps = prec
print(nsum(lambda n: 1/fac(n), [0, inf]))
print(N(E,prec))
print(nsum(lambda n: 1/n**2, [1, inf]))
print(N(zeta(2),prec))
#2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427
#2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427
#1.644934066848226436472415166646025189218949901206798437735558229370007470403200873833628900619758705
#1.644934066848226436472415166646025189218949901206798437735558229370007470403200873833628900619758705

（開発環境：Google Colab）

この考え方はProject Euler Problem 722: Slowly Converging Seriesを解くのに役に立ちます。