非正規化数 - EAGLE 雑記

浮動小数点数には非正規化数(denormal number, subnormal number)というものが存在することを知った．
http://ja.wikipedia.org/wiki/%E9%9D%9E%E6%AD%A3%E8%A6%8F%E5%8C%96%E6%95%B0
これの存在は以下のようなコードで確認できる．

#include <stdio.h>
#include <float.h>

double f(void)
{
  double d = 1.0;
  double e;
  while (d > 0) {
    e = d;
    d /= 2.0;
  }
  return e;
}

int main(void)
{
  printf("%e\n", DBL_MIN);
  printf("%e\n", f());
  return 0;
}

これを普通に gcc でコンパイルして実行したところ，自分の環境では

2.225074e-308
4.940656e-324

という出力が得られた．

一般に非正規化数の演算は遅い．
このような非正規化数を発生させない手段として，SSE においては DAZ (Denormals Are Zeros) や FTZ (Flush To Zero) というビットフラグが MXCSR レジスタに存在する．
インラインアセンブラで STMXCSR/LDMXCSR 命令を直接使ってフラグをセットすることももちろんできるけど，ここではイントリンシック命令を使うことにする．
FTZ は xmmintrin.h で定義されている _MM_SET_FLUSH_ZERO_MODE() で制御できる．
これで FTZ の有効・無効でどれくらい差があるか簡単に調べてみると，

#include <stdio.h>
#include <float.h>
#include <xmmintrin.h>
#include <sys/time.h>

double t(void)
{
  struct timeval tv;
  gettimeofday(&tv, NULL);
  return tv.tv_sec + tv.tv_usec*1e-6;
}

double f(void)
{
  double d = 1.0;
  double e;
  while (d > 0) {
    e = d;
    d /= 2.0;
  }
  return e;
}

void run(double x)
{
  double d = 0;
  double s = t();
  int i;
  for (i = 0; i < 100000000; i++) {
    d += x;
  }
  double e = t();
  printf("%e, %g sec\n", d, e-s);
}

int main(void)
{
  _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
  run(f());

  _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_OFF);
  run(f());
  return 0;
}

ON のときは約0.11秒で計算し結果が2.225074e-300だったのに対し，OFF のときは計算に約5.8秒かかり結果は4.940656e-316だった．

また，GCC においては -funsafe-math-optimizations 付きでコンパイルすると DAZ と FTZ を有効化するようなコードを生成するようだ．
実際に最初のコードを gcc -funsafe-math-optimizations でコンパイルしてやってみたところ，

...
00000000004006e0 <set_fast_math>:
  4006e0:	0f ae 5c 24 fc       	stmxcsr -0x4(%rsp)
  4006e5:	81 4c 24 fc 40 80 00 	orl    $0x8040,-0x4(%rsp)
  4006ec:	00
  4006ed:	0f ae 54 24 fc       	ldmxcsr -0x4(%rsp)
  4006f2:	c3                   	retq
  4006f3:	90                   	nop
...

というような命令列が追加されていて，どこからかここに飛んでセットしているようだった．0x8040は6ビット目(DAZ)と15ビット目(FTZ)だけが立っているビット列．
これを実行してみるとたしかに f() が DBL_MIN と等しい値を返していた．

DAZ は FTZ と同様に pmmintrin.h 内の _MM_SET_DENORMALS_ZERO_MODE() で制御できる．ただし SSE3 以上が必要．

FTZ と DAZ の違いがわからない…
FTZ のみを有効にした場合と DAZ のみを有効にした場合で計算結果は変わらないし計算時間もほとんど変わらなかった*1．

参考

http://software.intel.com/en-us/articles/x87-and-sse-floating-point-assists-in-ia-32-flush-to-zero-ftz-and-denormals-are-zero-daz/

*1:微妙に DAZ のみのほうが遅かったけど