先週の記事の続きです。

前回示したコードは、以下のようなものでした。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
        for( int i = 0; i < 10; ++i ) {
                uint64_t t = n - m;
                bool q = t < m;
                n = q ? m : t;
                m = q ? t : m;
                if( m == 0 ) { return n; }
        }
        return gcd_impl( m, n % m );
}

uint64_t gcd( uint64_t n, uint64_t m ) {
        return n > m ? gcd_impl( n, m ) : gcd_impl( m, n );
}

このコードの機械語命令列

gcd_impl関数のループ一回分の機械語命令列は、以下のようになっています（実際にはループアンローリングされており、レジスタ割り付けだけが異なる命令列が3回周期で現れます）。

sub    rdi, rsi    # t = n - m;
cmp    rdi, rsi    # q = t < m;
mov    rax, rdi
cmovb  rax, rsi    # n = q ? m : t;
cmovnb rdi, rsi    # m = q ? t : m;
test   rdi, rdi    # m == 0
jz     Label       # if( m == 0 )

ループ一回分は7命令になっているわけですが、だからといって7 cycleかかるわけではありません。近年のCPUでは、これくらいの命令列は並列に実行されます。

今回の命令列で、並列実行のボトルネックになるのは依存の連鎖です。具体的には、t = n - mを計算し、その結果を使ってq = t < mを計算し、その結果によって次ループのnとmが算出される部分です。

qはtが算出されなければ計算できませんし、nとmはqが算出されなければ計算できません。また、次のループのtはnとmが算出されなければ計算できません。

つまり、どんなに並列に実行できるCPUを持ってきても、この部分は各々1cycleかかってしまいます。

よってこのコードでは、ループ一回分に3cycleかかることになります。

なお、mov命令が気になるかもしれませんが、最近のCPUでは実質的に0cycleで実行されると考えてよいようです。

互減法をさらに高速化する

互減法の一ステップでは、以下のように計算が進みます。

GCD(n, m) = if n-m >= m then GCD(n-m, m) else GCD(m, n-m)

これを拡張し、以下のようにしてみます。

GCD(n, m) = if n-Km >= m then GCD(n-Km, m) else if n-m >= m then GCD(n-m, m) else GCD(m, n-m)

もしn-Km >= mであれば、互減法のKステップを一気に進めることができるため、互減法の効率が上がることが予想されます。また、それにより、アルゴリズム全体に占める除算の割合が減少し、高速化すると予想されます。ただし、以下の点に注意する必要があります。

途中結果n-Kmが負になる可能性がある。
Kmがオーバーフローする可能性がある。

1.については、移項することで回避可能です。 2.については、前処理でmを小さくしておけば回避可能です。

改良されたメインループ

メインループの内部を以下のようにしてみます。メインループの前に前処理を置くことでうまくmは小さくしているとして、m*Kはオーバーフローしないと仮定します。

uint64_t t = n - m;
uint64_t s = n - m * K;
bool q = t < m;
bool p = t < m * K;
n = q ? m : t;
m = q ? t : m;
if( m == 0 ) { return n; }
n = p ? n : s;

n = p ? n : sが追加されたため、さらに依存の連鎖が伸び、このループ一回には最低でも4cycleかかるようになりました。もし、m*Kの計算に4cycle以上かかる場合、そのせいで依存の連鎖が伸び、ループ一回にかかるサイクル数はもっと長くなります。あるいは、m*Kの計算が他の命令の実行の裏で並列に実行できない場合も、ループ一回にかかるサイクル数はもっと長くなります。

最近のIntel社製CPUであれば、MUL命令は3cycleで完了するため、m*Kの計算のせいで依存の連鎖が伸びるということはありません（Kの値によってはコンパイラ最適化によりLEA命令を駆使した命令列となり、もっと速く計算されます）。また、m*Kの計算は、他の命令の実行の裏で行うことができるため、そのせいでループ一回にかかるサイクル数が長くなるということはありません。

元のコードのループ一回には3cycleかかっていたところ、4cycleかけることで、Kステップ一気に実行できる可能性が追加されたことになります。 nとmの大きさに大きな差がなかった場合でも、少なくとも1ステップは進むため、それほど損したことにはなりません。

`K`として何をとるか

Kとして2や3などの小さな数をとれば、n-Km >= mである可能性がそれなりに高いため、改良により互減法のステップを効率よく進めることができるようになります。 Kとして大きな数をとる場合、n-Km >= mである可能性が低く、意味のある仕事をできる確率（pがfalseになり、n = sが実行される確率）が下がります。しかし、意味のある仕事をできた場合、互減法のステップを一気に進めることができます。

このトレードオフを考慮に入れ、いくつかの実験を行った結果、以下にようにすることでstd::gcdより二倍以上高速なgcd関数を作れることがわかりました。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
        constexpr uint64_t K = 5;
        for( int i = 0; i < 80; ++i ) {
                uint64_t t = n - m;
                uint64_t s = n - m * K;
                bool q = t < m;
                bool p = t < m * K;
                n = q ? m : t;
                m = q ? t : m;
                if( m == 0 ) { return n; }
                n = p ? n : s;
        }
        return gcd_impl( m, n % m );
}

uint64_t gcd_pre( uint64_t n, uint64_t m ) {
    for( int i = 0; i < 4; ++i ) {
        uint64_t t = n - m;
        bool q = t < m;
        n = q ? m : t;
        m = q ? t : m;
        if( m == 0 ) { return n; }
    }
    return gcd_impl( n, m );
}

uint64_t gcd( uint64_t n, uint64_t m ) {
        return n > m ? gcd_pre( n, m ) : gcd_pre( m, n );
}

gcd_pre関数は、m*Kがオーバーフローしないようにするための前処理です。互減法を4回繰り返すことで、mは大きくとも0xffff'ffff'ffff'ffff / 5であると保証できます。

Kとして5をとったことにより、32のような大きな商が立つ場合でも、6回程度ループを回れば互除法の一ステップを行ったことになります。これ以上大きな商が立つ場合は、除算を行った方がよさそうです。そのようなことの発生頻度が低いため、単純な互減法であれば10回程度で除算に切り替えたほうがよかったところ、80回まで増やすことが可能になりました。なお、gcd_implのループ回数80回は目安であり、必ずしも最適な値とは限りません。

このコードの機械語命令列

互減法をさらに高速化する

改良されたメインループ

Kとして何をとるか

`K`として何をとるか