拡張ユークリッドの互除法を解説してみた

先日、こんな記事を書いた。

拡張ユークリッドの互除法を自己流で実装してみた。 - いろはすの競プロ反省会

これだけではあまりに不親切なので、ここに書いたコードの解説をしていこうと思う。

※上の記事に書いたコードは自己流で書いたものなので、蟻本に載っている拡張ユークリッドの互除法のコードとは異なります。予めご了承ください。また、ユークリッドの互除法については既知とします。

まず、拡張ユークリッドの互除法とは何かというと、0でない整数a,bに対して次の式を満たす整数x,yの組を求めるアルゴリズムである。

$xa+yb=gcd(a,b)$ ・・・(1)

ここで、 $gcd(a,b)$ はa,bの最大公約数である。これはベズーの等式として知られる式であり、整数解が必ず存在する。

アルゴリズムの説明に入ろう。まず、上の式において $a=a_0, b=b_0, x=x_0, y=y_0$ とおく。

$x_0a_0+y_0b_0=gcd(a_0,b_0)$ ・・・(2)

ここで、

$a_0=(a_0/b_0)b_0+a_0\%b_0$ ・・・（３）

（便宜上、剰余の記号として'%'を用いた）であるから、(3)を(2)に代入し、

$(x_0(a_0/b_0)+y_0)b_0+x_0(a_0\%b_0)=gcd(a_0,b_0)$ ・・・(4)

を得る。ここで、

$x_1=x_0(a_0/b_0)+y_0,$ ・・・(5)

$y_1=x_0,$ ・・・(6)

$a_1=b_0,$ ・・・(7)

$b_1=a_0\%b_0$ ・・・(8)

とすると、式(4)は

$x_1a_1+y_1b_1=gcd(a_0,b_0)$ ・・・(9)

となる。ここで、式(7),(8)の変換をよく見ると、これはユークリッドの互除法で用いる変換と全く同じ変換である。従って、

$gcd(a_0,b_0)=gcd(a_1,b_1)$ ・・・(10)

これを(9)に代入し、

$x_1a_1+y_1b_1=gcd(a_1,b_1)$ ・・・(11)

を得る。

式(11)と式(2)を見比べると、何と全く同じ形をしている！　以下、同様に

$x_{i+1}=x_i(a_i/b_i)+y_i,$ ・・・(12)

$y_{i+1}=x_i,$ ・・・(13)

$a_{i+1}=b_i,$ ・・・(14)

$b_{i+1}=a_i\%b_i$ ・・・(15)

とすれば、次々に変換を続けていくことができる。

式(14),(15)の変換がユークリッドの互除法の変換と同じであることに注意すると、この操作はあるi=nで終了し、

$b_n=0,$

$a_n=gcd(a_n, b_n)$

となる。このとき、i=nに対応する等式は

$x_na_n=gcd(a_n, b_n)=a_n$

であるから、その解は（例えば）、

$x_n=1,$ ・・・(16)

$y_n=0$ ・・・(17)

である。あとは、ここから(12),(13)の逆変換

$x_i=y_{i+1},$ ・・・(18)

$y_i=x_{i+1}-y_{i+1}(a_i/b_i)$ ・・・(19)

を用いてこれまでの変換を逆に辿れば、元の式(2)の解が求まる。この変換式(18),(19)こそがアルゴリズムの肝であり、それを実装したのが以下のコードである。

//xa+yb=gcd(a,b)の整数解を求める
typedef pair<int,int> P;
P extgcd(int a, int b){
  if(b==0)
    return P(1, 0);

  P p = extgcd(b, a%b);
  int x = p.first;
  int y = p.second;

  return P(y, x - y*(a/b)); //逆変換(18),(19)
}