Fast inverse square rootなるものをご存知だろうか。これは、1.0 / sqrt(x)を高速に近似計算する命令である。
Fast inverse square root - Wikipedia
この計算は非常によく使う。例えばベクトルがあるとする。このベクトルを規格化したい。つまり長さを1にしたい。どうすればよいか。 当然、各要素をベクトルの長さで割れば良い。これ以上簡単なことはない。
inline std::array<float, 3> normalize(std::array<float, 3> const& v) { const auto l = length(v); return {{v[0] / l, v[1] / l, v[2] / l}}; }
だが、一般に除算は他の演算に比べて遅い。依然、「アセンブリ解読」においてコンパイラが必死で割り算を消し去っていたことは記憶に新しい。 なので、（即値割り算をあれだけいじくり回すなら、多分コンパイラはオプションによっては上のコードでも勝手にやるだろうが）割り算をなくそう。
inline std::array<float, 3> normalize(std::array<float, 3> const& v) { const auto inv_l = 1.0 / length(v); return {{v[0] * inv_l, v[1] * inv_l, v[2] * inv_l}}; }
ところで、length()の実装はどうなっているだろうか。
inline float length(std::array<float, 3> const& v) { return std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); }
まあこんな感じだろう。ということはだ。
const auto inv_l = 1.0 / std::sqrt(...);
ここに、fast inverse square rootが出てきている。 この手の演算は、CGやゲーム・科学技術計算で死ぬほど出てくる。精度が要求される場合は微妙だが、近似計算で間に合う場合はこの計算はできるだけ早く行いたい。
という背景を思えば、以下の謎の関数が存在する理由もわかってもらえるのではないだろうか。
float Q_rsqrt( float number ) { long i; float x2, y; const float threehalfs = 1.5F; x2 = number * 0.5F; y = number; i = * ( long * ) &y; // evil floating point bit level hacking i = 0x5f3759df - ( i >> 1 ); // what the fuck? y = * ( float * ) &i; y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration // y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed return y; }
Wikipediaからの孫引きで、それによると引用元は"quake3-1.32b/code/game/q_math.c". Quake III Arena. id Software. Retrieved 2017-01-21."だ。 途中で謎のbit演算とマジックナンバーが出現しており、WTFどころではない。
とまあこんな手法まで登場するほどこの処理は本当に使用頻度が高くそれにしてはコストも高いので、intel製のCPUには専用命令を搭載している（rsqrtss）。
この命令の中身は恐らく上記のような形になっているのだろう。これが使える場合は普通にこれを使ったほうがよいので、マクロによる分岐が無限に生じていく。可読性は消え去る。
で、思ったことがある。素人が適当に最適化するよりも、無言でコンパイルオプションを-Ofastにしたら十分速くなるとよく言われる。
さらにClangはLLVMの恩恵を十全に活かして、まるでコードの意味を理解しているかのようなアグレッシブな最適化をするということで有名だ。
Clangがこの専用命令を知っている場合、1/sqrt(x)はその命令を使う形に勝手に最適化できないだろうか。
使うにしてもこれは近似計算なので、オプションで何も言わなければ使わないことも考えられる。どこまで教えれば使うのだろう。
というわけで試してみる。手元のマシンはhaswell世代のCPUで、Ubuntu 16.04、clang-4.0だ。以下のようなコードを書く。
#include <math.h> float rsqrt(float x) { return 1.0 / sqrtf(x); }
とりあえず上記の最適化がなされるか確かめるため、いつものやつでいく。
$ clang-4.0 -std=c99 -ffast-math -march=native -mtune=native -Ofast rsqrt.c -c -o rsqrt_clang.o $ objdump -d rsqrt_clang.o > rsqrt_clang.asm
というわけで出たのが以下のコードだ。
0000000000000000 <rsqrt>: 0: vrsqrtss %xmm0, %xmm0,%xmm1 4: vmulss %xmm1, %xmm0,%xmm0 8: vfmadd213ss 0x0(%rip),%xmm1,%xmm0 11: vmulss 0x0(%rip),%xmm1,%xmm1 19: vmulss %xmm0, %xmm1,%xmm0 1d: retq
少なくとも、専用命令が使われている。では-march=nativeと-mtune=nativeを外してみよう。
0000000000000000 <rsqrt>: 0: rsqrtss %xmm0, %xmm1 4: mulss %xmm1, %xmm0, 8: mulss %xmm1, %xmm0, c: addss 0x0(%rip),%xmm0 14: mulss 0x0(%rip),%xmm1, 1c: mulss %xmm1, %xmm0 20: retq
どうやら効いていたのはFMAだけだったっぽい。-ffast-mathを外してみる。この辺りでこの命令は使われなくなるのでは……と思ったが、先のアセンブリコードと完全に一致した。
で、-Ofastを-O3にしたところ、ついにrsqrtssが消えた。
0000000000000000 <rsqrt>: 0: movaps %xmm0,%xmm1 3: xorps %xmm0,%xmm0 6: sqrtss %xmm1,%xmm0 a: ucomiss %xmm0,%xmm0 d: jnp 1c <rsqrt+0x1c> f: push %rax 10: movaps %xmm1,%xmm0 13: callq 18 <rsqrt+0x18> 18: add $0x8,%rsp 1c: movss 0x0(%rip),%xmm1 23: 24: divss %xmm0,%xmm1 28: movaps %xmm1,%xmm0 2b: retq
もしかしてOfastだと-ffast-mathが勝手に追加されてるのかな（うろ覚え）と思い、-Ofast -ffast-mathでやってみると、復活した。man clangをチラ見したもののどのフラグが自動で立つかは見当たらなかった。
やはり-ffast-mathによる最適化を行えば、rsqrtssが使われるらしい。恐らく理由は、これが近似計算だからだろうと思う。
コンパイラも結構やってくれるんだな、ということがわかった。ちなみにgcc-7.2.0も-ffast-math -O3で同等のバイナリを吐いた。
ここまでずっと32bit浮動小数点数でやってきたが、64bitでやると普通にsqrtを取って割るコードになる。恐らく対応する命令がないのだろう、とは思うが、どうなのだろう。
