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

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

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

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

$$n−1∗n=1mod p n^{-1} * n=1\mod p$$

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

$$n−1∗n−1=0mod p n^{-1}* n-1=0\mod p$$

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

$$n−1∗n−1=m∗p n^{-1}* n -1=m* p$$

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

$$n−1∗n+m∗p=1 n^{-1} * n +m* p=1$$

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

---

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

$$x∗a+y∗b=GCD(a,b) x*a+y*b=GCD(a,b)$$

---

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

---

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

$$n−1∗n+m∗p=1 n^{-1} *n +m*p=1$$

---

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

表記について

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

$$ai∗xi+bi∗yi=1 a_i*x_i+b_i*y_i=1$$

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

$$a1∗x1+b1∗y1=1 a_1*x_1+b_1*y_1=1$$

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

その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とか出てきてよくわからないかもしれません。これは$aa$と$bb$が次のように表せる事実を利用しています。

$$1∗a+0∗b=a0∗a+1∗b=b \begin{aligned} 1*a+0*b=a\\ 0*a+1*b=b \end{aligned}$$

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

$$GCD(a,b)=x∗a+y∗b GCD(a,b)=x*a+y*b$$

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

---

$${x0∗a+y0∗b=r0 (i=0)x1∗a+y1∗b=r1 (i=1)(xi−2−xi−1∗⌊ri−2ri−1⌋)∗a+(yi−2−yi−1∗⌊ri−2ri−1⌋)∗b=ri−2mod ri−1=ri (2≤i) \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}$$

---

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

$$1∗97+0∗53=970∗97+1∗53=53 \begin{aligned} &1*97+0*53=97\\ &0*97+1*53=53\\ \end{aligned}$$

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

$$r2=(x0−x1∗⌊r1r0⌋)∗a+(y0−y1∗⌊r1r0⌋)∗b=(1−0∗⌊9753⌋)∗97+(0−1∗⌊9753⌋)∗53=1∗97+(−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}$$

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

$$r3=(x1−x2∗⌊r2r1⌋)∗a+(y1−y2∗⌊r2r1⌋)∗b=(0−1∗⌊5344⌋)∗97+(1−(−1)∗⌊5344⌋)∗53=(−1)∗97+2∗53=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}$$

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

$$r4=5∗97+(−9)∗53=(1−(−1)∗4)∗97+(−1−2∗4)∗53=8r5=−6∗97+11∗53=(−1−5∗1)∗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番目で$r4mod c5=0r_4\mod c_5=0$で割り切れているので、5番目のときの右辺$11$が$GCD(97,53)GCD(97,53)$であり、$(x5,y5)=(−6,11)(x_5,y_5)=(-6,11)$が拡張ユークリッドの互除法の答えになります。

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

i<=2の式について

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

$$1∗a+0∗b=a0∗a+1∗b=b \begin{aligned} 1*a+0*b=a\\ 0*a+1*b=b \end{aligned}$$

$amod ba\mod b$は$aa$と$bb$の商$q=⌊ab⌋q=\lfloor \frac{a}{b}\rfloor$を使って次のように表せます。

$$amod b=a−q∗b a\mod b=a-q*b$$

これは上述の式を使うと

$$a−q∗b=1∗a+0∗b−(0∗a+1∗b)∗q a-q*b=1*a+0*b-(0*a+1*b)*q$$

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

$$1∗a+0∗b−(0∗a+1∗b)∗q=(1−0∗q)∗a+(0−1∗q)∗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$番目を$i−2i-2$番目と$i−1i-1$番目から求める時、$ri−2mod ri−1=0r_{i-2}\mod r_{i-1}=0$なら$ri−1r_{i-1}$の右辺が$GCD(a,b)GCD(a,b)$となっています。

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

$$(1−0∗q)∗a+(0−1∗q)∗b=(x0−x1∗q)∗a+(y0−y1∗q)∗b (1-0*q)*a+(0-1*q)*b=(x_0-x_1*q)*a+(y_0-y_1*q)*b$$

とします。$q=⌊ri−2ri−1⌋q=\lfloor \frac{r_{i-2}}{r_{i-1}}\rfloor$なので

---

$$(xi−2−xi−1∗⌊ri−2ri−1⌋)∗a+(yi−2−yi−1∗⌊ri−2ri−1⌋)∗b=ri−2mod ri−1=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$$

---

と一般化できます。これは$(2≤i)(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);
}

次の４行は

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

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

$${x0∗a+y0∗b=r0 (i=0)x1∗a+y1∗b=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}$$

次の４行は

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

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

$$⌊ri−2ri−1⌋ri−2mod ri−1xi−2−xi−1∗⌊ri−2ri−1⌋yi−2−yi−1∗⌊ri−2ri−1⌋ \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=97$と$b0=53b_0=53$の例でやってみます。まずはやることは、$a=ba=b$,$b=amod bb=a\mod b$を繰り返すだけです。ユークリッドの互除法と同じですね。

$$x0∗97+y0∗53=GCD(97,53) x_0*97+y_0*53=GCD(97,53)$$

このとき、

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

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

$$x∗53+y∗44=1 \begin{aligned} x*53+y*44=1\\ \end{aligned}$$

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

$$x∗53+y∗44=1x∗44+y∗9=1x∗9+y∗8=1x∗8+y∗1=1x∗1+y∗0=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}$$

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

$$x4=y5=0y4=x5−y5∗q=1−0∗81=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+1−yi+1∗⌊aibi⌋ \begin{aligned} x_i&=y_{i+1}\\ y_i&=x_{i+1}-y_{i+1}*\lfloor \frac{a_i}{b_i}\rfloor \end{aligned}$$

適用を繰り返すことで、与式$97∗x0+53∗y0∗197*x_0+53*y_0*1$の$x0x_0$と$y0y_0$を求められそうです。実際に適用してみると

$$x4=y5=0y4=x5−y5∗⌊a4b4⌋=1−0∗⌊81⌋=1x3=y4=1y3=x4−y4∗⌊a3b3⌋=0−1∗⌊98⌋=−1x2=y3=−1y2=x3−y3∗⌊a2b2⌋=1−(−1)∗⌊449⌋=5x1=y2=5y1=x2−y2∗⌊a1b1⌋=−1−5∗⌊5344⌋=−6x0=y1=−6y0=x1−y1∗⌊a0b0⌋=−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)mod 53=197*(-6)\mod 53=1$なので逆元となっていますね。

式の解説

$$xi=yi+1yi=xi+1−⌊aibi⌋∗yi+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$の例を考えます。

$$x0∗a0+y0∗b0=1 x_0*a_0+y_0*b_0=1$$

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

$$x0∗(q∗b0+r)+y0∗b0=1 x_0*(q*b_0+r)+y_0*b_0= 1$$

です。このとき

$$x0∗(q∗b0+r)+y0∗b0=b0∗(x0∗q+y0)+r∗x0 x_0*(q*b_0+r)+y_0*b_0=b_0*(x_0*q+y_0)+r*x_0$$

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

$$b0∗(x0∗q+y0)+r∗x0=a1∗x1+b1∗y1=1 b_0*(x_0*q+y_0)+r*x_0=a_1*x_1+b_1*y_1=1$$

ですね！これにより、

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

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

$$y0=x1−x0∗qx0=y1 \begin{aligned} y_0&=x_1-x_0*q\\ x_0&=y_1 \end{aligned}$$

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

$$x0=y1y0=x1−y1∗q \begin{aligned} x_0&=y_1\\ y_0&=x_1-y_1*q \end{aligned}$$

であり、これは一般化した式そのものです。さらに$a1=b0,a1=a0mod b0a_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=amod bb=a\mod b$を繰り返して、$b=0b=0$で再帰を打ち切っています。$b=0b=0$まではこの２行しか実行されません。

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

以下と対応しています。

$$⌊aibi⌋xi=yi+1yi=xi+1−⌊aibi⌋∗yi+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$という定数ですが、数式にする都合で定数にはしていません。