最大公約数をもっと高速に求める番外編（その2）

前回のコードを最適化していきます。

主要なループ*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倍ほど早く動作します。