gccの最適化指示-Ofastは危険

はじめに

gccの最適化指示である-Ofastをお気軽に使ってる記事を見掛けたので注意喚起的なやつです。

-Ofastとは何ぞや

gccのドキュメントから引用
https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html#index-Ofast

-Ofast
Disregard strict standards compliance. -Ofast enables all -O3 optimizations. It also enables optimizations that are not valid for all standard-compliant programs. It turns on -ffast-math, -fallow-store-data-races and the Fortran-specific -fstack-arrays, unless -fmax-stack-var-size is specified, and -fno-protect-parens. It turns off -fsemantic-interposition.

要は規格を無視した最適化をやっちゃうよ、ということ。

実害が出る例

#include <stdio.h>

double hoge(const double v[4])
{
    return v[0] * v[1] * v[2] * v[3];
}

double v[4] = {1e-100, 1e-100, 1e200, 1e200};

int main(void)
{
    printf("%g * %g * %g * %g = %g\n", v[0], v[1], v[2], v[3], hoge(v));
}

$10^{-100}\times10^{-100}\times10^{200}\times10^{200}$ を計算しています。

これをx86-64版のgcc 12.2を使用してコンパイル、実行すると

-O3を指定してコンパイルしたのの実行結果

1e-100 * 1e-100 * 1e+200 * 1e+200 = 1e+200

-Ofastを指定してコンパイルしたのの実行結果

1e-100 * 1e-100 * 1e+200 * 1e+200 = inf

-Ofastを指定してコンパイル、実行した場合には計算途中でオーバーフローしたのかinf(無限大)となりました。

どうしてこうなった?

普通に考えると $10^{-100}\times10^{-100}\times10^{200}\times10^{200}$ は先ず $10^{-100}$ に $10^{-100}$ を掛けて $10^{-200}$ になって、次に $10^{200}$ を掛けて $10^{0}=1$ になって、最後に $10^{200}$ を掛けるので答えは $10^{200}$ になる筈です。IEEE754の倍精度は $1.7976931348623157\times10^{308}$ 位が上限値なのでオーバーフローには余裕があります。

計算を行ってる関数hoge()のコンパイラが吐いたコードを見てみましょう。

-O3を指定したコンパイル結果

hoge:
        movsd   (%rdi), %xmm0    # v[0]の値を%xmm0に読み出す
        mulsd   8(%rdi), %xmm0   # v[1]の値を%xmm0に掛ける
        mulsd   16(%rdi), %xmm0  # v[2]の値を%xmm0に掛ける
        mulsd   24(%rdi), %xmm0  # v[3]の値を%xmm0に掛ける
        ret                      # 返り値は%xmm0

Cのソースプログラムの意図通りのコード生成がされています。

次いで-Ofastを指定した場合は

-Ofastを指定したコンパイル結果

hoge:
        movsd   (%rdi), %xmm0    # v[0]の値を%xmm0に読み出す
        movsd   16(%rdi), %xmm1  # v[2]の値を%xmm1に読み出す
        mulsd   8(%rdi), %xmm0   # v[1]の値を%xmm0に掛ける
        mulsd   24(%rdi), %xmm1  # v[3]の値を%xmm1に掛ける
        mulsd   %xmm1, %xmm0     # %xmm1の値を%xmm0に掛ける
        ret                      # 返り値は%xmm0

-O3を指定した場合とコードの内容が異なります。-O3を指定した場合と比べて命令数が増えているため一見遅そうな印象ですが

        movsd   (%rdi), %xmm0    # v[0]の値を%xmm0に読み出す

と

        movsd   16(%rdi), %xmm1  # v[2]の値を%xmm1に読み出す

は依存関係がないためスーパースカラーにより並列実行ができ、

        mulsd   8(%rdi), %xmm0   # v[1]の値を%xmm0に掛ける

と

        mulsd   24(%rdi), %xmm1  # v[3]の値を%xmm1に掛ける

についても同様なので実行速度的にはこちらの方が有利のようです。

問題の箇所は

        movsd   16(%rdi), %xmm1  # v[2]の値を%xmm1に読み出す

        mulsd   24(%rdi), %xmm1  # v[3]の値を%xmm1に掛ける

で、$10^{200}\times10^{200}$ の計算となるためオーバーフローが生じたようです。
浮動小数点演算に於いて、計算の途中でオーバーフローが生じないよう計算の順序をプログラマが考えることは普通にありますが、-Ofastの前にはそういった努力が無駄になる可能性があります。また、計算の順序が変わることで結果への影響が生ずる可能性もあるため、-Ofastの使用は余程の注意が必要でしょう。

ちなみに、「計算の順序を強制する様カッコで括ればいんじゃね?」と思って試してみたのですが

#include <stdio.h>

double hoge(const double v[4])
{
    return ((v[0] * v[1]) * v[2]) * v[3];
}

double v[4] = {1e-100, 1e-100, 1e200, 1e200};

int main(void)
{
    printf("((%g * %g) * %g) * %g = %g\n", v[0], v[1], v[2], v[3], hoge(v));
}

-Ofastを指定してコンパイルしたのの実行結果

((1e-100 * 1e-100) * 1e+200) * 1e+200 = inf

無駄な努力でした。

おわりに

おわりです。