daruma3940の日記

理解や文章に間違い等あればどんなことでもご指摘お願いします

拡張ユークリッド互除法について

f:id:daruma3940:20160520223745p:plain
ユークリッド互除法について書くのじぇ
ユークリッドの互除法について書くといったけどユークリッドの互除法についてはさらっと流して
拡張ユークリッド互除法について書きたいのじぇ
蟻本の拡張ユークリッドについての記述がよくわからなかったので自分でわかるようにした備忘録なのじぇ
f:id:daruma3940:20160521003616p:plain
まりちゃがこういうこといいだすのは大体よくわかってないことについてなので鵜呑みにしないほうがいいわよ
f:id:daruma3940:20160520223745p:plain

まずはユークリッド互除法について
自然数a,bに対して最大公約数を求める方法をユークリッド互除法というのじぇ
自然数a,bの最大公約数と毎回書くのは面倒なのでgcd(a,b)と書くことにするのじぇ。
ここでaをbで割った商をp余りをqとすると
a=b*p+qと書けることは明らかなのじぇ
b,qの最大公約数をgcd(b,q)とすると
b=x*gcd(b,q),q=y*gcd(b,q)のようにかけるので
a=gcd(b,q)*(x*p+y)と変形できるののじぇ。
gcd(b,q)は明らかにbを割り切り、上の式からaも割り切るのでgcd(b,q)はgcd(a,b)も割り切ることができなければならないのじぇ

ここで逆に
q=a-b*qとすると先ほどと同様にgcd(a,b)はgcd(b,q)を割り切れなければならないのじぇ
つまり
gcd(a,b)=X*gcd(b,q)
gcd(b,q)=Y*gcd(a,b)
となり、最大公約数は少なくとも1(>0)より大きいはずなのでX=Y=1となり
gcd(a,b)=gcd(b,q)となるのじぇ

aをbで割った余りをqとしたのでC言語的にq=a%bと書くとするのじぇ.
上の処理を何回も繰り返すと第二引数はどんどん減っていき最終的にゼロとなるのじぇ
gcd(a,b)=gcd(b,a%b)=gcd(a%b,b%(a%b))=.......=gcd(x,0)
0とxの最小公倍数はxなのでこれでgcd(a,b)は求まるのじぇこれがユークリッドの互除法なのじぇ

コードにするとこんな感じ(from 蟻本)
aをbにbをa%bにどんどん変えていく感じ

int gcd(int a,int b){
    if(b==0){return a;}
    return gcd(b,a%b);
}

f:id:daruma3940:20160520223745p:plain
これを拡張した拡張ユークリッドの互除法を利用することで
ax+by=gcd(a,b)となる整数解x,yをを求めることができるのじぇ

ちなみに右辺がなんでgcd(a,b)じゃないといけないのかについてはここを見てほしいのじぇ
mathtrain.jp

f:id:daruma3940:20160520223745p:plain
ユークリッド互除法のaをbに、bをa%bにどんどん変えていってどんどん係数を簡単にしていってgcd(a,b)を求めるのに似た方法を用いて
ax+by=gcd(a,b)もどんどん簡単にしていくのじぇ
f:id:daruma3940:20160520223745p:plain
最大公約数を求める方法と似てるなーーと思いながら読み進めてほしいのじぇ
f:id:daruma3940:20160520223745p:plain
とりあえず一回簡単にしてみるのじぇ
bx'+(a%b)y'=gcd(a,b)

aがbになってbがa%bになってxがx'にyがy'になっているのじぇ
gcd(b,a%b)=gcd(a,b)なので右辺は変わらんのじぇ

さっきのと似てるなーーーー


そしてこれをどんどん繰り返していくのじぇ
そうするとyの係数bはどんどん減少していって最終的にゼロになるのじぇ

ax'+0y'=gcd(a,b)

gcd(a,b)を求めるときbがゼロになったらgcd(a,b)=aだったのを思い出してほしいのじぇ
ということはx'の値は1でないといけないのじぇ

a*1+0*y=gcd(a,b)

0がかかっているのでy'の値はなんでもいいのじぇ。
ここではy'=0ということにしておくのじぇ
f:id:daruma3940:20160520223745p:plain
さてここまでやってきたことはax+by=gcd(a,b)をどんどん簡単にしていくことだったのじぇ

(a,b,x,y)

(b,a%b,x',y')
.
.
.

(a,0,1,0)
f:id:daruma3940:20160520223745p:plain
ax+by=gcd(a,b)となる整数解x,yを求めるには今度はこの逆をたどって最初の式に戻らなければならないのじぇ

(a,b,x,y)
↓  
(b,a%b,x',y')
.
.
.
↓  
(a,0,1,0)

f:id:daruma3940:20160520223745p:plain
その方法をここで説明するのじぇ

ax+by=gcd(a,b)

bx'+(a%b)y'=gcd(a,b)
にした場合のもとに戻し方を考えるのじぇ。

a%bはプログラミングの整数の割り算の扱い的に
a%b=a-(a/b)*b
とかけるので
bx'+(a%b)y'=gcd(a,b)は
bx'+(a-(a/b)*b)y'=gcd(a,b)となって
aとbについてまとめると
ay'+b(x'-(a/b)*y')=gcd(a,b)
となるのじぇ。
これを
ax+by=gcd(a,b)
と比較すると
x=y'
y=x'-(a/b)*y'

となるのでこのようにしてx,yの値を順にもとに戻していくことで
ax+by=gcd(a,b)となるx,yを求めることができるのじぇ
f:id:daruma3940:20160520223745p:plain
コードにするとこんな感じ(from 蟻本)

int extgcd(int a,int b,int& x,int& y){
    int d=a;//b==0のときgcd=a
    if(b!=0){
        //ここでgcdが求まる(b=0になる)まで潜る
        //さかのぼりながらx yを更新してゆく
        d=extgcd(b,a%b,y,x);

        //y=x'-(a/b)y'
        //xは1つ前のy yは一つ前のx
        y-=(a/b)*x;
    }
    else{
        //gcdが求まるとき(b==0のとき) x=1,y=0にする
        x=1;y=0;
    }
    return d;//gcdを返す
}


f:id:daruma3940:20160520223745p:plain

yの値をなんで0にしたのかについてはこれはほんとに何でもよくて他の値にすれば
他のax+by=gcd(a,b)を満たすx,yが求まるだけなのじぇ。
そもそもax+by=gcd(a,b)の解は一意ではないので...

f:id:daruma3940:20160520223745p:plain

ここからは完全に余談なのじぇ

さっき紹介したサイトを見ると
ax+by=gcd(a,b)の解き方は

この方程式の解を一つ求めてそれを(x0,y0)として

a(x1)+b(x2)=0
(x1:=x-x0 , y1:=y-y0 とした)

の形にすると
これは線形方程式なのでこの一般解はすぐに求まる。
そこからx=x1+x0,y=y1+y2の関係を用いて
x,yの一般式を求めることができると書いてあるのじぇ。
ここで(x0,y0)はさっき述べた拡張ユークリッド互除法で求めることができるのでこれで一般解を求めることができるのじぇ


f:id:daruma3940:20160520223745p:plain
これは非同次微分方程式を解くときと同じ考え方なのじぇ。

y''+ay'+y=f(x)
のを満たす解を一つ見つけてy_1とすると
(y-y_1)''+a(y-y_1)'+(y-y_1)=0 という形にできて
これは線形微分方程式なので簡単にとけて、その解をy_2とすると
y=y1+y2の関係から非同次微分方程式の一般解も求まるのじぇ

f:id:daruma3940:20160520223745p:plain
なるほど~~~~って感じなのじぇ