目次
Rumpの例題
1988年、IBMドイツ基礎研究所のSiegfried M. Rumpは電子計算機による計算誤差を示すための興味深い数値計算問題を考えた。
これをPythonで計算してみた結果は次の通り。
>>> a=77617
>>> b=33096
>>> c=333.75*b**6+a**2*(11*a**2*b**2-b**6-121*b**4-2)+5.5*b**8+a/(2*b)
>>> print(c)
1.1805916207174113e+21
このPythonによる計算結果はe+21（10の21乗）であることから、Rumpの論文にある結果1.17…とは全く異なる数値です。
真値は-0.827396…ですから、通常の電子計算機による計算結果は符号も値も異なるという非常に興味深い結果であることがわかります。電子計算機だから正しいというのは神話であるということです。
本連載「キャサリン・ジョンソン｜軌道計算でNASA宇宙計画を支えた黒人女性数学者」でも紹介しましたが、1962年の宇宙船フレンドシップ7による軌道計算について宇宙飛行士ジョン・グレンが電子計算機の計算を信用できず「キャサリンに計算させない限り飛ばない」と言った逸話が思い出されます。
電子計算機による計算の失敗事例 その1｜湾岸戦争（1991年）
1962年の話は昔話だろうと思いきや、実はそうではありません。電子計算機の計算誤差による深刻な事件はその後も続きました。
1991年の湾岸戦争の最中、米軍兵舎を狙って発射されたスカッドミサイルをパトリオットミサイルによって迎撃するのが失敗に終わり兵士28名の命が失われました。
システム起動からの経過時間計測の計算に計算誤差が影響したことが原因とされています。前回紹介した「Pythonの計算精度「0.1問題」」です。10進数の0.1は、2進数では0.0001100110011…となり有限桁では表現できません。
電子計算機による計算の失敗事例 その2｜アリアン5ロケット事故（1996年）
ESA（欧州宇宙機関）のアリアン5ロケットは10年近くの歳月と70億ドルの費用を投じて開発されました。1996年にロケットは離陸直後に爆発し打ち上げは失敗に終わりました。
原因は、慣性座標システムのソフトウェアのエラー。ロケットの水平速度に関連する64ビット浮動小数点を16ビット符号付整数へ変換する際、数字が16ビット符号付整数として保存できる最大値の32,768を超えてしまい、エラーを引き起こしてしまったということです。
精度保証付き数値計算｜decimalモジュールの場合
そこで必要になるのが「精度保証付き数値計算」と呼ばれる技術です。前回紹介したdecimalモジュールはその1つです。さっそくRumpの例題を計算してみます。
>>> from decimal import Decimal, getcontext
>>> getcontext().prec = 50
>>> a = Decimal('77617')
>>> b = Decimal('33096')
>>> c = Decimal('333.75')*b**6+a**2*(11*a**2*b**2-b**6-121*b**4-2)+Decimal('5.5')*b**8+a/(2*b)
>>> print(c)
-0.8273960599468213681411650954798162919990331157844
真値が得られました。
ところがdecimalモジュールによる任意精度演算は高度な技術計算には不向きです。sin(0.1)を計算させてみます。
>>> from decimal import Decimal, getcontext
>>> getcontext().prec = 50
>>> import math
>>> math.sin(0.1)
0.09983341664682815
>>> math.sin(Decimal('0.1'))
0.09983341664682815
sin(0.1)の真値は
0.09983341664682815230681419
ですから小数点以下17桁は正しい結果です。しかしmath.sin(Decimal('0.1'))としても精度を50桁に指定したのにもかかわらず、小数点以下17桁しか得られません。
精度保証付き数値計算｜mpmathモジュールの場合
mpmathはPythonの任意精度浮動小数点演算パッケージです。sin(0.1)を50桁の精度で計算するには次のようにします。
>>> from mpmath import mp
>>> mp.dps = 50
>>> mp.pretty = True
>>> mp.sin(mp.mpf('0.1'))
0.099833416646828152306814198410622026989915388017982
2行目のmp.dpsは10進精度です。mp.precで2進精度を指定できます。
3行目のmp.prettyはデフォルトがFalseです。Falseだとprint()に渡してあげないと、きれいな数値出力になりません。mp.pretty = Trueとすることでprint()に渡すことなくきれいな数値出力が得られます。
Pythonで三角関数を扱うにはいくつかのモジュールが必要です。ここではmp.sin()が使えます。そこに0.1をmp.mpf('0.1')として渡すことで任意精度浮動小数点演算が実行されます。
sin(0.1)の真値（小数点以下50桁）は
0.099833416646828152306814198410622026989915388017982
ですから、ジャストな値であることが確認できます。
Rumpの例題を計算してみます。
>>> from mpmath import mp
>>> mp.dps = 50
>>> mp.pretty = True
>>> a = mp.mpf('77617')
>>> b = mp.mpf('33096')
>>> 333.75*b**6+a**2*(11*a**2*b**2-b**6-121*b**4-2)+5.5*b**8+a/(2*b)
-0.82739605994682136814116509547981629199903311578438
真値が得られました。
mp.mpf()は実数を扱う場合に用います。mp.mpc()で複素数、mp.matrix()で行列を引数に渡すことができます。
UNIX bcによる任意精度演算
このようにPythonで精度保証付き数値計算、任意精度浮動小数点演算を行うにはdecimalモジュールまたはmpmathモジュールが必要になります。使い勝手も簡単ではなく、慣れが必要です。計算結果が真値であるかどうかを確かめるために、私はUNIXのbcプログラムを使っています。
bcの任意精度演算は楽で便利です。bcによるRumpの例題
>>> scale = 50
>>> a=77617
>>> b=33096
>>> 333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/(2*b)
>>> -.82739605994682136814116509547981629199903311578439
モジュール不必要、精度を指定するだけで真値が得られます。このようにPythonとbcの両方を使いながら精度を確認しています。
