ユークリッドの互除法

ユークリッドの互除法の計算量

gcd(a,b)としたときに $Θ(log max(a,b) )$
より簡単にN=max(a,b)とすると、 $Θ(logN)$ です。

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

なぜlogが出てくるのか

例えば、計算量にlogが出てくるものといえば二分探索があります。区間を半分、半分、半分、...みたいにしていくものは大体logがつくと思っています。これは長さNの区間を1/2倍,1/2倍,1/2倍...としていくときに1になるまで1/2倍する回数が logN 回くらいだからです。

なんでlogN回くらいなのか

これは 1をN以上にするまでに何回2倍するかと考えるとわかりやすくなると思います。

1 → 2 → 4 → 8 →...
と2倍していくといずれNを超えると思います。

つまり、

2 ^x ≥ N

となる x を求めたいわけです。

ところで、logの定義は覚えていますか？ logというのは

a^b= c

というときに、aのなんちゃら乗がcとなるようなbを求めるものです。
例えば、a=2,c=8だったときに

2^b=8

となるようなbを、logを使って求めることができます。

b=log₂8

といった感じです。この場合のbは

b=log₂2³

b=3*log₂2

b=3

となって、bが3と求まります。

logを復習したところでさっきの式を見てみましょう。

2 ^x ≥ N

この x は

x ≥ log₂ N

こんな風に求めることができます。つまり、1を2倍、2倍、2倍、...という操作をlog₂N回行えばNを超える数になります。計算量としてはlog₂N回なのですがわかりやすい形に整えるため、底の変換公式を利用します。

log₂N

底の変換公式を使うと、

log₂N = ${\frac{log_e N}{log_e 2}}$

log₂N = ${\frac{1}{log_e 2}\times log_e N}$

とすることができ、定数をくくりだせました。¹

これで1をNが超えるまで2倍、2倍、2倍...としていったときの操作回数はlogN回であることが分かったので、Nを1/2倍、1/2倍、1/2倍...としていったときの操作回数もlogN回です。

話を戻して...

ではユークリッドの互除法では本当にlogが出てくるのでしょうか。

ユークリッドの互除法では

(a, b)の最大公約数は(b, a%b)の最大公約数と等しい という性質を利用しています。²

ここで2つほど疑問が生まれます。

どうしてこの性質が成り立つのか
この性質が成り立ったからどうなんだ

どうしてこの性質が成り立つのかというのは、ひとまずあとに置いておきます。先にこの性質が成り立つことによって計算量がlogになることを示します。

ユークリッドの互除法の計算量がlogになる説明。

前提として、(a, b)は a≥b であるものとします。毎回max(a,b)とするのが面倒だからです。

そして、aをbで割ったときの商を q , 余りを r と文字におきます。
(ちなみに、商(quotient),余り(remainder)の頭文字から来ています)

すると、

a = qb + r

という式でaを表すことができます。

a≥b という条件から、aの中にはbは必ず1つ以上あるので q≥1 になります。よって、先ほどの式からqを抜いた式を考えると

a ≥ b + r

のようになります。

また、商は余りよりも大きいため

b > r

これを変形すると、

b + r > 2r

a ≥ b + r
b + r > 2r

の二つをくっつけて、

a ≥ b + r > 2r

真ん中をなくすと、

a ≥ 2r $\frac{a}{2} ≥ r$

という関係式が成り立ちます。このことから、計算をするたびにざっくり半分以下、半分以下、半分以下、...となっていき、計算回数はおよそΘ(loga)回になります。

では、先ほど後回しにした

(a, b)の最大公約数は(b, a%b)の最大公約数と等しい が成り立つ理由を説明していきます。

(a, b)の最大公約数を G
(b, a%b)の最大公約数を g
という文字におきます。 G = gであることを示せれば勝ちです。

aをbで割ったときの商を q , 余りを r と文字におきます。

先ほどのGとgはこのように言い換えることもできて、
(a, b)の最大公約数を G
(b, r)の最大公約数を g

a=qb+r

という式でaを表せます。この式の右辺はbとrから成り立っているのでgで割り切ることができます。よって、左辺であるaも割り切ることができるので、gはaの約数です。

このことからaとbの公約数がgであることが分かります。その中でも最大の公約数がGであるので

G ≥ g

が成り立ちます。

次に、

a-qb=r

という先ほどの式を変形したものを考えます。この式の左辺はaとbから成り立っているのでGで割り切ることができます。
よって、右辺であるrも割り切ることができるので、Gはrの約数です。

このことからbとrの公約数がGであることが分かります。その中でも最大の公約数がgであるので

g ≥ G

が成り立ちます。

G ≥ g
g ≥ G

この2の式から

G=g

が成り立ちます。

これで、

(a, b)の最大公約数は(b, a%b)の最大公約数と等しい この性質が正しいことであることが分かりました。

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

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

これもユークリッドの互除法と同じでextgcd(a,b,x,y)としたときに $Θ(log max(a,b) )$

int extgcd(int a, int b, int &x, int &y){
    if(b==0){
        x=1;
        y=0;
        return a;
    }
    int q = a/b;
    int r = a%b;
    int s, t;
    int d = extgcd(b, r, s, t)
    x = t;
    y = s - q * x;
    return d;
}

こちらの記事を見てもらうとより理解が深まると思います。
拡張ユークリッドの互除法〜一次不定方程式 ax + by = c の解き方〜

拡張ユークリッドの互除法はなにをするのか？

ユークリッドの互除法はaとbを与えるとその最大公約数を求めるものでした。
拡張ユークリッドの互除法では、aとbの最大公約数に加えて次のようなことを求めてくれます。

ax + by = gcd(a, b)となるxとy

ぱっと見てしっかり理解するのは難しいので次のような例を見ていきます。

111x + 30y = 12

まず、gcd(111, 30)は3です。
拡張ユークリッドの互除法で求めることができるのがax+by=gcd(a,b)のxとyなので、111x+30y=3を先に求めます。

int x,y;
extgcd(111, 30, x, y);

こうするとxとyという変数に条件を満たす数がそれぞれ代入されます。
今回の例では

x=3,y=-11

と求まります。

では、extgcdの内部では何をしているのでしょうか。

extgcdのやっていること

一言で言えば、式変形をプログラムで動かしています。

を求めるために式変形をします。

a = qb + r ax + by = gcd(a,b) (qb + r)x + by = gcd(a,b) qbx + rx + by = gcd(a,b) (qx + y)b +rx = gcd(a,b)

s=qx + y

t=x

というようにs,tとおくことで

bs + rt = gcd(a,b)

のように表現できます。

次に、sとtでxとyを表します。

x=t

s=qx+y s-qx=y y=s-qx y=s-qt

s,tが求りさえすればxとyが求まって、それぞれ値を代入できます。

s,tは求まるのか

結論から言うと、s,tは求まります。
式変形やプログラムをよく見てみましょう。

式変形
- (qx + y)b +rx = gcd(a,b)
プログラム

int extgcd(int a, int b, int &x, int &y){
    if(b==0){
        x=1;
        y=0;
        return a;
    }
    int q = a/b;
    int r = a%b;
    int s, t;
    int d = extgcd(b, r, s, t)
    x = t;
    y = s - q * x;
    return d;
}

注目してほしいのはaとbの移り変わりです。

(a,b) → (b,a%b)

という風に変化しています。これはまさしくユークリッドの互除法です。
そして、終了条件であるif(b==0)というのもユークリッドの互除法の終了条件と同じです。
ユークリッドの互除法でb==0となるときはaが最大公約数でした。

ax+by=gcd(a,b)

のbが0だったとき

ax=gcd(a,b)

aが最大公約数なので

gcd(a,b)x=gcd(a,b)

よって、 x=1

yはここで0としておきます。
数学的にはyはなんでも良く、b==0のときにy=100やy=-50などにしても問題なく求まります。
しかし、y=0とすることによって|x|+|y|が数ある解の中で最小のものになります。

よってb==0で終了し、s,tは求まります。

実装

これまで拡張ユークリッドの互除法について説明してきましたが、実装をしていきたいと思います。といっても今まで出てきたコードを少しまとめるだけです。

int extgcd(int a, int b, int &x, int &y){
    if(b==0){
        x=1;
        y=0;
        return a;
    }
    int q = a/b;
    int r = a%b;
    int s, t;
    int d = extgcd(b, r, s, t)
    x = t;
    y = s - q * x;
    return d;
}

qとrをa/b,a%bに置き換えます。

int extgcd(int a, int b, int &x, int &y){
    if(b==0){
        x=1;
        y=0;
        return a;
    }
    int s, t;
    int d = extgcd(b, a%b, s, t)
    x = t;
    y = s - (a/b) * x;
    return d;
}

s,tの変数をyとxに置き換えます。

int extgcd(int a, int b, int &x, int &y){
    if(b==0){
        x=1;
        y=0;
        return a;
    }
    int d = extgcd(b, a%b, y, x)
    y = y - (a/b) * x;
    return d;
}

これで完成です。

まとめ&あとがき

ユークリッドの互除法はΘ(log max(a,b))

上式でのeはネイピア数を表しています。↩
%は割った余りを表しています。↩