2007年2月28日
浮動小数点演算ではまった話
浮動小数点演算のありがちな問題ではまりました。 いろいろ調べているうちに x86 特有のちょっとおもしろい 現象に遭遇したので紹介したいと思います。
パーセンテージの計算
簡単な C のプログラムでパーセンテージを計算しようと思い、 次のようなコードを書きました。int x, y; ... int a = (double)x / y * 100;
int a = x * 100 / y としないのは、 x が大きい場合に x * 100 が オーバーフローを起こす (INT_MAX を越える) ためです。
このコードは一見、期待通りに動いていたのですが、 しばらく使っていると、手元の環境では x = 53, y = 100 のときに a は 53 ではなく 52 になることに気づきました。 これは次の理由によります。
式の最初の (double)53 / 100 の計算結果 0.53 は double (手元の環境では 64ビット) では正確に表現できないため、 0.529999... と近似値になります。その結果、 53.0 / 100 * 100 は 52.9999... となり、 最終的に int に変換する際に少数点以下を切り捨てると 52 が得られるというわけです。
この問題は double で正確に表現できない数 0.53 を発生させたのが 原因ですから、次のように(double)x * 100 / y のように先に100をかけてやれば 問題ありません。 double の仮数部は 52ビットあるので、 32ビットの int の値に 100をかけても誤差は生じません。
int x, y; ... int a = (double)x * 100 / y;
変数にいったん入れると結果が変わる問題
というわけで、元の問題は解決したのですが、ついでに いろいろ実験をしていると、おもしろいことに気づきました。
次のコードをコンパイルして実行すると、手元の環境 (32ビットの x86 + gcc 4.1.0) では 52 53 と表示されます (volatile はコンパイラの最適化を防ぐために入れています)。
#include <stdio.h>
int main () {
volatile double x = 53.0;
volatile int a = x / 100.0 * 100.0;
volatile double y = x / 100.0 * 100.0;
volatile int b = y;
printf("%d %d\n", a, b);
return 0;
}
a と b の計算はやっていることはどちらも基本的に同じですから、 結果に違いが出るのは意外です。計算結果を一度 y に代入するだけでなぜこのような違いが出るのでしょうか。
gcc -S でコンパイルしてアセンブリコードを見ると、 a と b を 計算するコードにほとんど違いはなく、唯一の違いは fstpl, fldl 命令を使っているストア・ロードを行っている部分だけです (太字部分)。
#APP
/* a のコード */
#NO_APP
fldl -16(%ebp) /* メモリの内容 (53.0) をレジスタにロード */
fldl .LC1 /* メモリの内容 (100.0) をレジスタにロード */
fdivrp %st, %st(1) /* 53.0 / 100.0 を計算 */
fldl .LC1 /* メモリの内容 (100.0) をレジスタにロード */
fmulp %st, %st(1) /* 0.5299999... * 100.0 を計算 */
fnstcw -42(%ebp) /* コントロールワードをメモリにストア */
movzwl -42(%ebp), %eax /* それを EAX にロード (上位16ビットは0) */
orw $3072, %ax /* AX の中身に 0xc00 のフラグをたてる */
movw %ax, -44(%ebp) /* AX をメモリにストア */
fldcw -44(%ebp) /* コントロールワードを設定 */
fistpl -48(%ebp) /* 52.9999... を切り捨ててメモリにストア */
fldcw -42(%ebp) /* コントロールワードを元に戻す */
movl -48(%ebp), %eax /* メモリの内容 (52) を EAXにロード */
movl %eax, -20(%ebp) /* EAX の内容をメモリ (変数a) にストア */
#APP
/* b のコード */
#NO_APP
fldl -16(%ebp)
fldl .LC1
fdivrp %st, %st(1)
fldl .LC1
fmulp %st, %st(1)
fstpl -32(%ebp) /* 計算結果をメモリにストア (y に代入) */
fldl -32(%ebp) /* ふたたび浮動小数点レジスタにロード */
fnstcw -42(%ebp)
movzwl -42(%ebp), %eax
orw $3072, %ax
movw %ax, -44(%ebp)
fldcw -44(%ebp)
fistpl -48(%ebp)
fldcw -42(%ebp)
movl -48(%ebp), %eax
movl %eax, -36(%ebp) /* EAX の内容をメモリ (変数b) にストア */
#APP
これで原因がわかってきました。 Binary Hacks の Hack #98 の 「x86が持つ浮動小数点演算命令の特殊性」にあるように、 x86 では 浮動小数点の演算は 80 ビットのレジスタで行われます。 80ビットの計算結果を double サイズ (手元の環境では64ビット) のメモリに戻すときには 64 ビットに丸められます。 この仕様により、次のような違いが生まれます。
- a の場合
- 53.0 / 100.0 * 100.0 を 80ビットで計算して 52.9999... が発生し、それを切り捨てた 52 が a に代入される。
- b の場合
- 53.0 / 100.0 * 100.0 を 80ビットで計算して 52.9999... が発生し、y に代入される。 このとき 64ビットへの丸めが行われ、y には丸められた値 53.0 が代入される。 53.0 を切り捨てした 53 が b に代入される。
このように、メモリにストアされるときに行われる 80ビットから 64ビットへの暗黙的な変換のおかげで、 b の方はたまたま 53 という正確な結果が得られています。 実際、上のアセンブリコード内の太字の部分のコードをコメントア ウトすると、実行結果はどちらも 52 になります。
浮動小数点演算の途中の結果を変数に代入すると結果が変わる というのはなかなか不可解な挙動です。 レジスタとメモリで共通のサイズ・形式で浮動小数点数を格納するプロセッサ ではこのようなことは起きないはずです。
What Every Computer Scientist Should Know About Floating-Point Arithmetic の D.11.1 にこの「変数にいったん入れると結果が変わる問題」に関する解説が載っています。
80ビットの浮動小数点数の中身を見る
浮動小数点レジスタの中身を 64 ビットに丸めずに、 80ビットのままメモリにストアするには、 fstpl ではなくfstpt 命令を使います (首藤さんに教えてもらいました)。
それでは、 53.0 / 100.0 * 100.0 を 80ビットで計算した結果を見てみましょう。
#include <stdio.h>
void print_float80(char *p) {
int num_bits = 0;
for (int i = 9; i >= 0; --i) {
const char c = p[i];
for (int j = 7; j >= 0; --j) {
printf("%d", (c & (1 << j)) != 0);
++num_bits;
if (num_bits == 1 || num_bits == 16) {
printf(" ");
}
}
}
printf("\n");
}
int main(int argc, char **argv) {
char result[10];
double x = 100;
double y = 53;
double z = 100;
// y / x * z == 53 / 100 * 100
__asm__("fldl %1\n\t"
"fldl %2\n\t"
"fldl %3\n\t"
"fdivp\n\t"
"fmulp\n\t"
"fstpt %0\n\t"
: "=m" (result)
: "m" (z), "m" (x), "m"(y));
print_float80(result);
return 0;
}
上のプログラムをコンパイルして実行すると、次の結果が表示されます。
0 100000000000100 1101001111111111111111111111111111111111111111111111111111111111
最初の 0 は符合ビット、次の15ビットは指数部、最後の64ビットは仮数部です。 この場合、符合は正、指数部は 5 (2の5乗 = 32)、仮数部は 1.656249999999999999891579782751449556599254719913005828857421875 となります。仮数部の値に指数部の値をかけると 52.9999999999999999965305530480463858111761510372161865234375 という値が得られます。 確かに 52.9999... という値になっていることがわかります。 仮数部の精度は64ビットなので信頼できるのは最初の19桁程度です。
IEEE 754 の 32ビットと 64ビットの浮動小数点数には仮数部の前に 暗黙の 1 (暗黙ビット) が置かれることになっていますが、80ビット の浮動小数点数には 歴史的事情により 暗黙ビットは存在しません。 浮動小数点数の表現については コンピュータの構成と設計 や Write Great Code に詳しく解説されています。
ちなみに、2進数を扱う計算には irb を使うと便利です。
% irb >> 0b100000000000100 => 16388 >> 0b100000000000100-(2**14-1) => 5 >> '%08b' % 5 => "00000101"
まとめ
パーセンテージの計算で生じた切り捨ての問題、 変数にいったん入れると結果が変わる問題、 および 80ビットの浮動小数点数の中身を見る方法、 を紹介しました。
追記
光成さんによるgcc の long double はなぜ12バイトなのかという記事もどうぞ。
追記その2
x87 の丸めの方向は fldcw 命令で変更することができます。glibc なら _FPU_SETCW を使うと便利です。上記のコードを以下のように変更すると、53 53 という結果に変わります。
int main () {
// この2行を追加
fpu_control_t control_word = _FPU_RC_UP | _FPU_IEEE;
_FPU_SETCW(control_word);
volatile double x = 53.0;
volatile int a = x / 100.0 * 100.0;
volatile double y = x / 100.0 * 100.0;
volatile int b = y;
printf("%d %d\n", a, b);
return 0;
}