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

なんだか最近最大公約数を高速に求めたい読者（競技プログラマ？）が多いようなので、書いておきます。

最大公約数を求める改良版アルゴリズムとして、Stein's algorithm（Binary GCD）というものがあります。

ナイーブな実装では、普通の互除法よりも遅いですが、ctz（末尾に何ビット0が続くか）を分岐なしに高速に計算できれば実用的になります。

Steins; GCD - Qiitaでは、表引きを用いたctz（本文中ではntz）実装が使われています。 このctz実装は、C++11のconstexpr関数として実装できていることから明らかなように、ポータブルではありますが、乗算とメモリアクセスがあるため低速です。

高速性を重視する場合、ポータビリティは失われますが、__builtin_ctzllが使えます。なお、__builtin_ctzllはx86ではBSF命令にコンパイルされます。

生成されるコードの重要な部分に分岐命令が含まれないよう、注意深くコーディングした以下のコードは、std::gcdの3.6倍ほど高速です。

uint64_t gcd_stein_impl( uint64_t x, uint64_t y ) {
    if( x == y ) { return x; }
    const uint64_t s = x > y ? x - y : y - x;
    const int n = __builtin_ctzll( s );
    const uint64_t t = x > y ? y : x;
    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 );
}