非正規化数と flush to zero フラグについて

非正規化数とは

たいていの計算機で浮動小数点数は IEEE754 に定める単精度（binary32)や倍精度（binary64) 形式で表現されます。これからは半精度（binary16）が流行るかもしれませんが。
いずれにしても浮動小数点数の全ビットパターンは、次の5種類のいずれかになります。

非正規化数は、有限数の絶対値最小の数（指数部 000...01, 仮数部ゼロ、Fortran では TINY 関数で得られる）よりも小さい絶対値を表します。

指数部	仮数部
111...1	000...0	Infinity (無限大)
111...1	ゼロ以外	NaN（非数）
ゼロと最大の間	任意	有限数
000...0	ゼロ以外	非正規化数
000...0	000...0	ゼロ

単精度では TINY(1.0) = 1.17549435E-38, 倍精度では TINY(1d0) = 0.222507385851-307 がその境目ですね。

非正規化数の例

ビットパターン例を挙げたほうがよいのでプログラム動かしましょう。これは倍精度の例です。

subroutine printz(a, label)
  implicit none
  real(8), intent(in):: a
  character(*), intent(in):: label
  integer(1):: buf(8)
  character(8):: judge
  buf = transfer(a, buf)
  buf(1:8) = buf(8:1:-1)
  if (a == 0.0d0) then
    judge = 'ISZERO'
  else
    judge = 'NONZERO'
  endif
  judge = adjustr(judge)
  write(*, '(a16,": ",8z3.2," (",g24.17,")",a8)') label, buf, a, judge
end subroutine

program test
  implicit none
  real(8):: a, b, c, d, e
  real(8):: kakezan
  call printz(4.0_8, '4')
  call printz(2.0_8, '2')
  call printz(1.0_8, '1')
  call printz(0.5_8, '1/2')
  a = tiny(1.0_8)
  b = a * 3.0_8
  d = a * (1.0_8 + 1.0_8 / rrspacing(1.0_8))
  c = b * 0.5_8
  call printz(b, '3 tiny(1d0)')
  call printz(c, '3/2 tiny(1d0)')
  call printz(d, 'tiny(1d0)+')
  call printz(a, 'tiny(1d0)')
  c = c - a
  e = c * (1.0_8 + 2.0_8 / rrspacing(1.0_8))
  call printz(e, '1/2 tiny(1d0)+')
  call printz(c, '1/2 tiny(1d0)')
  call printz(d - a, 'minimum')
end program

実行結果

               4:  40 10 00 00 00 00 00 00 (  4.0000000000000000    ) NONZERO
               2:  40 00 00 00 00 00 00 00 (  2.0000000000000000    ) NONZERO
             1.5:  3F F8 00 00 00 00 00 00 (  1.5000000000000000    ) NONZERO
               1:  3F F0 00 00 00 00 00 00 (  1.0000000000000000    ) NONZERO
             1/2:  3F E0 00 00 00 00 00 00 ( 0.50000000000000000    ) NONZERO
     3 tiny(1d0):  00 28 00 00 00 00 00 00 ( 0.66752215755216041-307) NONZERO
   3/2 tiny(1d0):  00 18 00 00 00 00 00 00 ( 0.33376107877608021-307) NONZERO
      tiny(1d0)+:  00 10 00 00 00 00 00 01 ( 0.22250738585072019-307) NONZERO
       tiny(1d0):  00 10 00 00 00 00 00 00 ( 0.22250738585072014-307) NONZERO
  1/2 tiny(1d0)+:  00 08 00 00 00 00 00 01 ( 0.11125369292536012-307) NONZERO
   1/2 tiny(1d0):  00 08 00 00 00 00 00 00 ( 0.11125369292536007-307) NONZERO
         minimum:  00 00 00 00 00 00 00 01 ( 0.49406564584124654-323) NONZERO

普通の値（4〜1/2）の挙動から、指数部は十六進表現で先頭3文字、ビットでいうと先頭12ビット（そのうち一番最初のビットは符号なので、11ビット）であることが知れます。
2の冪乗（4, 2, 1）の場合は仮数部は全部ゼロで、その中間の値は仮数部の値が大きくなっていくわけです。

tiny(1d0) より小さな数はこの指数部がゼロになっていて、非正規化数 denormalized number となります。

仮数部の最下位ビットを1に下のが「2の冪のすぐ隣の値」ですが、見比べることで、 tiny(1d0) の直上でも、その下のオーダーでも、ゼロ近辺でも、 0.00000000000000005-307 = 5.0e-324 程度の差は表現できることがわかります。

どうして発生するか

普通に四則演算で発生します。

• 加減算：値が tiny(1.0) や tiny(1d0) 未満の違いの値の減算（符号違いで絶対値がそういう差になる加算）
• 乗除算：普通の有限値に対して、絶対値が小さくなるような除算・乗算
• 関数： exp, hypot, sin, tan, atan, atan2, sinh, tanh, asinh, atanh などがアンダーフローを起こしうるとされているので、結果が非正規化数となりえます。

面白いのは cos や log ではアンダーフローが発生しないことですね。これは結果がゼロになる引数がゼロではないので、そのへんでは入力値のすべてについて確認すれば、非正規化数が発生しないことを断言できるわけです。

どうして嫌がられるか

さて、この非正規化数、継子扱いされるのですが、演算の種類によるわけです。

• 加減算のためには、「表現できる値の差」が一様なので、あったほうがよい
• 逆数を取るとオーバーフロー例外または無限大になるので、除算では怖い

演算の結果、非正規化数が発生することを、アンダーフローといいます。アンダーフローの場合に特別な処理をしようとすると、 IEEE754 の範囲では、アンダーフロー例外を発生されるという手があります。ありますが、それって 普通の UNIX では SIGFPE シグナルを受けてプロセスが死ぬ、ということなので、死んでほしいわけではない場合には嬉しくないわけですね。

FTZ (flush to zero) フラグと DAZ (denormal as zero) フラグ

そんなわけで MMX 以後の x86 プロセッサには、FTZ フラグというのができました。演算結果が非正規化数になるときに、ゼロに変えてしまうというものです。MMX のコントロールワードのフラグを立てることにより MMX/SSE/AVX 命令がそういう動きをします。言い換えると、8087 由来の命令には効きません。そもそも x86 の浮動小数点数処理には 8087 系の命令と MMX 系の命令があるわけですが、その古いほうですね。まあ使わないと思うけど。
さらに確か SSE3 以後の x86 には DAZ フラグというのができています。オペランドが非正規化数の場合にゼロ扱いするというものです。

似たような flush to zero 機能が ARM aarch64 にもあるらしいです。NEON てのが SIMD 命令のことですね。SIMD じゃない昔風の浮動小数点命令があったのかなかったのか記憶が定かでありません。

x86 + C 言語

gcc だったら pmmintrin.h というので操作することができます。DAZ だけを有効にすることもできます。

#include <pmmintrin.h>

_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);

intel c でも使えます。

Aarch64 の場合はどうするのかわかりません。

x86 + intel fortran

intel fortran の場合にはいちおう -ftz オプションで FTZ が有効になるのですが、デフォルトの -O2 に含まれるので、どっちかというと FTZ を無効にしたい場合に -O0 または -no-ftz を打ち込む感じです。

アセンブリ吐いて見ると、FTZ 有効時にはメインプログラムの冒頭に何やらコードが入ってフラグを立てています。所詮はプロセッサのフラグですので、その後で上記 pmmintrin.h でほげると上書きすることもできるはずです。

x86 + gfortran

gfortran の場合は -ffast-math で FTZ が有効になります。何か他に悪いことがあるかも知れません。

実行結果（FTZ）

               4:  40 10 00 00 00 00 00 00 (  4.0000000000000000    ) NONZERO
               2:  40 00 00 00 00 00 00 00 (  2.0000000000000000    ) NONZERO
             1.5:  3F F8 00 00 00 00 00 00 (  1.5000000000000000    ) NONZERO
               1:  3F F0 00 00 00 00 00 00 (  1.0000000000000000    ) NONZERO
             1/2:  3F E0 00 00 00 00 00 00 ( 0.50000000000000000    ) NONZERO
     3 tiny(1d0):  00 28 00 00 00 00 00 00 ( 0.66752215755216041-307) NONZERO
   3/2 tiny(1d0):  00 18 00 00 00 00 00 00 ( 0.33376107877608021-307) NONZERO
      tiny(1d0)+:  00 10 00 00 00 00 00 01 ( 0.22250738585072019-307) NONZERO
       tiny(1d0):  00 10 00 00 00 00 00 00 ( 0.22250738585072014-307) NONZERO
  1/2 tiny(1d0)+:  00 00 00 00 00 00 00 00 (  0.0000000000000000    )  ISZERO
   1/2 tiny(1d0):  00 00 00 00 00 00 00 00 (  0.0000000000000000    )  ISZERO
         minimum:  00 00 00 00 00 00 00 00 (  0.0000000000000000    )  ISZERO

aarch64 + TCS fortran

TCS fortran の場合は -Kfz で aarch64 の FTZ が有効になります。