拡張ユークリッドの互除法は逆元の計算に利用できるアルゴリズムですね。逆元を求めれば、mod上で割り算をすることができます。この記事では逆元を計算することにフォーカスして、２種類の実装から拡張ユークリッドの互除法を理解することを目指します！ここの記事ではユークリッドの互除法については解説していません。
拡張ユークリッドの互除法
\mod p上での整数nに対する逆元は以下のように表すことができます。pは素数とします。
逆元の定義はnにかけたとき、1になる数なのでこれは明らかです。次に右辺を移項します。
これにより、n* n^{-1}-1がpで割り切れることになります。つまり任意の整数mを使って
と表すことができます。これを式変形すると
となります。拡張ユークリッドの互除法は
a,bを0でない整数とする。以下を満たすx,yが存在し、計算することが出来る。
というアルゴリズムでした。n^{-1} * n+m* p=1を見ると、GCD(n,p)=1ならn^{-1}とmを求めるのに利用できそうです。つまり、拡張ユークリッドの互除法で逆元が計算できます。
n,pを0でなく互いに素な整数とする。以下を満たすn^{-1},mが存在し、計算することが出来る。
以下ではこのn^{-1},mを計算するアルゴリズムを実装します！以下では、x*a+y*b=GCD(a,b)の表記を倣ってプログラムを書いています。例えば、x=n^{-1}でy=mです。
表記について
以下のような表記が登場します。
この表記は変形の過程で現れる i番目の方程式 を表しています。読んでください。例えばi=1なら
です。なお、iは０から開始します。
その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とか出てきてよくわからないかもしれません。これはaとbが次のように表せる事実を利用しています。
\begin{aligned}
1*a+0*b=a\\
0*a+1*b=b
\end{aligned}
なぜこのような形にするのかといえば、拡張ユークリッドの式x*a+y*b=GCD(a,b)に形を合わせるためです。この方法では右辺ではユークリッドの互除法を行い、左辺は等式が成立するように式を変形させます。変形を繰り返して式が
となったときのxとyが答えです。変形方法は以下のような式に表せます。
\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\leq i)の式が何やら複雑ですが、ここではまだわからなくても問題ありません。
具体例をあげましょう。a=97,b=53のとき0番目の方程式と１番目の方程式は
\begin{aligned}
&1*97+0*53=97\\
&0*97+1*53=53\\
\end{aligned}
です。2番目を求めるためには(2\leq i)の式を利用します。
\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}
となります。これは97\mod 53=44をx*a+y*bの形で表すことができたことになります！
もう一度やってみましょう。今回は3番目を求めます。変わらずa=97,b=53なので
\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}
となります。これは53\mod 44=9をx*a+y*bの形で表すことができたことになります！
これをr_i=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番目でr_4\mod c_5=0で割り切れているので、5番目のときの右辺1がGCD(97,53)であり、(x_5,y_5)=(-6,11)が拡張ユークリッドの互除法の答えになります。
実際に97*(-6)\mod 53=1なので逆元となっていますね。
i<=2の式について
(2\leq i)の式について解説します。最初、aとbは次のように表しました。
\begin{aligned}
1*a+0*b=a\\
0*a+1*b=b
\end{aligned}
a\mod bはaとbの商q=\lfloor \frac{a}{b}\rfloorを使って次のように表せます。
これは上述の式を使うと
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
と表せます。つまり、2番目は0番目と1番目から計算できます。同様に3番目を1番目と2番目から計算し....を繰り返して変形し続ければ、i番目をi-2番目とi-1番目から求める時、r_{i-2}\mod r_{i-1}=0ならr_{i-1}の右辺がGCD(a,b)となっています。
上の式を一般化することを考えます。r_0=a,r_1=bとし
(1-0*q)*a+(0-1*q)*b=(x_0-x_1*q)*a+(y_0-y_1*q)*b
とします。q=\lfloor \frac{r_{i-2}}{r_{i-1}}\rfloorなので
(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)の式です。
実装との対応
プログラムは次のような内容でした。
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;
次の初期条件を設定しています。
\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;
次の式と対応しています。
\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);
}
また、a_0=97とb_0=53の例でやってみます。まずはやることは、a=b,b=a\mod bを繰り返すだけです。ユークリッドの互除法と同じですね。
このとき、
\begin{aligned}
a&=53\\
b&=97\mod 53=44
\end{aligned}
となります。よって、１番目の方程式は
\begin{aligned}
x*53+y*44=1\\
\end{aligned}
です。これをb=0まで繰り返します。
\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=0になりました。ここで再帰が終わるので解を復元します。まず、１つ前の４番目の方程式のx_4とy_4を求めることを考えます。i番目のx,yをそれぞれx_i,y_iとすると以下で求まります。
\begin{aligned}
x_4&=y_5=0\\
y_4&=x_5-y_5*q=1-0*\frac{8}{1}=1
\end{aligned}
これを一般化すると次のようになります。
\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*x_0+53*y_0*1のx_0とy_0を求められそうです。実際に適用してみると
\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)です。実際に97*(-6)\mod 53=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=0の例を考えます。
ここでa_0=q*b_0+rとして、a_0を使わずに等式を成立させることを考えます。ここでq=\lfloor a_0/b_0 \rfloor,r=a_0-q*b_0です。つまり、商と余りです。
です。このとき
x_0*(q*b_0+r)+y_0*b_0=b_0*(x_0*q+y_0)+r*x_0
この式をc_1として書き換えると
b_0*(x_0*q+y_0)+r*x_0=a_1*x_1+b_1*y_1=1
ですね！これにより、
\begin{aligned}
x_1&=x_0*q+y_0\\
y_1&=x_0
\end{aligned}
であったことが分かります。これをx_0,y_0についての方程式にすると
\begin{aligned}
y_0&=x_1-x_0*q\\
x_0&=y_1
\end{aligned}
ですね。結局x_0=y_1なので、
\begin{aligned}
x_0&=y_1\\
y_0&=x_1-y_1*q
\end{aligned}
であり、これは一般化した式そのものです。さらにa_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=b,b=a\mod bを繰り返して、b=0で再帰を打ち切っています。b=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);
以下と対応しています。
\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}
まとめ
ユークリッドの互除法は比較的理解しやすいアルゴリズムだと思うのですが、拡張ユークリッドの互除法はアルゴリズムとして記述すると複雑になりがちだと思います。この記事を通して自分のアルゴリズムへの理解が随分深まりました。
間違いや質問がありましたらコメントまでよろしくおねがいします！
Discussion
こんにちは。とても勉強になりました。
めちゃくちゃ些細で釈迦に説法ではありますが、
冒頭の部分
m*p を左辺に持ってきているので -m * p が良さそうですかね。
m' =-m とおくのでも良いかもです。
そもそもmは任意の整数なのであれですが..