最大公約数を求める非常に有名なアルゴリズムとして、ユークリッドの互除法が知られています。
このアルゴリズムは、二つの入力を十進法で表したとき、小さいほうの桁数の高々五倍の回数の剰余演算で最大公約数が求まることが知られています（ラメの定理）*1。二つの入力が符号なし64bit整数であれば、高々92回です。これは最悪ケースであり、多くの場合はかなり少ない回数で済みます。実験してみると、平均的には37回程度のようです。
しかし、剰余演算は非常に遅い演算です。Intelの最近のCPUでは、21cycle～83cycleかかるとされています。ここから計算すると、二つの符号なし64bit整数の最大公約数を求めるには、1μ秒弱の時間がかかることになります。
実際、以下のコードを実行してみると、乱数生成を除いて7秒ほどかかります。つまり一回当たり0.42μ秒ということで、だいたい計算と合います。
#include <iostream>
#include <random>
#include <chrono>
// 16M
std::pair<uint64_t, uint64_t> arr[0x100'0000];
int main() {
std::mt19937 mt;
std::uniform_int_distribution<uint64_t> dist( 1, 0xffff'ffff'ffff'ffff );
for( auto& [e1, e2] : arr ) {
e1 = dist( mt );
e2 = dist( mt );
}
const auto start = std::chrono::system_clock::now();
for( auto& [e1, e2] : arr ) {
std::gcd( e1, e2 );
}
const auto finish = std::chrono::system_clock::now();
std::cout << std::chrono::duration_cast<std::chrono::nanoseconds>(finish - start).count() / 1e9 << " sec" << std::endl;
}
コンパイル：clang++ -std=c++17 -march=native -O3
環境：Windows subsystem for Linux、clang5.0.1、Intel(R) Core(TM) i5-7200U@2.50GHz（ターボブースト時3.1GHz）
最大公約数を求めるだけの問題に出会うことはそれほどありませんが、アルゴリズムの内部で最大公約数を求める必要があり、ユークリッドの互除法を部品的に使うアルゴリズムはそれなりにあります。そういったアルゴリズム群の（定数倍）高速化のために、gcd関数を高速化することを考えます。
互除法の仕組み
まずは、ユークリッドの互除法で最大公約数（greatest common divisor: GCD）が求まる理由を復習します。以下では、GCD(n, m)と書いたらnとmの（数学的事実としての）最大公約数を意味することとし、gcd(n, m)と書いた場合は、（帰納的な定義を持つ）最大公約数を求める関数の定義か、それにnとmを代入したものを意味することとします。
まず、GCD(n, m) = GCD(n-m, m)であることが知られています（これの証明は大学入試問題として有名ですが、ここでは割愛します*2）。そのため、gcd(n, m) = if n < m then gcd(m-n, n) else gcd(n-m, m); gcd(0, m) = mといった形で最大公約数を求める関数を、帰納的に定義することができます。これは互減法というアルゴリズムであり、これを使うとコードサイズを小さくできる可能性があります（下記参照。golfer向けであり、実際にプログラムサイズを小さくするのに役に立つとは限りません）。しかし、最悪の入力、つまりgcd(99999,1)などに対して99999回の引き算が発生するため、低速です。
http://golf.shinh.org/reveal.rb?Modular+Inverse/tails_1484205185&c
2019/03/04 01:50 追記：（注意）上記コードは、互減法そのものではありません。モジュラ逆数を求めよという問題に対する解答であり、いわゆる拡張ユークリッドの互除法を使うのが一般的であるところ、互除法ではなく互減法を用いてコードを縮めた例です。
このコードを提出したtailsさんから、混乱を招く恐れがあるとの指摘をいただいたため、追記しました。
GCD(n, m) = GCD(n-m, m)を繰り返し適用すると、GCD(n, m) = GCD(n-m, m) = GCD(n-2m, m) = ... = GCD(mod(n, m), m)となることがわかります。nとmは対称なので、n≧mと仮定しても一般性を失いません。この条件のもと、gcd関数は、gcd(n, m) = gcd(m, mod(n, m)); gcd(n, 0) = nのように定義することができます。これがユークリッドの互除法です。
互除法の各ステップの解析
いくつか手で実験をしてみればわかりますが、[1, N]の区間から一様ランダムに取り出した二つの整数n, mに対する互除法の場合、商は1が立ち剰余はn-mになる、というパターンが多いことがわかります。
これの理由は単純で、商として2以上が立つには、片方がもう一方の二倍以上である必要があるからです。そのようになる確率は、二数の選び方が一様ランダムであれば1/2です。
互除法を進めていくと確率は片寄っていきますが、商として1が立つパターンが頻発することは変わりません。また、商が1でなくとも、商は小さいことが多く、大きな商の出現は稀です。一つ一つの数に対して確率が小さいのはもちろんのこと、それらすべての確率の総和も小さいです。
N=264-1として実験してみると、1が立つのは41%、2が立つのは17%、3が立つのは9%、4が立つのは6%、5が立つのは4%、6が立つのは3%、までで八割を占めています。そのあとを追ってみると、13以下が90%、26以下が95%、128以下が99%です。一様ランダムな場合、商はその逆数に比例した確率で出現します。互除法を進めていくにつれ、一様ランダム性が崩れる（区間が減少するのに加え、偶然小さい数が出現する確率が大きくなる）ため解析は複雑になります。しかし、小さい数が出現すると互除法が速やかに終了するため、全体的な確率には大きな影響を及ぼさないようにも思われます。結局のところ、商の出現確率はその逆数に比例していると考えても、大きく外れてはいないでしょう。
剰余演算は遅いので避ける
剰余演算を使うことで単純な引き算と比べて高速化できる、というのが互除法でした。しかし、剰余演算は非常に遅い演算であり、できる限り避けたいところです。
互除法の各ステップにおいて、商として1が立つ場合、数十サイクルもかけて除算を行うのは非効率的であり、1cycleで終わる引き算を行えばよいことになります。
先ほどの解析によれば、商として1が立つことは多い割合で発生するため、数十サイクルのロスを頻繁に行っていることがうかがえます。
商として何が立つかは実際に割り算を行わないといけないため、そのような効率化はできないように思えるかもしれません。
しかし、m != 0の時、`n >= m && m > n - mがn / m == 1と等価なため、以下のようにgcd関数を作ることができます。
uint64_t gcd_impl( uint64_t n, uint64_t m ) {
if( m == 0 ) { return n; }
uint64_t t = n - m;
if( m > t ) {
return gcd_impl( m, t );
} else {
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 );
}
しかし、これを使っても速度が数%改善される程度にとどまります。これは、商として1が出るパターンの出現頻度は四割であり、かつそれは予測不可能に出現するため、分岐予測が全く当たらないのが原因です。
当てられない分岐は遅いので避ける
高速化のテクニックに、当てられない分岐を取り除くというものがあります。
m == 0も予測不可能に出現しますが、これなしというわけにはいかないこと、出現は最後に一回だけであり影響は軽微であることから、この分岐を取り除くことは考えません。
n > mは、入力の最初に50%の確率で発生し、これは絶対に当てられません。
よって、m > tの分岐を取り除くことを考えます。その前にちょっと準備をします。
予測できない分岐なしに互減法を実装する
予測できない分岐を取り除くため、x86の命令セットには、cmov (conditional move)命令が用意されています。これは分岐なしにcond ? a : b;のような計算をすることができる命令です。昔のIntelはさぼっていたのかこの命令のレイテンシが2cycleでしたが、最近では1cycleとなっているため、さらに気軽に使えるようになりました。cmov命令を使う欠点は、分岐すれば計算する必要がなかったかもしれないのに、aもbも計算しないといけなくなる点です。
互減法は以下のように書けます。
while( m != 0 ) {
uint64_t t = n - m;
bool q = m > t;
n = q ? m : t;
m = q ? t : m;
}
このコードを書いたとしても、必ずしもcmov命令が使われるとは限りません。実際g++でコンパイルすると条件分岐を使ったコードにコンパイルされてしまいます。一方clang++でコンパイルするとcmov命令が使われました。
互除法と互減法を組み合わせる
先ほどの実験では、商として1が現れることが多く、特に128以下の商が立つ状況は99%の確率で発生することを明らかにしました。そこで、互除法ではなく互減法の方が速いのではないか、と思いたくなりますが、これは誤りです。一様ランダムに選ばれた符号なし64bit整数同士の互除法は実験してみると平均37ステップかかることがわかります。つまり、三回に一回は128を超える商がどこかのステップで出現することとなり、そのような場合に引き算を続けるのは非常に非効率的です。
一方、多くの場合に商が1であり、そのような場合は除算より引き算の方が一桁速いことを考えると、互減法が優れている局面は確かに存在することもわかります。
よって、うまく互除法と互減法を切り替える、または組み合わせるなどの手法が有効そうであることがわかります。
ここで、切り替えるために条件分岐を使ってしまうのは、やはり当てられない分岐となるためによくなさそうです。実際、先ほどの失敗は当てられない分岐を作ってしまったことによるものでした。
分岐せずにうまく切り替える方法はあるのでしょうか？
キーとなるアイデアは、「分岐予測に失敗して時間を浪費するくらいなら、その時間に無駄になるかもしれない仕事をしておけば、その仕事が無駄ではなかった場合に高速化する」というものです。このアイデアはニュートン・ラフソン法を用いている平方根関数を高速化するときにも使いました。これは、「収束判定という予測できない分岐をするかわりに、どんな入力に対しても収束する回数だけ回す」というものです。すでに収束している場合は無駄な加速を行っていることになりますが、それでも分岐予測に失敗するよりは速いのです。
このアイデアを使うと、以下のようなコードになります。
uint64_t gcd_impl( uint64_t n, uint64_t m ) {
for( int i = 0; i < 10; ++i ) {
uint64_t t = n - m;
bool q = m > t;
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 );
}
ここで、ループの回数10回は適当であり、最適化の余地があります。このコードは、互減法を10回行った後、互除法的な除算を一回やる、というのを繰り返すものになっています。商が小さければ互減法の部分でnとmが小さくなってくれるため、トータルの除算回数が減ることになり、高速化します。一方、商が大きくなる場合でも、10回しか引き算せずに除算を行うようになっているため、延々と引き算をし続ける、といった非効率的な状況にはまることを避けられます。
なお、10回のループに関してはループアンローリングが行われることを期待しています。もしそうでなくても、少なくとも当てられる分岐にはなっています。よって、この関数の中には、if( m == 0 ) { return n; }という避けようのない分岐を除いて、予測不可能な分岐はありません。
先ほどのアイデアがどう使われているかを別の言葉でも説明してみます。
もし商が大きい場合、10回の引き算しても意味はなく、なるべく早く除算モードに切り替えることが高速化につながります。しかし、商が大きいという予測不可能な分岐をしてモードを切り替えるくらいなら、10回程度の引き算という無駄な仕事をした方が高速です。
逆に（運よく）商が小さい場合、引き算は無駄な仕事とならず、高速化に貢献します。
まとめ
最大公約数を求める有名なアルゴリズム、ユークリッドの互除法には除算が含まれており、一般的なCPUで計算すると遅くなりがちです。互減法を部分的に組み合わせて高速化を図ろうとしましたが、予測不可能な分岐が追加されたため、それほど高速化はしませんでした。そこで、無駄かもしれない仕事と引き換えに予測不可能な分岐を除去する方法を考えました。その手法を使うと、元のナイーブな互除法に比べて1.5倍ほど高速に動作するプログラムができました。
