機械イプシロン

少し前にTwitterで、ε-δ論法に関連して、計算機イプシロンなるものがあるよ、というツイートを見かけた。 そして何の気なしにWikipedia で計算機イプシロン&note{machine-epsilon-ja-wikipedia:計算機イプシロン - Wikipedia, 2014-03-18版, 2017-09-18閲覧};について引いたら、MSDN&note{MSDN-machine-epsilon:データ型定数 (CRT), VS2005対応, 2017-09-18閲覧};&note{MSDN-machine-epsilon-vs2015:Data Type Constants, VS2015対応, 2017-09-18閲覧};の説明が誤っていると書いてあった。

まじかよと思って調べてみた。

検証環境 Visual Studio 2013 Professional Update 5 OS Windows 7 64bit Ultimate

結論†

MSDNの定義の通りに解釈すると、計算機イプシロンの意味をなさない定義になってしまう。 MSDNの定義の通りの値だと、FLT_EPSILONより小さい値となってしまう

C/C++におけるfloat†

• 計算機イプシロンの解説をする前に、floatをはじめとする浮動小数点数のフォーマットについておさらいしておく。
• C/C++におけるfloat は32bit幅であり、その内部はbitごとに符号化されている。具体的には
  • 31(MSB) 符号ビット
  • 30:23 指数部
  • 22:0 仮数部
• という具合に、MSB(最上位ビット)から順番に1:8:23 の割合でそれぞれ役割が決まっている。
  • 最上位ビットの符号ビットは、いわゆるsigned int と同じ挙動なので戸惑うことはないと思う。
  • その次の指数部と仮数部は、浮動小数点の浮動小数点らしさを決めている。
  • 指数部は、その数の絶対的な大きさを表す。
  • 仮数部は、数の細かい部分を表す。
  • 合わせて、（指数部）ｘ（仮数部）の掛け算で浮動小数点は表される。
• 仮数部は23ビットあり、0から始まって数が大きくなるほど仮数部も大きくなり、最大0x7fffffまで増えると、その次は指数部が1つ大きくなり、仮数部はまた0に戻る。
• なので、浮動小数点と言いつつも、基本的には離散的にしか実数を表せないことに着目しよう。

計算機イプシロンの定義†

• さて、Wikipediaによれば、計算機イプシロンの定義は、

  「1より大きい最小の数」と1との差のことである

• と定義されている。&note{machine-epsilon-ja-wikipedia};
• 一方で、MSDNによると計算機イプシロンの定義は

  (1+x != 1)となる最小のx

• とある。&note{MSDN-machine-epsilon};
• 一見すると、同じように見えるのだが、実は厳密には違うという点を記しておきたい。

16進数での表記†

• 先に説明したように、floatも16進数で表記すれば離散値で表すことができる。
• このとき、1.0は、16進数表記で 0x3f800000 となる。念のため、先程の表記に合わせると符号ビットが0（正の数）、指数部が0x7f(2^0)、仮数部が0x0(1.0)となる。
• これらの掛け算&mimetex(+1 \times 2^0 \times 1.0 = 1.0);より、1.0が表される
• 0x3f800000 の次に大きい数は、0x3f800001 であり、10進数では結果は1.0000001192092896となる。
• 1との差分は&mimetex(1.192092896^{-7});(0.0000001192092896)である。実際にCでFLT_EPSILONを出力してみると1.1920928955e-007である事がわかる。
• 試しに、どれぐらい小さい値なのか、ちょっと検証コードを書いてみよう。

#geshi(c++){{ return 0; }}

• なおコード全体はこちら。

• コードは、読めばわかると思うが、1 + x が1と等しくなるまで、xをどんどん1/2ずつ小さくしながら繰り返す。
• 結果、xが0.00000011920928960、つまりFLT_EPSILONのときは、加算結果が1ではなくなっているが、xがその半分、0.000000059604644775のときには1+x = 1となっている。 それぞれ、約1/8388608=1/(2^23)と、1/16777216=1/(2^24) である。

• このコードだけ見ると、うっかり「MSDNの説明も合ってるんじゃないの？」と思ってしまうが、実は厳密な意味では誤りである。なぜなら、1+x != 1 となる最小のxは、FLT_EPSILONより小さい数があるからである。

• 今度はもうひとつ検証コードを書いてみよう。
• FLT_EPSILONは正の実数であるので、1.0と0の間に存在する。この間の数を2分探索して、1+x != 1 となる最小の数を調べるコードである。以下の通り。

#geshi(c++){{ return 0; }}

• この実行結果が以下の通り。

1.00000011920928960e+000        5.9604651881e-008       0x33800001

• この通り、FLT_EPSILONの約半分、0.000000059604651881 を加算すると1ではなくなる。 これがなんで起きるのかは、すでに冒頭のWikipediaにも書かれているが、丸め込みである。
• つまり、加算した結果がfloatで表せるちょうどの値にないときは丸め込みが発生し、FLT_EPSILON の次に小さい値、0.000000059604651881を加算すると、結果は1.00000011920928960eに丸め込められる。

• このあたりの丸めの挙動も、本当はround to even、round to zero など、丸め込みの方法があるのだが、今はそのあたりは割愛する。（多分IEEE754に書かれてる気がするんだが、真面目に調べてない）

ULPについて†

• FLT_EPSILONはなるべく実数に近い話を想定して書かれているのでこんがらがるのだが、要は16進数で考えると話が早い。
• 以下に表を記すが、

• という具合に、16進数で表した場合に、いくつずれているか、という問題に帰結できる。
• この16進数で表した際にいくつずれているか、というのを表す単位として、ULP(Unit in last place)がある。
• このULPがいわゆるfloatが表せる細かさの限度なのだが、この細かさは指数部により、絶対値が変わる。そのため、1を基準にして、ULPの大きさが決められている。

1 + x != 1 と 1 - x  != 1

• ULPは、指数部によって、絶対値が変わると述べたばかりだが、1 - x != 1 となる最小の数xは、FLT_EPSILONの半分である。これは、1-xの実際の計算結果は1より小さいため、1よりさらに粒度が細かくなり、結果FLT_EPSILONの半分になる。

なにゆえ、1なのか？†

• もともとのツイートは、ε−δ論法から発展して、FLT_EPSILON に話が発展している。とすると、FLT_EPSILONより、小さい数が存在するなら、FLT_EPSILON の意味は一体何なのか？という話になる。
• 現に、本当にfloatで表せる最小の数として、マクロFLT_MINが定義されており、こちらは 1.17549435082228750e-038 と、文字通り桁違いに小さい値がある。
• じゃあ、なんでFLT_EPSILONは1基準なのかという話なのだが、多分FLT_MINでは使い勝手が悪いからだろうと思う。
• 元来、FLT_EPSILONは十分小さい、という意味で使われている数字である。つまり、例えば計算結果が十分小さいかどうか、という判定に使われる。
• これがFLT_MINだと、本当に最小の数字なので、それより小さい誤差なんてのはありえない。なので、十分に近いかどうかなんて判定できない。というか判定が意味をなさない。
• なので、指数のレンジのちょうど真ん中である、1を基準にFLT_EPSILONを決めたのではないか。（ココらへんは筆者の予想、感想。裏付けなし）

1.1920928955e-007       0x34000000
1.00000000000000000e+000        5.9604644775e-008       0x33800000

最新の20件