コラム 人と星とともにある数学 数学

Pythonで数値計算 mpmathライブラリ編|Pythonで数学を学ぼう! 第34回

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の両方を使いながら精度を確認しています。

  • この記事を書いた人
  • 最新記事

桜井進(さくらいすすむ)様

1968年山形県生まれ。 サイエンスナビゲーター®。株式会社sakurAi Science Factory 代表取締役CEO。 (略歴) 東京工業大学理学部数学科卒、同大学大学院院社会理工学研究科博士課程中退。 東京理科大学大学院非常勤講師。 理数教育研究所Rimse「算数・数学の自由研究」中央審査委員。 高校数学教科書「数学活用」(啓林館)著者。 公益財団法人 中央教育研究所 理事。 国土地理院研究評価委員会委員。 2000年にサイエンスナビゲーターを名乗り、数学の驚きと感動を伝える講演活動をスタート。東京工業大学世界文明センターフェローを経て現在に至る。 子どもから大人までを対象とした講演会は年間70回以上。 全国で反響を呼び、テレビ・新聞・雑誌など様々なメディアに出演。 著書に『感動する!数学』『わくわく数の世界の大冒険』『面白くて眠れなくなる数学』など50冊以上。 サイエンスナビゲーターは株式会社sakurAi Science Factoryの登録商標です。

あわせて読みたい

-コラム, 人と星とともにある数学, 数学
-