g++のOfastオプションとNaNの話

C++のプログラムがNaNの周りで不可解な挙動をしたので追ってみたという記事です

結論から言うと、-Ofast（-ffast-math）は有限の浮動小数点数しか扱わないので気を付けようということです*1

何もわかってない部分も多いのでツッコミ歓迎です :bow:

背景と概要

最近レイトレの高速化を少しやっていて、その過程でsimdを使ってみたがその結果をうまく取得することができなかった

取得したい結果は_mm256_cmp_pdの大小比較（_CMP_LE_OS）の結果で、仕様によるとtrueのときは0xFFFFFFFFFFFFFFFFが入り、falseの時は0x0が入る

gccではresults[0]のようにして結果の__mm256d型の各要素をdouble型として取得できるので、results[0] != 0やstd::signbit(results[0])のようにして結果の取得ができると思った*2が、_mm256_cmp_pdの結果と先ほどの判定が一致せず困った

NOTE: 僕の環境はIntelなので今後の内容は主にIntel® 64 and IA-32 Architectures Software Developer Manuals を参照してます

浮動小数点数

今回の話では0xFFFFFFFFFFFFFFFF（全bitが1）と0x0（全bitが0）を内部に持つdouble型が問題になる

一般に、浮動小数点数の指数部のbitが全て1だと無限大やNaNなど特殊な扱いになることが多く、この辺りはアーキテクチャにもよるらしい（wiki）

マニュアルを読むと、僕の環境では0x0は+Zeroで、0xFFFFFFFFFFFFFFFFはNaN*3に属するとわかる

通常、これらを区別するには0との比較（bitがどこかに立ってればfalse）やstd::signbit（最上位bitが立ってればtrue*4）を使えば十分なはずである

// compile option: -std=c++20 -O3

#include <iostream>
#include <bit>
#include <cstdint>
#include <cmath>

int main() {
    // 最適化を防ぐために入力を噛ませている
    // 実行時に 0 を与えて mask64 == UINT64_MAX にしている
    int shift; std::cin >> shift;
    uint64_t mask32 = (1LL<<32) - (1<<shift);
    uint64_t mask64 = (mask32 << 32) + mask32;

    double d = std::bit_cast<double>(mask64);
    std::cout << d << std::endl;               // output: nan
    std::cout << (d==0) << std::endl;          // output: 0
    std::cout << std::signbit(d) << std::endl; // output: 1
}

コンパイラオプション

今回のプログラムは高速化を目的としていたこともあり、高速化用のコンパイラオプションとして-Ofast -march=nativeをつけていた

この-Ofastは-O3にコンパイラオプションを更に追加したもので、その中でも今回は-ffast-mathが関係している

-ffast-math自体も複数のコンパイラオプションから構成されている*5が、具体的には--ffast-mathの中でも(d==0)のほうは-ffinite-math-onlyオプションで結果が変わり、std::signbit(d)のほうは-fno-signed-zerosオプションで結果が変わるようになる

よって-Ofast環境下では、NaNやInfinityを含むdouble型を操作をする場合にはstd::bit_castなど用いて*6bit表現を介することで処理するのが（おそらく）正攻法になる

つまりは今回のようにNaNに対して何か操作をしようとした時点で未定義動作なので何も保証はないのだが、ただの0との比較やstd::signbit判定において何を最適化する部分があるのか気になったのでもう少し深追いしてみた

// compile option: -std=c++20 -O3 -ffast-math

#include <iostream>
#include <bit>
#include <cstdint>
#include <cmath>

int main() {
    // 最適化を防ぐために入力を噛ませている
    // 実行時に 0 を与えて mask64 == UINT64_MAX にしている
    int shift; std::cin >> shift;
    uint64_t mask32 = (1LL<<32) - (1<<shift);
    uint64_t mask64 = (mask32 << 32) + mask32;

    double d = std::bit_cast<double>(mask64);
    std::cout << d << std::endl;               // output: nan
    std::cout << (d==0) << std::endl;          // output: 1
    std::cout << std::signbit(d) << std::endl; // output: 0
}

アセンブラの比較（d==0のケース）

#include <iostream>
#include <bit>
#include <cstdint>
#include <cmath>

int main() {
    // 最適化を防ぐために入力を噛ませている
    // 実行時に 0 を与えて mask64 == UINT64_MAX にしている
    int shift; std::cin >> shift;
    uint64_t mask32 = (1LL<<32) - (1LL<<shift);
    uint64_t mask64 = (mask32 << 32) + mask32;

    double d = std::bit_cast<double>(mask64);
    std::cout << d << std::endl;               // output: nan

    // output: 1 if compiled with -ffinite-math-only else 0
    std::cout << (d==0) << std::endl;
}

上記のコードを-O3オプションでコンパイルしたものと-O3 -ffinite-math-onlyでコンパイルしたものを比較してdiffをとった

intel記法で、xmm0とedxはゼロクリアされ、xmm1にdが入っており、edxの値が直後に出力される状況である

# -O3
<       mov     eax, 0
<       ucomisd xmm1, xmm0
<       setnp   dl
<       cmovne  edx, eax
# -O3 -ffinite-math-only
>       comisd  xmm1, xmm0
>       sete    dl

-O3 -ffinite-math-only（後者）について

comisd  xmm1, xmm0 : xmm1とxmm0の比較を行う
sete    dl : ZFが立っている場合edxの最下位byteを1にする（つまり今回だとedxを1にする）

かなり直感的に処理をしていることがわかる

ただし、comisdの仕様を見ると以下のように記述されており、UNORDEREDの場合でもゼロフラグが立つため、xmm1にnanが与えられた場合には常に1が返ってきて、非直感的な結果となるのが確認できる

Operation
COMISD (All Versions)
RESULT := OrderedCompare(DEST[63:0] <> SRC[63:0]) {
(* Set EFLAGS *) CASE (RESULT) OF
  UNORDERED: ZF,PF,CF := 111;
  GREATER_THAN: ZF,PF,CF := 000;
  LESS_THAN: ZF,PF,CF := 001;
  EQUAL: ZF,PF,CF := 100;
ESAC;
OF, AF, SF := 0; }

-O3（前者）について

mov     eax, 0 : eaxに0を代入
ucomisd xmm1, xmm0 : xmm0とxmm1を比較
setnp   dl : PFが0のときedxの最下位byteを1にする
cmovne  edx, eax: ZFが0のときedxにeaxを代入する（つまりedxを0にする）

先ほどと比べてもう少し何かをしてそうなことがわかる

COMISDがGREATER_THANまたはLESS_THANのときに限ってZF=0となるので、最後のcmovne命令より、xmm1 > xmm0 or xmm1 < xmm0のときはedx=0である

残りのケースであるCOMISDがUNORDEREDのときとEQUALの区別としてsetnpが使われており、EQUALの場合のみPF=0となるのでxmm0=xmm1の時に限ってedx=1となるコードになっていることが確認できる

比較まとめ

xmm1がnanのケースを丁寧に弾くか弾かないかで結果と命令数が変わるというだけのことだった

ちなみに、comisdとucomisdの違いはsource operandがQNaNの時に例外シグナルを飛ばすかどうかで、comisdのみシグナルを飛ばす（SNaNの時は両方シグナルを飛ばす）

-ffinite-math-only側でcomisdが使われてるのはNaNなんか知るかという意志を感じる（NaNをサポートしてないので姿勢としては正しそう？）

アセンブラの比較（std::signbit(d)のケース）

#include <iostream>
#include <bit>
#include <cstdint>
#include <cmath>

int main() {
    // 最適化を防ぐために入力を噛ませている
    // 実行時に 0 を与えて mask64 == UINT64_MAX にしている
    int shift; std::cin >> shift;
    uint64_t mask32 = (1LL<<32) - (1LL<<shift);
    uint64_t mask64 = (mask32 << 32) + mask32;

    double d = std::bit_cast<double>(mask64);
    std::cout << d << std::endl;               // output: nan

    // output: 0 if compiled with -fno-signed-zeros else 1
    std::cout << std::signbit(d) << std::endl;
}

上記のコードを-O3オプションでコンパイルしたものと-O3 -fno-signed-zerosでコンパイルしたものを比較してdiffをとった

intel記法で、xmm1にdが入っており、edxの値が直後に出力される状況である

# -O3
>       mov  rcx, rbx
<       movmskpd        edx, xmm1
<       and     edx, 1
# -O3 -fno-signed-zeros
>       xor     edx, edx
>       mov  rcx, rbx
>       pxor    xmm0, xmm0
>       comisd  xmm0, xmm1
>       seta    dl

-O3（前者）について

mov  rcx, rbx: diffの対応の意味で入れたが、関係ないので無視する
movmskpd        edx, xmm1: xmm1の最上位bitが立っていればedxの最下位bitを1とする
and     edx, 1: edx&1を計算してedxに格納する

すごく単純で何も問題なさそうだが、このmovmskpdがどうやらsimd演算命令で、（多分）命令として遅いのが問題っぽい

仕様読んだところ、この演算はxmm1がの128bit目も見に行くとのことで、128bit目の情報が下から2番目のbitに格納されるので最後にand`をしている

-O3 -fno-signed-zeros（後者）について

xor     edx, edx: edxをゼロクリアしている
mov rcx, rbx: 関係ない処理。データハザード回避のために迷い込んでいると勝手に予想している
pxor    xmm0, xmm0: xmm0をゼロクリアしている
comisd  xmm0, xmm1: xmm0とxmm1の比較
seta    dl: CF=0かつZF=0のときedxの最下位byteを1にする

comisdでCF=0かつZF=0というのはxmm0 > xmm1の時のみ成立するので、つまりはxmm1=nanだとedx=0となる

それ自体はそれでいいのだが、なぜ-fno-signed-zerosという+Zeroと-Zeroを区別しないというオプションが関係あるのかと考えると、確かにxmm1=-0.0のとき、本来は最上位bitが立っているのでedx=1となってほしいが、xmm0=0.0=-0.0=xmm1よりedx=0となることがわかる

普通のケースでは0未満かどうかを見れば最上位bitが立っているかどうかはわかるので確かに良さそうで、命令数は増えたものの重い命令はないから速い（ということなんだろう）

比較まとめ

-fno-signed-zerosオプションの有無によってnanの扱いも変わるなど、思わぬところにも挙動に違いが出ることがわかった（未定義なので何されても仕方ない）

また、当たり前だが単純に命令数だけではなく命令の重さもちゃんと考慮して最適化してそうなことも伺えた

あと、素人目にはsimdでしか最上位bit取得する命令がないの不思議な感じがする

まとめ

-ffast-mathは結構怖いので使うときは気を付けたほうがよさそう（今回触れたもの以外でも怖そうな最適化をしている）

以下の記事でも、パフォーマンスのチェックに-ffast-mathを使うが、実際にはそのままオプションを足すのではなく、差分と安全性を考慮したうえでソースコードに手を加えるようにしていると書かれている*7 kristerw.github.io

*1:X上ではO3オプションといいましたが、Ofastでした

*2:なるべくcast系の命令を挟みたくなかった

*3:細かくは QNaN だがこの記事ではあまり関係ない

*4:https://en.cppreference.com/w/cpp/numeric/math/signbit にあるように、NaNに対しても使える

*5:https://kristerw.github.io/2021/10/19/fast-math/ がわかりやすい

*6:union型使ってましたが、https://x.com/t5nat/status/1713359396372377625?s=20 で bit_cast の存在を教えてもらいました、ありがとうございます

*7:僕は面倒なので気を付けて Ofast 使うかという気持ちになっているが...