Pythonライブラリ mpmath 詳解：任意精度演算の世界へようこそ 🐍✨

1. 高精度計算

ライブラリの核となる機能です。`mp.dps` や `mp.prec` を設定することで、必要な桁数で計算を実行できます。これにより、標準の浮動小数点数では発生しうる桁落ちや丸め誤差の問題を回避できます。

from mpmath import mp, sin

# 標準 float での計算 (桁落ちの例)
x_float = 1e-8
result_float = (1 - (1 - x_float)**2) / x_float
print(f"Float result (x=1e-8): {result_float}") # 本来 2 に近い値になるはず

# mpmath での計算
mp.dps = 30
x_mpf = mp.mpf('1e-8')
result_mpf = (1 - (1 - x_mpf)**2) / x_mpf
print(f"mpmath result (x=1e-8, dps=30): {result_mpf}") # より正確な値

# 円周率を100桁計算
mp.dps = 100
print(f"Pi (100 digits): {mp.pi}")

# sin(1e-20) の計算
mp.dps = 50
small_angle = mp.mpf('1e-20')
print(f"sin(1e-20) with dps=50: {mp.sin(small_angle)}") # ほぼ 1e-20 に近い値
print(f"Approximation (1e-20): {small_angle}")
        

2. 広範な数学関数

`math` や `cmath` モジュールにある基本的な関数（`sqrt`, `exp`, `log`, `sin`, `cos`, `tan`, 三角関数の逆関数、双曲線関数など）はもちろん、特殊関数も豊富に実装されています。

• ガンマ関数 (`gamma`)、ゼータ関数 (`zeta`)
• 誤差関数 (`erf`, `erfc`)
• ベッセル関数 (`besselj`, `bessely`, `besseli`, `besselk`)
• エアリー関数 (`airyai`, `airybi`)
• 直交多項式 (ルジャンドル、チェビシェフ、エルミートなど)
• 超幾何関数 (`hyp0f1`, `hyp1f1`, `hyp2f1`, `hyper`)
• ランベルトのW関数 (`lambertw`)
• など多数

これらの関数もすべて任意精度で計算できます。

from mpmath import mp, gamma, zeta, erf, besselj

mp.dps = 30

print(f"Gamma(5.5): {mp.gamma(5.5)}")
print(f"Zeta(2) (should be pi^2/6): {mp.zeta(2)}")
print(f"pi^2/6: {mp.pi**2 / 6}")
print(f"erf(1): {mp.erf(1)}")
print(f"Bessel J0(2.5): {mp.besselj(0, 2.5)}")
        

3. 複素数演算

`mpc` 型を用いて、任意精度の複素数演算が可能です。実数と同様に、基本的な演算や数学関数が利用できます。

from mpmath import mp, mpc, exp, log, sqrt

mp.dps = 30

z1 = mpc(1, 1) # 1 + j
z2 = mpc('2+3j')

print(f"z1 + z2: {z1 + z2}")
print(f"z1 * z2: {z1 * z2}")
print(f"exp(z1): {mp.exp(z1)}")
print(f"log(z2): {mp.log(z2)}")
print(f"sqrt(-1): {mp.sqrt(mpc(-1))}") # sqrt(-1+0j)
        

4. 数値積分 (`quad`)

高精度な数値積分（定積分）を行うための `quad` 関数を提供しています。積分区間には無限大 (`mp.inf`, `-mp.inf`) も指定できます。様々なオプションで積分アルゴリズムを調整することも可能です。

from mpmath import mp, quad, exp, sin, inf

mp.dps = 50

# ∫[0, inf] exp(-x^2) dx = sqrt(pi)/2 の計算
integrand = lambda x: exp(-x**2)
result, error_est = mp.quad(integrand, [0, mp.inf], error=True)
print(f"Gaussian integral [0, inf]: {result}")
print(f"Expected sqrt(pi)/2: {mp.sqrt(mp.pi)/2}")
print(f"Estimated error: {error_est}")

# ∫[0, 1] sin(x)/x dx の計算 (Si(1))
integrand2 = lambda x: sin(x)/x
result2 = mp.quad(integrand2, [0, 1])
print(f"Integral of sin(x)/x from 0 to 1: {result2}")
        

5. 数値微分 (`diff`)

関数の数値微分を計算する `diff` 関数も用意されています。高階微分も計算可能です。

from mpmath import mp, diff, sin, cos, tan

mp.dps = 30
x_val = mp.pi / 4

# sin(x) の x=pi/4 での1階微分 (cos(pi/4))
d1 = mp.diff(mp.sin, x_val)
print(f"d/dx sin(x) at pi/4: {d1}")
print(f"cos(pi/4): {mp.cos(x_val)}")

# tan(x) の x=pi/4 での3階微分
d3 = mp.diff(mp.tan, x_val, n=3) # n=3 で3階微分
print(f"d^3/dx^3 tan(x) at pi/4: {d3}")
# 検算: tan'''(x) = 2 sec^2(x) (3 tan^2(x) + 1)
expected_d3 = 2 * (1/mp.cos(x_val)**2) * (3*mp.tan(x_val)**2 + 1)
print(f"Expected value: {expected_d3}")
        

6. 根探索 (`findroot`)

方程式 f(x) = 0 の根を高精度で求めるための `findroot` 関数があります。初期値を与えることで、ニュートン法などのアルゴリズムを用いて根を探索します。

from mpmath import mp, findroot, cos, sin

mp.dps = 50

# cos(x) = x の根を探す
# f(x) = cos(x) - x = 0 とする
f = lambda x: cos(x) - x

# 初期値 0.5 から探索
root = mp.findroot(f, 0.5)
print(f"Root of cos(x) = x: {root}")

# 検算: cos(root) と root がほぼ等しいか確認
print(f"cos(root): {mp.cos(root)}")
print(f"Difference |cos(root) - root|: {abs(mp.cos(root) - root)}")
        

7. 行列演算 (`matrix`)

`matrix` クラスを用いて、任意精度の要素を持つ行列を作成し、基本的な線形代数演算（和、差、積、転置、逆行列、行列式、固有値問題など）を行うことができます。ただし、大規模な行列計算には NumPy/SciPy の方が効率的な場合が多いです。mpmath の強みは要素を高精度で扱える点にあります。

from mpmath import mp, matrix

mp.dps = 25

A = mp.matrix([[1, mp.pi], [mp.e, 2]])
B = mp.matrix([[mpf('0.5'), 1], [0, mpf('1.5')]])

print("Matrix A:")
print(A)
print("Matrix B:")
print(B)

# 行列の和
print("A + B:")
print(A + B)

# 行列の積
print("A * B:")
print(A * B)

# 逆行列
try:
    A_inv = A**-1 # または A.inv()
    print("Inverse of A:")
    print(A_inv)
    # 検算 A * A_inv が単位行列になるか
    print("A * A_inv:")
    print(A * A_inv)
except Exception as e:
    print(f"Could not compute inverse of A: {e}")


# 行列式
det_A = mp.det(A)
print(f"Determinant of A: {det_A}")
        

8. 区間演算 (`iv`)

基本的な区間演算もサポートされています (`mpmath.iv`)。これは、数値計算における結果の範囲を厳密に評価したい場合に有用です。

from mpmath import iv, sin, pi

# 精度設定 (区間演算用)
iv.dps = 20

# 区間の作成
x_interval = iv.mpf([3.141, 3.142])
print(f"Interval x: {x_interval}")

# 区間に対する関数評価
sin_x_interval = iv.sin(x_interval)
print(f"sin(x) for x in {x_interval}: {sin_x_interval}")

# 円周率を含む区間
pi_interval = iv.pi()
print(f"Interval for pi: {pi_interval}")
print(f"sin(pi) interval: {iv.sin(pi_interval)}")
        

コンテキスト管理 (`workdps`, `workprec`)

計算の一部だけ一時的に精度を変更したい場合、`with` 文と `mp.workdps` または `mp.workprec` を使うと便利です。`with` ブロックを抜けると、精度は自動的に元に戻ります。これは、特定の関数呼び出しだけ高精度で行いたい場合などに有効です。

from mpmath import mp, sqrt, workdps

mp.dps = 15 # グローバル精度は15桁
print(f"Global dps: {mp.dps}")
print(f"sqrt(2) with global dps: {mp.sqrt(2)}")

# with ブロック内だけ精度を100桁に変更
with mp.workdps(100):
    print(f"  Inside with block, dps: {mp.dps}")
    sqrt2_high_prec = mp.sqrt(2)
    print(f"  sqrt(2) with dps=100: {sqrt2_high_prec}")

# with ブロックを抜けると精度は元に戻る
print(f"Outside with block, dps: {mp.dps}")
print(f"sqrt(2) with global dps again: {mp.sqrt(2)}")

# extradps/extraprec で現在の精度に加算することも可能
with mp.extraprec(50): # 現在の精度 + 50ビット
     print(f"  Inside extraprec(50), prec: {mp.prec}")
     print(f"  sqrt(2) with extra 50 bits: {mp.sqrt(2)}")
              

パフォーマンスに関する考慮事項

任意精度計算は、ハードウェアレベルで最適化された標準の倍精度浮動小数点演算に比べて、必然的に計算コストが高くなります。精度 (`dps` や `prec`) を高く設定すればするほど、計算時間とメモリ使用量は増加します。

• 必要最小限の精度: 計算に必要な精度を見極め、過度に高い精度を設定しないようにしましょう。
• gmpy2 の利用: 可能であれば gmpy2 をインストールすることで、特に高精度域でのパフォーマンスが大幅に向上します。
• アルゴリズムの選択: mpmath は多くの関数で高精度に適した効率的なアルゴリズムを採用していますが、問題によっては計算量が増大することもあります。
• NumPy/SciPy との使い分け: 大規模な配列や行列に対する標準的な数値計算（精度がそれほど要求されない場合）は、NumPy や SciPy の方が高速です。高精度が必要な部分のみ mpmath を使う、といった使い分けが有効です。

他のライブラリとの連携

mpmath は、特に数式処理ライブラリ SymPy と密接に連携しています。SymPy は内部で mpmath を利用して、シンボリックな式を高精度で数値評価します。

import sympy
from sympy import symbols, integrate, N

# SymPy でシンボリックな式を定義
x = symbols('x')
expr = sympy.exp(-x**2)

# SymPy の integrate を使ってシンボリックに積分
integral_expr = integrate(expr, (x, 0, sympy.oo))
print(f"Symbolic integral: {integral_expr}")

# N() を使って mpmath で高精度数値評価 (ここでは 50 桁)
numeric_value = N(integral_expr, 50)
print(f"Numeric value (dps=50): {numeric_value}")

# mpmath を直接使った場合と比較
from mpmath import mp, quad, exp, inf, sqrt, pi
mp.dps = 50
mp_value = quad(lambda t: exp(-t**2), [0, inf])
print(f"mpmath direct value (dps=50): {mp_value}")
print(f"sqrt(pi)/2: {sqrt(pi)/2}")
              

NumPy 配列と mpmath オブジェクトを直接組み合わせることは通常できませんが、リスト内包表記や `numpy.vectorize` などを使って、NumPy 配列の各要素に mpmath の関数を適用することは可能です。ただし、パフォーマンスには注意が必要です。