Python/Numpyで任意精度（arbitrary-precision）、多倍長精度（multiple-precision）の演算を行う。

準備

mpmathとgmpy2をインストールします。

pip install gmpy2
pip install mpmath

mpmathは単体でもPython用の任意精度演算ライブラリとして動作します。デフォルトでは内部でPythonの組み込みint型を使用しますがgmpy2がインストールされている場合はこれを自動的に検出して内部で使用します。ドキュメントによるとgmpy2をバックエンドとして使用することで特に100桁以上を扱う場合において計算を大幅に高速化できます。これを使わない手は無いのでmpmath+gmpy2は実質セットで使われると考えてよいでしょう。

次を確認しましょう。

>>> import mpmath
>>> mpmath.libmp.BACKEND
'gmpy'

gmpyと出ていればgmpy2がmpmathによって正しく検出されています。

使用例

from mpmath import mp

mp.dps = 100

x = mp.mpf('1234.56789123456789123456789123456789')
y = mp.mpf('9876.54321987654321987654321987654321')
print(x + y)
print(x * y)
print(x / y)

これで10進で約100桁の精度で演算が行えます。精度を指定するにはmp.dpsもしくはmp.precに整数を代入します。mp.dpsは10進の桁数を指定できます。mp.precには浮動小数点の仮数部のビット数を指定できます。64bit浮動小数点の仮数部は53bitで、これは10進約15桁にあたります。

どちらも大きい数を指定するほど精度は高くなりますが、計算機による計算であることには変わりないので少しの誤差は出ます。

mp.mpfで実数、mp.mpcで複素数、mp.matrixで行列を表すことができます。値の表現には文字列を用いることに注意が必要です。数値で値を表現すると一度組み込み型の精度に丸められたのちにmp.mpfに食わせることになるので精度が落ちます。例をみてわかる通り、普通の四則演算は特に意識せずに行うことが出来ます。

mp.piや、mp.sin、mp.cosなどのnumpyにある数学の定数や関数の任意精度版は普通に用意されています。必要に応じて調べていけば特に困ることは無いでしょう。

numpy配列

mpmathは自身でもmp.matrixという2次元配列型を持っていますが、配列を扱う上ではnumpyのn次元配列型であるnumpy.ndarrayを使う方が便利なことがあります。mp.mpfなどには四則演算が定義されているので、dtype=objectのnumpy.ndarrayを使えば強引に普通のnumpy配列と同じよな操作で任意精度演算が実現可能です。

import numpy as np
import matplotlib.pyplot as plt

from mpmath import mp

N = 10

A = np.array([[
    mp.mpf('0.9')**abs(i-j) for i in range(N)
] for j in range(N)])

b = np.arange(N)*mp.mpf('0.1')/N

print(A @ b)
print(A @ A @ b)
print(A @ A @ A @ b)

np.linspaceなどは動きませんが、broadcastを使って同じようなことが出来ています。scipy.linalgにある関数などをはじめ多くの関数はdtype=objectの任意精度型配列に対して動きません。mpmathにはmp.matrixに対する直接的なscipy.linalg関数の代替が用意されているので必要に応じてこれらを使う必要があります。

mp.matrixは辞書型を用いたスパース表現によって実装されています。当然のように任意精度型は普通の浮動小数点型よりもサイズが大きいため、2次元配列の要素全てをメモリ上に置いているとサイズが大きく計算も非効率になります。行列のスパース表現としては、CRS (Compressed row storage)などのデータ形式がありますがドキュメントを見る限りmp.matrixは辞書型を使うようです。任意精度型の素の演算が普通より重いこともあり、効率的な計算のためには必要に応じて任意精度型の演算と普通のnumpy.ndarrayによる行列計算を明示的に切り替えることが必要でしょう。

matplotlibへの入力

matplotlibにdtype=objectのnumpy.ndarrayを正しいサイズで渡せばプロットされます。

import numpy as np
import matplotlib.pyplot as plt

from mpmath import mp

N = 100
x = np.arange(N)*2*mp.pi/N
y = np.array([mp.sin(xi) for xi in x])

plt.plot(x, y)
plt.show()

mp.sinはnp.sinと違ってarray-likeを喰わせても要素ごとに関数の入力として扱ってくれません。なので、この例のようにリスト内包表記などを使って1個1個の要素に対してmp.sinを当てていく必要があります。一方、matplotlibの関数には正しいサイズのarray-likeを渡すと任意精度型でも正しくプロットしてくれています。

matplotlibは可視化に用いるため、直前で任意精度演算から組み込み型に変換されることは問題ありません。もしmatplotlibに直接渡せない場合は明示的に直前で組み込み型に変換する必要がありましたが、この操作は自動手的に行われるようです。