拡張ユークリッドの互除法を2つの実装から理解する

公開:2020/10/10
更新:2020/10/10
8 min読了の目安(約8000字TECH技術記事

拡張ユークリッドの互除法は逆元の計算に利用できるアルゴリズムですね。逆元を求めれば、modmod上で割り算をすることができます。この記事では逆元を計算することにフォーカスして、2種類の実装から拡張ユークリッドの互除法を理解することを目指します!ここの記事ではユークリッドの互除法については解説していません。

拡張ユークリッドの互除法

modp\mod p上での整数nnに対する逆元は以下のように表すことができます。ppは素数とします。

n1n=1modp n^{-1} * n=1\mod p

逆元[1]の定義はnnにかけたとき、11になる数なのでこれは明らかです。次に右辺を移項します。

n1n1=0modp n^{-1}* n-1=0\mod p

これにより、nn11n* n^{-1}-1ppで割り切れることになります。つまり任意の整数mmを使って

n1n1=mp n^{-1}* n -1=m* p

と表すことができます。これを式変形すると

n1n+mp=1 n^{-1} * n +m* p=1

となります。拡張ユークリッドの互除法は


a,ba,bを0でない整数とする。以下を満たすx,yx,yが存在し、計算することが出来る。

xa+yb=GCD(a,b) x*a+y*b=GCD(a,b)

というアルゴリズムでした。n1n+mp=1n^{-1} * n+m* p=1を見ると、GCD(n,p)=1GCD(n,p)=1ならn1n^{-1}mmを求めるのに利用できそうです。つまり、拡張ユークリッドの互除法で逆元が計算できます


n,pn,pを0でなく互いに素な整数とする。以下を満たすn1,mn^{-1},mが存在し、計算することが出来る。

n1n+mp=1 n^{-1} *n +m*p=1

以下ではこのn1,mn^{-1},mを計算するアルゴリズムを実装します!以下では、xa+yb=GCD(a,b)x*a+y*b=GCD(a,b)の表記を倣ってプログラムを書いています。例えば、x=n1x=n^{-1}y=my=mです。

表記について

以下のような表記が登場します。

aixi+biyi=1 a_i*x_i+b_i*y_i=1

この表記は変形の過程で現れる i番目の方程式 を表しています。読んでください。例えばi=1i=1なら

a1x1+b1y1=1 a_1*x_1+b_1*y_1=1

です。なお、iiは0から開始します。

その1 非再帰

まずは非再帰でユークリッドの互除法を実装します。飛ばして再帰の実装まで進んでも構いません。プログラムから見ていきます。

pair<int,int> extGCD1(int a,int b){
    int x=1;
    int y=0;
    int nx=0;
    int ny=1;
    while(a%b!=0){
        int q=a/b;
        int r=a%b;
        int tmpx=x-q*nx;
        int tmpy=y-q*ny;

        a=b;b=r;
        x=nx;y=ny;
        nx=tmpx;ny=tmpy;
    }
    return pair<int,int>(nx,ny);
}

何やら1とか0とか出てきてよくわからないかもしれません。これはaabbが次のように表せる事実を利用しています。

1a+0b=a0a+1b=b \begin{aligned} 1*a+0*b=a\\ 0*a+1*b=b \end{aligned}

なぜこのような形にするのかといえば、拡張ユークリッドの式xa+yb=GCD(a,b)x*a+y*b=GCD(a,b)に形を合わせるためです。この方法では右辺ではユークリッドの互除法を行い、左辺は等式が成立するように式を変形させます。変形を繰り返して式が

GCD(a,b)=xa+yb GCD(a,b)=x*a+y*b

となったときのxxyyが答えです。変形方法は以下のような式に表せます。


{x0a+y0b=r0 (i=0)x1a+y1b=r1 (i=1)(xi2xi1ri2ri1)a+(yi2yi1ri2ri1)b=ri2modri1=ri (2i) \begin{cases} x_0*a+y_0*b=r_0\ & (i=0)\\ x_1*a+y_1*b=r_1\ & (i=1)\\ (x_{i-2}-x_{i-1}*\lfloor \frac{r_{i-2}}{r_{i-1}}\rfloor)*a+(y_{i-2}-y_{i-1}*\lfloor \frac{r_{i-2}}{r_{i-1}}\rfloor)*b=r_{i-2}\mod r_{i-1}=r_i\ & (2\leq i) \end{cases}

(2i)(2\leq i)の式が何やら複雑ですが、ここではまだわからなくても問題ありません。[2]
具体例をあげましょう。a=97,b=53a=97,b=53のとき0番目の方程式と1番目の方程式は

197+053=97097+153=53 \begin{aligned} &1*97+0*53=97\\ &0*97+1*53=53\\ \end{aligned}

です。2番目を求めるためには(2i)(2\leq i)の式を利用します。

r2=(x0x1r1r0)a+(y0y1r1r0)b=(109753)97+(019753)53=197+(1)53=44 \begin{aligned} r_2&=(x_{0}-x_{1}*\lfloor \frac{r_{1}}{r_{0}}\rfloor)*a+(y_{0}-y_{1}*\lfloor \frac{r_{1}}{r_{0}}\rfloor)*b\\ &=(1-0*\lfloor \frac{97}{53}\rfloor)*97+(0-1*\lfloor \frac{97}{53}\rfloor)*53\\[10pt] &=1*97+(-1)*53=44 \end{aligned}

となります。これは97mod53=4497\mod 53=44xa+ybx*a+y*bの形で表すことができたことになります!
もう一度やってみましょう。今回は3番目を求めます。変わらずa=97,b=53a=97,b=53なので

r3=(x1x2r2r1)a+(y1y2r2r1)b=(015344)97+(1(1)5344)53=(1)97+253=9 \begin{aligned} r_3&=(x_{1}-x_{2}*\lfloor \frac{r_{2}}{r_{1}}\rfloor)*a+(y_{1}-y_{2}*\lfloor \frac{r_{2}}{r_{1}}\rfloor)*b\\ &=(0-1*\lfloor \frac{53}{44}\rfloor)*97+(1-(-1)*\lfloor \frac{53}{44}\rfloor)*53\\[10pt] &=(-1)*97+2*53=9 \end{aligned}

となります。これは53mod44=953\mod 44=9xa+ybx*a+y*bの形で表すことができたことになります!
これをri=0r_i=0まで繰り返します。

r4=597+(9)53=(1(1)4)97+(124)53=8r5=697+1153=(151)97+(2(9)1)53=1r6=0 \begin{aligned} r_4&=5*97+(-9)*53=(1-(-1)*4)*97+(-1-2*4)*53=8\\ r_5&=-6*97+11*53=(-1-5*1)*97+(2-(-9)*1)*53=1\\ r_6&=0 \end{aligned}

6番目でr4modc5=0r_4\mod c_5=0で割り切れているので、5番目のときの右辺11GCD(97,53)GCD(97,53)であり、(x5,y5)=(6,11)(x_5,y_5)=(-6,11)が拡張ユークリッドの互除法の答えになります。

実際に97(6)mod53=197*(-6)\mod 53=1なので逆元となっていますね。

i<=2の式について

(2i)(2\leq i)の式について解説します。最初、aabbは次のように表しました。

1a+0b=a0a+1b=b \begin{aligned} 1*a+0*b=a\\ 0*a+1*b=b \end{aligned}

amodba\mod baabbの商q=abq=\lfloor \frac{a}{b}\rfloorを使って次のように表せます。

amodb=aqb a\mod b=a-q*b

これは上述の式を使うと

aqb=1a+0b(0a+1b)q a-q*b=1*a+0*b-(0*a+1*b)*q

これをaとbごとにくくると、

1a+0b(0a+1b)q=(10q)a+(01q)b 1*a+0*b-(0*a+1*b)*q=(1-0*q)*a+(0-1*q)*b

と表せます。つまり、2番目は0番目と1番目から計算できます。同様に3番目を1番目と2番目から計算し....を繰り返して変形し続ければ、ii番目をi2i-2番目とi1i-1番目から求める時、ri2modri1=0r_{i-2}\mod r_{i-1}=0ならri1r_{i-1}の右辺がGCD(a,b)GCD(a,b)となっています。

上の式を一般化することを考えます。r0=a,r1=br_0=a,r_1=bとし

(10q)a+(01q)b=(x0x1q)a+(y0y1q)b (1-0*q)*a+(0-1*q)*b=(x_0-x_1*q)*a+(y_0-y_1*q)*b

とします。q=ri2ri1q=\lfloor \frac{r_{i-2}}{r_{i-1}}\rfloorなので


(xi2xi1ri2ri1)a+(yi2yi1ri2ri1)b=ri2modri1=ri (x_{i-2}-x_{i-1}*\lfloor \frac{r_{i-2}}{r_{i-1}}\rfloor)*a+(y_{i-2}-y_{i-1}*\lfloor \frac{r_{i-2}}{r_{i-1}}\rfloor)*b=r_{i-2}\mod r_{i-1}=r_i

と一般化できます。これは(2i)(2\leq i)の式です。

実装との対応

プログラムは次のような内容でした。

pair<int,int> extGCD1(int a,int b){
    int x=1;
    int y=0;
    int nx=0;
    int ny=1;
    while(a%b!=0){
        int q=a/b;
        int r=a%b;
        int tmpx=x-q*nx;
        int tmpy=y-q*ny;

        a=b;b=r;
        x=nx;y=ny;
        nx=tmpx;ny=tmpy;
    }
    return pair<int,int>(nx,ny);
}

次の4行は

int x=1;
int y=0;
int nx=0;
int ny=1;

次の初期条件を設定しています。

{x0a+y0b=r0 (i=0)x1a+y1b=r1 (i=1) \begin{cases} x_0*a+y_0*b=r_0\ & (i=0)\\ x_1*a+y_1*b=r_1\ & (i=1)\\ \end{cases}

次の4行は

int q=a/b;
int r=a%b;
int tmpx=x-q*nx;
int tmpy=y-q*ny;

次の式と対応しています。

ri2ri1ri2modri1xi2xi1ri2ri1yi2yi1ri2ri1 \begin{aligned} &\lfloor \frac{r_{i-2}}{r_{i-1}}\rfloor\\ &r_{i-2}\mod r_{i-1}\\ &x_{i-2}-x_{i-1}*\lfloor \frac{r_{i-2}}{r_{i-1}}\rfloor\\ &y_{i-2}-y_{i-1}*\lfloor \frac{r_{i-2}}{r_{i-1}}\rfloor \end{aligned}

その2 再帰

次は再帰でユークリッドの互除法を書いてみます。次のようになります。

pair<int,int> extGCD2(int a,int b) {
    if (b == 0) return pair<int,int>(1,0);
    pair<int,int> d=extGCD2(b, a%b);
    int q=a/b;
    int x=d.second;
    int y=d.first-q*d.second;
    return pair<int,int>(x,y);
}

また、a0=97a_0=97b0=53b_0=53の例でやってみます。まずはやることは、a=ba=b,b=amodbb=a\mod bを繰り返すだけです。ユークリッドの互除法と同じですね。

x097+y053=GCD(97,53) x_0*97+y_0*53=GCD(97,53)

このとき、

a=53b=97mod53=44 \begin{aligned} a&=53\\ b&=97\mod 53=44 \end{aligned}

となります。よって、1番目の方程式は

x53+y44=1 \begin{aligned} x*53+y*44=1\\ \end{aligned}

です。これをb=0b=0まで繰り返します。

x53+y44=1x44+y9=1x9+y8=1x8+y1=1x1+y0=1 \begin{aligned} &x*53+y*44=1\\ &x*44+y*9=1\\ &x*9+y*8=1\\ &x*8+y*1=1\\ &x*1+y*0=1 \end{aligned}

5番目の方程式で、b=0b=0になりました。ここで再帰が終わるので解を復元します。まず、1つ前の4番目の方程式のx4x_4y4y_4を求めることを考えます。ii番目のx,yx,yをそれぞれxi,yix_i,y_iとすると以下で求まります。

x4=y5=0y4=x5y5q=1081=1 \begin{aligned} x_4&=y_5=0\\ y_4&=x_5-y_5*q=1-0*\frac{8}{1}=1 \end{aligned}

これを一般化すると次のようになります。

xi=yi+1yi=xi+1yi+1aibi \begin{aligned} x_i&=y_{i+1}\\ y_i&=x_{i+1}-y_{i+1}*\lfloor \frac{a_i}{b_i}\rfloor \end{aligned}

適用を繰り返すことで、与式97x0+53y0197*x_0+53*y_0*1x0x_0y0y_0を求められそうです。実際に適用してみると

x4=y5=0y4=x5y5a4b4=1081=1x3=y4=1y3=x4y4a3b3=0198=1x2=y3=1y2=x3y3a2b2=1(1)449=5x1=y2=5y1=x2y2a1b1=155344=6x0=y1=6y0=x1y1a0b0=5(6)9753=11 \begin{aligned} x_4&=y_5=0\\ y_4&=x_5-y_5*\lfloor \frac{a_4}{b_4}\rfloor=1-0*\lfloor \frac{8}{1}\rfloor=1\\ \\[3pt] x_3&=y_4=1\\ y_3&=x_4-y_4*\lfloor \frac{a_3}{b_3}\rfloor=0-1*\lfloor \frac{9}{8}\rfloor=-1\\ \\[3pt] x_2&=y_3=-1\\ y_2&=x_3-y_3*\lfloor \frac{a_2}{b_2}\rfloor=1-(-1)*\lfloor \frac{44}{9}\rfloor=5\\ \\[3pt] x_1&=y_2=5\\ y_1&=x_2-y_2*\lfloor \frac{a_1}{b_1}\rfloor=-1-5*\lfloor \frac{53}{44}\rfloor=-6\\ \\[3pt] x_0&=y_1=-6\\ y_0&=x_1-y_1*\lfloor \frac{a_0}{b_0}\rfloor=-5-(-6)*\lfloor \frac{97}{53}\rfloor=11\\ \\[3pt] \end{aligned}

(x,y)=(6,11)(x,y)=(-6,11)です。実際に97(6)mod53=197*(-6)\mod 53=1なので逆元となっていますね。

式の解説

xi=yi+1yi=xi+1aibiyi+1 \begin{aligned} x_i&=y_{i+1}\\ y_i&=x_{i+1}-\lfloor \frac{a_i}{b_i}\rfloor*y_{i+1} \end{aligned}

この式はどこから出てきたのでしょう。i=0i=0の例を考えます。

x0a0+y0b0=1 x_0*a_0+y_0*b_0=1

ここでa0=qb0+ra_0=q*b_0+rとして、a0a_0を使わずに等式を成立させることを考えます。ここでq=a0/b0,r=a0qb0q=\lfloor a_0/b_0 \rfloor,r=a_0-q*b_0です。つまり、商と余りです。

x0(qb0+r)+y0b0=1 x_0*(q*b_0+r)+y_0*b_0= 1

です。このとき

x0(qb0+r)+y0b0=b0(x0q+y0)+rx0 x_0*(q*b_0+r)+y_0*b_0=b_0*(x_0*q+y_0)+r*x_0

この式をc1c_1として書き換えると

b0(x0q+y0)+rx0=a1x1+b1y1=1 b_0*(x_0*q+y_0)+r*x_0=a_1*x_1+b_1*y_1=1

ですね!これにより、

x1=x0q+y0y1=x0 \begin{aligned} x_1&=x_0*q+y_0\\ y_1&=x_0 \end{aligned}

であったことが分かります。これをx0,y0x_0,y_0についての方程式にすると

y0=x1x0qx0=y1 \begin{aligned} y_0&=x_1-x_0*q\\ x_0&=y_1 \end{aligned}

ですね。結局x0=y1x_0=y_1なので、

x0=y1y0=x1y1q \begin{aligned} x_0&=y_1\\ y_0&=x_1-y_1*q \end{aligned}

であり、これは一般化した式そのものです。さらにa1=b0,a1=a0modb0a_1=b_0,a_1=a_0\mod b_0になっていることも分かります。

実装との対応

プログラムは次のような内容でした。

pair<int,int> extGCD2(int a,int b) {
    if (b == 0) return pair<int,int>(1,0);
    pair<int,int> d=extGCD2(b, a%b);
    int q=a/b;
    int x=d.second;
    int y=d.first-q*d.second;
    return pair<int,int>(x,y);
}

まず、次の部分はa=ba=b,b=amodbb=a\mod bを繰り返して、b=0b=0で再帰を打ち切っています。b=0b=0まではこの2行しか実行されません。

if (b == 0) return pair<int,int>(1,0);
pair<int,int> d=extGCD2(b, a%b);

次の部分は

int q=a/b;
int x=d.second;
int y=d.first-q*d.second;
return pair<int,int>(x,y);

以下と対応しています。

aibixi=yi+1yi=xi+1aibiyi+1 \begin{aligned} &\lfloor\frac{a_i}{b_i}\rfloor\\ &x_i=y_{i+1}\\ &y_i=x_{i+1}-\lfloor \frac{a_i}{b_i}\rfloor*y_{i+1}\\ \end{aligned}

まとめ

ユークリッドの互除法は比較的理解しやすいアルゴリズムだと思うのですが、拡張ユークリッドの互除法はアルゴリズムとして記述すると複雑になりがちだと思います。この記事を通して自分のアルゴリズムへの理解が随分深まりました。

間違いや質問がありましたらコメントまでよろしくおねがいします!

脚注
  1. 正確には乗法逆元です。 ↩︎

  2. 実際には常にx0=1,y0=0,x1=0,y1=1x_0=1,y_0=0,x_1=0,y_1=1という定数ですが、数式にする都合で定数にはしていません。 ↩︎