浮動小数点数演算に関するメモ

float x86

とあるコードで正しく動くであろうマシンで得られた結果と、
開発中のツールから得られた結果を比較して、開発中のものを
テストするというものがあったのですが、浮動小数点の精度の
違いがもしかしたら出るのではないかということで。

x86の浮動小数点に関する問題

浮動小数点演算ではまった話 - bkブログ
80 bit float number problem


等々によく解説されているので詳しくはそちらを見てください。
x86の浮動小数点は特殊です。なんで x86マシンで計算した浮動小数点の結果と
他のアーキテクチャでの浮動小数点数の結果を比較して、挙動が正しいかを
判定するという場合は注意が必要です。

問題の回避方法

SSE2では 64ビット精度での各種演算の命令が使えます。なんでこちらの
方を使えば 80ビットとして扱われるということはなくなります。


上記のサイトのプログラムを SSE2命令を使うようビルドしてみます。

#include <stdio.h>

int main (void)
{
    double a, b, c;

    a = 1.0 / 1e307; // aはすごく小さな数
    printf("a= %25.20e\n", a);

    b= a / 3e16;   // aをある数で除して、同じ数を乗ずる
    b= b * 3e16;   // aに戻るか? 64bit演算だと戻らなくて正常
    printf("b= %25.20e\n", b);

    c= a / 3e16 * 3e16; // 内部表現が80bitだと戻ってしまう
    printf("c= %25.20e\n", c);

    return 0;
}

コンパイル

   % gcc -msse2 -mfpmath=sse double.c
   % ./a.out
   a= 1.00000000000000010695e-307
   b= 1.48219693752373963253e-307
   c= 1.48219693752373963253e-307 // 64ビット精度なので戻らない

   % gcc -mfpmath=387 double.c 
   % ./a.out
   a= 1.00000000000000010695e-307
   b= 1.48219693752373963253e-307
   c= 1.00000000000000010695e-307 // 80ビット精度なので戻る

gccの場合 "-mfpmath"オプションで浮動小数点演算について
命令を切り替えることができます。
sseを指定する場合は -msse2をあわせてつけてください。
80bit精度で行う場合は -mfpmath=387としてください。


gccでは x86_64向けにビルドされたものでは mfpmathオプションは
デフォルトで sse2になっているため、特に意識しなければ意図通り
64ビット精度で計算されます。386向け gccでは 387がデフォルトです。
(Macの場合は 32ビットでもデフォルトで SSE2が有効になっているようです。)


clangだとデフォルトは sse2のようです。ドキュメントを少し見た限り
80bit精度で計算する方法がわかりませんでした。まあ移植等々を考えると
その方針の方がいいかなとも思います。

最後に

浮動小数点演算は気をつけよう。