小数の計算精度の奇妙な現象

すごくどうでもいい話ではあります。 この文章を読んでも、多分何も得るものはありませんし、 実際にこんな現象が役に立つ場面がそうそうあるとは思えません。 それでもいいやという奇特で暇な方だけ、このまま読み進めてください。 また、以下の文章はfloatなどの浮動小数点数の表現についてある程度の知識があることを 前提としています。 別に大した知識は要らないと思うので、まあ気軽に読んでください。

※以下のプログラムは、Pentium4上のWinXP, VC7.1のReleaseビルド(←重要)でコンパイルしています。

前提として

ご存知の通り、C/C++言語には、浮動小数点数型と呼ばれる型が3種類用意されています。 float、double、long doubleです。 何で3種類あるかというと、精度がそれぞれ違うからです。当たり前です。 Windows上でのごくごく一般的な環境では、floatはIEEE754単精度浮動小数点数形式(以下単精度)、 doubleはIEEE754倍精度浮動小数点数形式(以下倍精度)として実装されているはずです。 どうしてこの形式が使われているかというと、 単純明快、Intel系のCPUはIEEE754形式の浮動小数点表現を使用しているからです。 CPUに合わせてデータの表現形式を決めるのは、それはそれは当然のことです。

例その1

さて、この表現形式を使用すると、何かにつけて誤差というものが発生します。 誤差にも色々な種類がありますが、今回のポイントとなるのはこんな誤差です。

001 
002 #include <stdio.h>
003 
004 int main() {
005   volatile float f = 1.0f;
006   f += 0.00000001f;
007   printf("f = %.20f", f);
008 }

floatは(doubleもですが)飛び飛びの値しか取れません。 floatの表現形式では、1.0の次の値は、1.0000001192092896ですから、0.00000001を足したところで、 四捨五入されて1.0になってしまいます。

f = 1.00000000000000000000

さて、先のプログラムには奇妙な点が一つあります。 既にお気づきのことと思いますが、volatileです。 volatileは「その変数を計算するたびにメモリから読み込む」 という指示をコンパイラに与える修飾子で、主にマルチスレッドなプログラムで使われます。 当然、上のプログラムでは何の意味も無い……ように思えるので、 volatileを外してコンパイル、実行してみます。

f = 1.00000000999999990000

どうでしょう。0.00000001(に近い値)が足されているではありませんか。 この値は、float型では表現不可能な値です。 つまり、6,7行目において、fはfloat型ではなくなっているのです。

例その2

もう一つ、例を挙げてみましょう。 float型を計算する時の誤差の例としてよく使われる、こんなプログラムです。

001 
002 #include <stdio.h>
003 
004 int main() {
005   volatile float f = 0.0f;
006   for(int i = 0; i < 1000000; ++i) {
007     f += 0.000001f;
008   }
009   printf("f = %.20f", f);
010 }

よくある「誤差が蓄積していって云々」なプログラムです。 期待を裏切らず、こんな出力を吐いてくれます。

f = 1.00903892517089840000

さて、ここでvolatileを外すとどうなるか、 大体予想はついていると思いますが、こうなります。

f = 0.99999999747524271000

精度が7桁ほど違っています。 先ほどの例と同じく、volatileを外した方が精度が良くなるのです。

原因

どうしてこんなことが起こるのか? その答えをきちんと理解するためには、アセンブリ言語まで辿らなくてはいけません。 そこまでやるのは面倒なので、少し大雑把な説明で済ませます。

答えは、IntelのCPUが浮動小数点数を計算するときの仕組みに潜んでいます。 Intelが発行している、 『IA-32 インテル® アーキテクチャ・ソフトウェア・デベロッパーズ・マニュアル、上巻:基本アーキテクチャ』 から引用してみましょう。

x87 FPUデータ・レジスタは、8 つの80 ビット・レジスタで構成される。 値は、これらのレジスタに拡張倍精度浮動小数点フォーマットで格納される。 実数値、整数値、またはパック形式BCD整数値を、 図4-3. に示すフォーマット(筆者注:各種整数・小数フォーマット)のいずれかでメモリから x87 FPUデータ・レジスタのいずれかにロードすると、 値が拡張倍精度浮動小数点フォーマットになっていない場合は値が自動的にこのフォーマットに変換される。 計算の結果をx87 FPU レジスタからメモリに戻すときは、 拡張倍精度浮動小数点フォーマットのままにしておくことも、 倍精度または単精度浮動小数点フォーマット、整数フォーマット、 またはパックドBCD整数フォーマットに変換することもできる。

少し分かり辛いというか訳分からないと思いますが、おもいっきり簡単に要約すると、 「メモリ上での表現形式が何であれ、小数の計算をする際には全部、『拡張倍精度浮動小数点数形式(以下拡張倍精度)』 で計算されるよ」ということです。 逆に言えば、計算結果をメモリに記録してしまうと、メモリ上での精度の低い表現形式に束縛されて、 数値の精度が落ちてしまうのです。 だから、volatileをつけると精度が落ちてしまったわけです。 また、メンバ変数やグローバル変数についても、 同様の現象(メモリ読み書きによる精度落ち)が発生する可能性があります (参考:とってもごはん)。 出来るだけローカル変数を利用するようにしましょう。

ちなみにこの拡張倍精度、80ビットの長さで、10進数にしておよそ18桁の精度を誇ります。 単精度の3倍近い精度です。 long doubleはこの拡張倍精度に対応させればいいと思うのですが、 何故かVCだとlong doubleはdoubleと同じ倍精度なのです。何でだろ?