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;
}