-ffinite-math-onlyでも動くNaN/infチェック

以前こういう記事を書いた。

in-neuro.hatenablog.com

まあNaNじゃないことを仮定しての高速化なんだからstd::isnan(x)を常にfalseだと思ってdead code eliminationするのは当たり前じゃんという気はするが、備忘録的に。

とはいえ、コンパイルオプションに-Ofastを付けられたとしてもNaNチェックをしたいことはある。どうしたらよいか。コンパイラにNaNチェックをしていると思われなければよい。要するに、浮動小数点数を使ってNaNチェックをするから駄目なのだ。直接ビットパターンを調べれば良い。

64bit浮動小数点数のレイアウトはWikipediaに載っている。とはいえ符号-指数-仮数の順であることを覚えていれば、符号が1bitなのは自明なので、仮数が52bitであることを覚えてさえいればあとは引き算ができればレイアウトがわかる。

ja.wikipedia.org

特殊な浮動小数点数は指数部が全て1の浮動小数点数として定義されている。仮数部が全部0ならinfで、それ以外はNaNになる。なので、指数部のビットが全て1であるかを確認すれば、isfinite相当のものは実装できる。

まず、doublestd::uint64_tに変換する。以下が鉄板のパターンで、これ以外だと大抵規格違反になる。

bool is_NaN_or_inf(const double x)
{
    std::uint64_t n;
    std::memcpy(reinterpret_cast<char*>(&n),
                reinterpret_cast<const char*>(&x),
                sizeof(double));
    // ...
}

ポインタや参照を直接変換すると、アライメントの条件が変わってしまうことがある。そのせいで、たいていの型の間のポインタ・参照変換は規格違反となる。アライメントはメモリアドレスに対する制限で、例えばある型は8の倍数のアドレスにしか置けないというような制限だ。だが、(unsigned|signed) charとstd::byteへの変換はexplicitに許されている。これはアドレスに対する制限の最小単位だからだ。

なのでdoublestd::uint64_tの領域の先頭のポインタをchar配列へのポインタだと思い込んで移し替えれば良い。これは合法にビットレベルキャストをするための頻出パターンなのでコンパイラは(可能な場合)普通にmov命令に変換する。なのでパフォーマンスへの影響などはあまり恐れる必要はない。

さて、このstd::uint64_tのビットパターンを確認する。特定のビットだけが問題なので、まずそこ以外をマスクして関係のないビットを全て0にする。その後、関係のあるビットの1を0に、0を1にひっくり返せば、関係のあるビットが全て1なら結果的に全てのビットが0になり、全体が0になる。逆に0が混じっていれば1のビットが残り、非ゼロになる。その実装は以下のようにできる。ここで、仮数部は見ていないのでinfとnanの両方が引っかかることに注意したい。

bool check_NaN_or_inf(const double x) noexcept
{
    std::uint64_t n;
    std::memcpy(reinterpret_cast<char*>(&n),
                reinterpret_cast<const char*>(&x),
                sizeof(double));
    // bin: 0111'1111'1111'0000'0000'...'0000
    // hex:    7    F    F    0    0 ...    0
    constexpr std::uint64_t mask = 0x7FF0000000000000;
    return ((n & mask) ^ mask) == 0;
}

C++11の範囲で書いたが、C++14から整数リテラルに区切り文字として'を使えるようになったので、上のように0が連続する場合は使ったほうがいい。0が一個多かったりしてもこれではわからないので。

-Ofastにしても消し去られてはいないことをチェックしよう。以下のリンクで普通のstd::isfiniteと見比べるとわかりやすい。std::isfinite版は即値を戻り値レジスタに入れてすぐreturnしている。std::isfiniteの戻り値にnotをつけないと結果が一致しない(今回実装した関数はNaNかinfの時にtrueを返すので)のにつけ忘れたままembedしてしまった。まあブログだしいいや。

godbolt.org

アセンブリをよく見ていると、実際memcpymovqになっていることがわかる。それと、ビット演算部分が少し最適化されて出力されている。このサイトはアセンブリの頻出命令はカーソルホバーすると解説が出てくるので活用してほしい。

NaNチェックの素直な実装であるreturn x != xよりはbit演算分だけ命令数が多くなるが、-Ofastでも使えるNaNチェックができた。あまりホットスポットでないところで使いましょう。

float版は以下の通り。

bool is_NaN_or_inf(const float x) noexcept
{
    std::uint32_t n;
    std::memcpy(reinterpret_cast<char*>(&n),
                reinterpret_cast<const char*>(&x),
                sizeof(float));
    // bin: 0111'1111'1000'0000'0000'...'0000
    // hex:    7    F    8    0    0 ...    0
    constexpr std::uint32_t mask = 0x7F800000;
    return ((n & mask) ^ mask) == 0u;
}