Pythonで数値計算　decimalモジュール編｜Pythonで数学を学ぼう！第33回

計算誤差｜電子計算機の宿命

電子計算機の計算は万能ではありません。計算には誤差があります。電子計算機の資源が有限だからです。計算は正確・精確な結果を出力できなければ意味がありません。

電子計算機にとって計算誤差はマシンの誕生の時からの宿命ともいえます。電子計算機の発展は、計算誤差をいかに克服するかという難題に対するチャレンジとともにあったといえます。

連載第28回「最小2乗法」で紹介したように18世紀の数学者ガウスは、天体観測につきまとう測定誤差を科学することで誤差論という新しい数学をつくりだしました。

これと同じように計算の誤差をいかに克服するかという難題に対し、様々な工夫が考案されてきました。とくにソフトの世界では、計算精度を大きくするためのアルゴリズムの研究が進化しました。

浮動小数点数と単精度・倍精度

電子計算機では小数データは浮動小数点数（floating point number）という形式で表されます。浮動小数点数では数を仮数と指数に分けて表現します。

数＝仮数 × 10^指数

そしてメモリにデータを格納する際の国際規約（IEEE754）が単精度・倍精度・4倍精度です。単精度binary32はCPUが16ビット・32ビット時代の主流で、現在は64ビット時代なので倍精度binary64が主流です。

倍精度binary64は、符号ビットが1ビット、仮数部が52ビット、指数部が11ビットの計64ビットの構成です。これによって扱える数の範囲は

10^-307.95 〜 10^307.95

となります。

FPU（Floating Point Unit：浮動点小数演算処理装置）の発展により倍精度や拡張倍精度の計算が可能になりました。現在のx64と呼ばれる新しい64ビットCPUでは、FPUを使わずSSE2と呼ばれる拡張命令で倍精度実数の計算を行っています。

浮動点小数の詳細は2進数で考える必要があるので込み入ってきます。それよりも実際のPythonの計算精度を確認してもらう方が実感があるでしょう。

Pythonの計算精度「0.1問題」

Pythonで0.1+0.1+0.1と0.3を比較してみます。検証に用いる演算子==は非常に便利です。両辺が等しい場合にTrue、等しくない場合にFalseを返します。

0.1 + 0.1 == 0.2
True
0.1 + 0.1 + 0.1 == 0.3
False
0.1 + 0.1 + 0.1 +0.1 == 0.4
True

このように、0.1+0.1+0.1と0.3の場合だけFalseと判定されます。両辺の差を計算してみます。

0.1+0.1+0.1-0.3
5.551115123125783e-17

計算誤差があることがわかります。この状況の詳細はPython公式ドキュメントにあります。
15. 浮動小数点演算、その問題と制限
https://docs.python.org/ja/3/tutorial/floatingpoint.html
を参照してみてください。

decimalモジュール

このようなPythonの計算誤差に対する対処方法は2つあります。decimalモジュールと多倍長演算ライブラリmpmathです。まずはdecimalモジュールの使い方を見ていきます。

from decimal import Decimal
Decimal('0.1')+Decimal('0.1')+Decimal('0.1') == Decimal('0.3')
True

今度は結果がTrueとなりました。
このように、decimalモジュールからDecimal関数をimportして、0.1をDecimal('0.1')とすることで誤差のない小数として扱えるようになります。

2の平方根を100桁の精度で計算する

decimalモジュールの演算精度のデフォルトは28桁です。getcontext()関数を用いることで演算精度を指定できます。Decimal型数値には平方根、指数、対数などを計算するメソッドが備わっています。

2の平方根を有効数字100桁で計算してみます。

from decimal import Decimal, getcontext
getcontext().prec = 100
print(Decimal('2').sqrt())
1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573

まずgetcontext()関数もインポートして、getcontext().precに演算精度（有効数字）100を代入します。２をDecimal('2')としてsqrt()メソッドをつけます。Decimal('2').sqrt()のままだとDecimal('1.414213…7641573')（これをDecimal型数値という）のように出力されます。print()に渡すことで数値だけが出力されます。

e²を計算するにはexp()メソッドを使います。

print(Decimal('2').exp())
7.389056098930650227230427460575007813180315570551847324087127822522573796079057763384312485079121795

2の自然対数log 2 = ln 2を計算するにはln()メソッドを使います。

print(Decimal('2').ln())
0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875

2の常用対数log 2を計算するにはlog10()メソッドを使います。

print(Decimal('2').log10())
0.3010299956639811952137388947244930267681898814621085413104274611271081892744245094869272521181861720

このようにdecimalモジュールにより誤差がない正確な演算を行うことができます。しかし複雑な数式に対しては使うことができません。

次回は多倍長演算ライブラリmpmathによる計算をみていきます。