◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇

浮動小数点数と数値計算の罠

発行：エスオーエル株式会社
https://www.sol-j.co.jp/

連載「X線CTで高精度寸法測定！？」
2020年10月14日号　VOL.113

平素は格別のお引き立てを賜り、厚く御礼申し上げます。
X線CTスキャンによる精密測定やアプリケーション開発情報などをテーマに、
無料にてメールマガジンを配信いたしております。

◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇◆◇



技術の仕事をしていると、何度も何度も計算します。
手計算、電卓、エクセル、Python、LISP、C++ など手段もいろいろ使います。

計算が好きかという質問には、
「好きと嫌いの重ね合わせ状態」という回答になります。

誰かがやってくれる計算、特にコンピュータがやってくれる計算を
自分がわざわざやりたいという気持ちにはなりません。

でも、計算しなければ見えないものが見えた時の感動を
直接味わいたいので、いろんな計算を自分でしてみたいとも思います。

（ガウスは、単純計算も好んで自分でやったという話を聞いたことがあります。）

特に仕事でする計算には、大きな責任が伴うので、手抜きはできません。


今回は、コンピュータが出力する数値を鵜呑みにすると、
せっかくの知的な楽しみを逃すことになりますというお話をします。

コンピュータの結果を信じると大変なことになりますよ、とまでは言いません。
もはやコンピュータの計算と無縁に生きることは困難な時代になりました。

せめて、
　・丸め誤差
　・桁落ち
　・打切り誤差
　・情報落ち
というものが存在して、重大事故が起きないように考えてくれている人々がいる
ので、感謝しましょうと言いたいです。


さて、浮動小数点数というものがあり、こいつとは長い付き合いです。
もう頻繁に使います。自分でプログラムを書くときだけでなく、
コンピュータが出す結果を吟味するときも、なんだかいつも気が付くと、
考えたり、調べたり、格闘したり、まとわりついたり、楽しんだり、しています。


　　昔、i386 ユーザーだった頃（中～高校生位）、
　　コプロセッサーの i387 にとても憧れました。
　　i386 は CPU ですが、i387 は 浮動小数点演算の専用回路（FPU）です。

　　今は、CTのボリュームデータ、カメラの画像データ、座標データなど、
　　扱うデータのほとんどに浮動小数点演算が関連しています。
　　今の CPU は、FPU を内蔵してしまっているのです。


それでは、コンピュータの結果を信じると大変なことになりますよ、
という脅しによく使われる例を計算してみます。

まず、

　　f(a,b) = (333.75-a^2)*b^6+a^2*(11*a^2*b^2-121*b^4-2)+5.5*b^8+a/(2*b)

という式を用意します。
ここに、

　　a = 77617.0
　　b = 33096.0

を代入して、計算するのですが、3種類の精度で計算してみます。
結果は、（ひとまず、どの精度でも 20桁まで出力させて、）

　　　　単精度（有効桁数は10進で約 7桁）： 1.1726039648056030273
　　　　倍精度（有効桁数は10進で約15桁）： 1.1726039400531786949
　　拡張倍精度（有効桁数は10進で約19桁）： 1.1726039400531786318

となりました。

　　「となりました」というのが自分にとって重要です。
　　「となるそうです」ではありません。
　　どこかの本の丸写しではなく、自分のコンピュータに計算させました。

ちなみに、エクセルに計算させると、「1.1726039400531800000」でした。


さて、精度を上げていくと、数字の一致する桁が増えるので、
「1.1726039」までは信頼できて、
倍精度の「1.1726039400531786」も使えそうだなと思ってしまいます。

でも、真の値は全く違います。

f(a,b) を手計算で頑張ってみます。
まず、a は素数です。
b は素因数分解できて、b = 33096 = 2^3 ×3 ×7 ×197 です。

桁落ちしないように、上手く順序を並び替えて、計算してみると、

　　f(a,b) = a/(2*b) - 2   = -54767/66192

まで進めることができます。
なんと、コンピュータで計算させると、正の値ですが、
本当は負の値だということが分かります。

続けて、-54767/66192 も計算したいのですが、
コンピュータはもはや信用できないし、手計算も面倒なレベルです。

そこで、昔習っていたソロバンの登場です。パチパチすると、答えは、

　　f(a,b) = -0.827396...

となりました。

いや～、恐ろしいですね。
というお話なのですが、この結果自体よりも、
（このメルマガを書くために）この結果を得るための労力がものすごかった。。


まず、目標として、3つの異なる精度で計算した数値が、
「精度を上げていくと、数字の一致する桁が増える」という結果を得なければ、
この話は成り立ちません。
3つの結果が暴れると、何か計算がおかしいと気づいてしまいます。
結果を比べても間違いに気づかず、逆に誤った自信を持ってしまうという話です。


この話の元ネタである Rump さんの論文オリジナルの式：

　　f = 333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/(2*b)

を使うと、成功しませんでした。使ったコンピュータも言語も違うからです。
本当に、よくこんな例を作ったなと思います。
浮動小数点演算の性質を熟知していないとできないと思います。

どうやら、浮動小数点演算で、いろいろ問題が起こらないように、
IEEE 754 という規格も頑張っているようです。

でも、この Rump さんの例題は、
IEEE 754 ならば回避できるという類ではないようです。

そこで、このメルマガを書くために、
またもや浮動小数点数と格闘する羽目になりました。


どうやら、自分がよく使っている VC++ では、
この結果を得られないということが分かりました。

VC++ では、doulbe型 と long double型 が同じもののようです。
（ sizeof() で確認すると同じ型サイズでした。）

そこで、急遽パソコンに g++ (GNU C++) と VS Code の環境を作りました。
（VS Code ユーザーではありませんでしたが、使ってみるととても良い！）

g++ では、「float」「double」「long double」が精度の異なる
浮動小数点数型になっています。
（ sizeof() で確認すると、それぞれ 4, 8, 16 で型サイズが違います。）

次に、Rump さんオリジナルの式を
意味が変わらないように変形して、欲しい結果が得られるようにしました。
ちょっと変形すると全く違った値になります。そして、それは真の値でもない。
とても不安定な計算だと分かります。

式を分解して各項を出力してみたり、浮動小数点演算の性質を考えたりと、
結構苦労して、ようやく完成しました。


結論としては、
　　測定には、誤差の評価が重要です。
　　同様にコンピュータの数値計算でも誤差の評価が重要です。
ということなのでしょうが、

いろいろあって、上手い一言で締めくくることができません。
感想：コンピュータに狙った特定の間違いをさせることは大変でした。


--
高野智暢