前回のコードを最適化していきます。
主要なループ*1のうち、最も時間がかかるのはctz部分（3サイクルレイテンシのBSF命令にコンパイルされます）です。
これを少しでも早く実行開始できれば、さらに高速に動作するはずです。
ここで、二の補数表現の特徴からctz(s) == ctz(-s)である点に着目します。
すると、ctz(x < y ? y - x : x - y)などとせずにctz(x - y)としても得られる結果は同じです。
こうすることで、1サイクルでも早くctzの計算に取り掛かれるようになります。
そのように式変形をすると、gccを使った場合は余計なおせっかい*2が発動してかえって遅くなります。 以下は、そのようなおせっかいをしないclang向けのコードです。
uint64_t gcd_stein_impl( uint64_t x, uint64_t y ) {
if( x == y ) { return x; }
const uint64_t a = y - x;
const uint64_t b = x - y;
const int n = __builtin_ctzll( b );
const uint64_t s = x < y ? a : b;
const uint64_t t = x < y ? x : y;
return gcd_stein_impl( s >> n, t );
}
uint64_t gcd_stein( uint64_t x, uint64_t y ) {
const int n = __builtin_ctzll( x );
const int m = __builtin_ctzll( y );
return gcd_stein_impl( x >> n, y >> m ) << ( n < m ? n : m );
}
このコードは、clang++-9 -O3でコンパイルした場合、std::gcdと比べて5.5倍ほど早く動作します。
