はじめに
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を使用してコンパイル、実行すると
1e-100 * 1e-100 * 1e+200 * 1e+200 = 1e+200
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()のコンパイラが吐いたコードを見てみましょう。
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を指定した場合は
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));
}
((1e-100 * 1e-100) * 1e+200) * 1e+200 = inf
無駄な努力でした。
おわりに
おわりです。