色々調べているうちに、「拡張ユークリッドの互除法」というアルゴリズムについて書かれている記事にたどり着きました。
- https://qiita.com/drken/items/b97ff231e43bce50199a
- https://drken1215.hatenablog.com/entry/2018/06/08/210000
これを使えばいいんだ！と思いましたが、どうにも理解が難しい部分がありました（素人なので）。それを小一時間悩んだ末に解明したのでここにメモしておきます。
拡張ユークリッドの互除法とは
先ほど書いたページにもちゃんと書いてあるのですが、今一度確認。
正の整数 $a,b$ について、それらの最大公約数を $d = \gcd(a,b)$ として
$$
ax + by = d \quad\cdots(1)
$$
となる $(x,y)$ を求めるアルゴリズムです。
$$
a = qb + r
$$
とすると、代入して
$$
(qb+r)x + by = d \Leftrightarrow b(qx+y)+rx=d
$$
となります。これで、$a,b$ についての問題を、$b, r$についての問題に置き換えることができます。ここで $qx+y=s, x=t$ とした方程式
$$
bs + rt = d \quad\cdots(2)
$$
の解 $(s,t)$ が求まったならば、元の問題の解は逆算して
$$
x = t, y=s-qx
$$
であるとわかります。
$(s,t)$が求まるようになるのか、という話ですが、式(1)から式(2)で扱う整数が$(a,b)$から$(b,r)$に変わるのは普通の(?)ユークリッドの互除法と同じで、行き着く先は$(d,0)$になります。このときの方程式は
$$
dx+0y=d \quad\cdots(3)
$$
で、これの解は $(x,y) = (1,0)$ です。（ここで$y$は、数学的には何でもいいけど、0にしとくのが無難）
プログラムにする
ここまでは理解できました。さて、これをコードにするにはどうすればいいのか。
先述のサイト https://qiita.com/drken/items/b97ff231e43bce50199a に載っているコードがこちらです
// 返り値: a と b の最大公約数
// ax + by = gcd(a, b) を満たす (x, y) が格納される
long long extGCD(long long a, long long b, long long &x, long long &y) {
if (b == 0) {
x = 1;
y = 0;
return a;
}
long long d = extGCD(b, a%b, y, x);
y -= a/b * x;
return d;
}
何をしてくれるプログラムなのかはわかるし、if(b==0)のところもわかります。先述の式(3)の処理です。
僕がつまづいたのはその後の
long long d = extGCD(b, a%b, y, x);
y -= a/b * x;
return d;
です。これを解き明かしていきます。と、その前にこの関数の使い方を理解しておいた方が良いと思います。
メイン関数を以下のようにしてaとbを入力してあげればx, y, dを教えてくれます。
int main(){
long long a, b;
cin>>a>>b;
long long x, y;
long long d = extGCD(a, b, x, y);
cout<<"x = "<<x<<endl;
cout<<"y = "<<y<<endl;
cout<<"d = "<<d<<endl;
}
//入力
4 13
//出力
x = -3
y = 1
d = 1
当然 $4\times (-3) + 13 \times 1 = 1$ です。
ポイントは、extGCDをやる前に、解を入れる変数（箱）を用意してやるというところですね。参照渡しについての説明は省きますが、先に箱を用意しておいて引数として渡せばextGCDがそこに入れてくれると思っておきます。
解読
まず、long long dですが。これは再帰の底（式(3)）で得られた最大公約数 $d$ を再帰の上の関数に送っているだけです。
問題は再帰しているところ。これをわかりやすくするために、ちょっと冗長な記述に書き換えてみます。if(b==0)のところは最初は省きます。
先ほどの数学的考え方に立ち戻ると、アルゴリズムの最初では $a=qb+r$ となる $q,r$ が必要です。よって
long long extGCD(long long a, long long b, long long &x, long long &y) {
long long q = a/b;
long long r = a%b;
}
とします。
そして、再帰によって $bs+rt=d$ となる $(s,t)$ を求めたいです。ということで $s,t$ を入れる箱を用意してあげて、再帰的に関数を実行します。
long long extGCD(long long a, long long b, long long &x, long long &y) {
long long q = a/b;
long long r = a%b;
long long s, t;
long long d = extGCD(b, r, s, t);
}
extGCDに(b, r, s, t)を入れれば、s, tに $bs+rt=d$ の解が入り、返り値は $b,r$ の最大公約数です。
さて、こうして求められた $(s,t)$ の値から、もともと欲しかった式(1)の解 $(x,y)$ を構築できるのでした。その手順はまず
$$
x=t
$$
さらにその $x$ を利用して
$$
y = s-qx
$$
ですので、その通りに書きます。（xとyを参照渡ししているので、extGCDの中でxとyを書き換えれば外にあるxとyを変更できます（雑な説明））
long long extGCD(long long a, long long b, long long &x, long long &y) {
long long q = a/b;
long long r = a%b;
long long s, t;
long long d = extGCD(b, r, s, t);
x = t;
y = s - q * x;
}
ついでに、 $d$ も返してあげます。if(b==0)の処理も書いてあげます。
long long extGCD(long long a, long long b, long long &x, long long &y) {
if (b == 0) {
x = 1;
y = 0;
return a;
}
long long q = a/b;
long long r = a%b;
long long s, t;
long long d = extGCD(b, r, s, t);
x = t;
y = s - q * x;
return d;
}
これで完成です。当然これでも正しく動作します。
ただ、無駄があるのでそこを切り詰めていきます。それによって最初に紹介した、素人だとぱっと見でわからないプログラム（悪く言っているわけではないです）ができます。
まず、qとrはわざわざ変数を作る必要が無いのでそれぞれ書き換えます。
long long extGCD(long long a, long long b, long long &x, long long &y) {
if (b == 0) {
x = 1;
y = 0;
return a;
}
long long s, t;
long long d = extGCD(b, a%b, s, t);
x = t;
y = s - a/b * x;
return d;
}
sとtの箱もいらないです。まず引数のtをxに変えても同じような動作になることがわかります。
long long extGCD(long long a, long long b, long long &x, long long &y) {
if (b == 0) {
x = 1;
y = 0;
return a;
}
long long s;
long long d = extGCD(b, a%b, s, x);
y = s - a/b * x;
return d;
}
内部にあるextGCDが直接 $t$ の値を x に入れてくれます。さらに、sをyに変えます。
long long extGCD(long long a, long long b, long long &x, long long &y) {
if (b == 0) {
x = 1;
y = 0;
return a;
}
long long d = extGCD(b, a%b, y, x);
y = y - a/b * x;
return d;
}
ここで、y = y - a/b * x;の右辺のy,xには、いわゆる $s,t$ の値が入っています。それを利用して左辺のyを、外側のextGCDが返したい値に書き換えます。
これは複合代入演算子を使って記述できるので
long long extGCD(long long a, long long b, long long &x, long long &y) {
if (b == 0) {
x = 1;
y = 0;
return a;
}
long long d = extGCD(b, a%b, y, x);
y -= a/b * x;
return d;
}
と書けます。これで最初のプログラムと同じになりました！
おわり
ナニやってんだこれ？と思っていたプログラムが解読できてスッキリしました。「拡張ユークリッドの互除法」に辿り着く人でもし、ナンダこれ？となった人がいて、その方の役に立てれば幸いです。
おしまい。