先週の記事の続きです。

最大公約数をもっと高速に求める（その3） - よーる

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

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 );
}

分岐を外していい場合

（その1）で示したコードでは、互減法と互除法のどちらか適切な方を選ぶ分岐をするより、10回の固定回ループの方が速いことを使って高速化しました。 これは、分岐予測ミスのペナルティ時間で10回くらいのループは実行できるということによっています。 ちなみに、10回というのをちゃんと最適化すると、8回か9回がベストなようです。

今回は、80回の固定回ループとなっているため、それと同じ論理では固定回数ループの方が速いとは言えません。 つまり、互減法と互除法のどちらか適切な方を選ぶ分岐をすることで高速化する可能性が残されています。

nがmの何倍大きいときに互除法が適切だと判断するべきかという問題は、実験してみると32倍あたりに最適解があるようです。 これはループ6回分に相当します。

（その1）で示したコードでは、ループ一周当たり3cycleで8周がベストという結果でした。つまり、分岐ミスペナルティと引き換えに無駄な仕事をする価値があるのは24cycleくらいだということになります。 今回のコードでは、ループ一周当たり4cycleで6周がベストという結果だったので、やはり24cycleくらいだということが支持されます*1

そこで、以下のようなコードとしてみます。

uint64_t gcd_impl( uint64_t n, uint64_t m ) {
    constexpr uint64_t K = 5
    if( m == 0 ) { return n; }
    for( ; n / 32 < m; ) {
        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 );
}

このようにすることで、さらに一割ほど高速化しました。

---

今週はちょっと忙しくて記事を書く時間が取れなかったのでちょっと少ないですが、ここまで！

ここまでで、std::gcd関数の2.4倍ほど速いgcd関数を作れています。