Web教材一覧>
アルゴリズム
(BOK大区分:1 基礎理論、中区分:2 アルゴリズムとプログラミング、中区分:2 アルゴリズム)
素数に関係するアルゴリズムで有名なユークリッドの互除法、それを拡張した拡張ユークリッドの互除法について学習します。
アルゴリズム、最大公約数、ユークリッドの互除法、拡張ユークリッドの互除法
最大公約数(GCD:Greatest Common Divisor)とは、2つの正整数aとbを素因数分解してときの共通項です。例えば
a=924(=2×2×3×7×11)
b=360(=2×2×2×3×3×5)
の最大公約数Gは
G=2×2×3=12
になります。
なお、最小公倍数(LCM:Least Common Multiple)は、aとbで割り切れる最小の整数のことです。
a=G×(7×11)
b=G×(2×3×5)
なので、最小公倍数Lは、
L=G×(7×11)×(2×3×5)
となります。すなわち、
L=a×b/G=924×360/12=27720
で計算できます。
ユークリッドの互除法は、最大公約数を計算する効率的なアルゴリズムです。
ユークリッドの互除法は、再帰関数を用いると簡潔な表現になります。
上のプログラムには、いくつかの改善点があります。
これらを考慮したプログラムを示します。
function euclid(n1, n2) {
var a, b, r; /* これらの変数は関数内でのみ有効。他のプログラムには無影響 */
a = Math.abs(n1); /* n1の絶対値 */
b = Math.abs(n2);
if (b > a) { /* b>a のときの置き換え */
r = a; a = b; b = r;
}
r = a % b;
while (r > 0) {
if (r > b-r) r = b-r; /* 剰余の修正 */
a = b;
b = r;
r = a % b;
}
return b;
}
aをbで割ったときの商をq、剰余をrとすると、
r=a-q×b
となります(a=924、b=360とすると、r=924-2×360=204)。
ここで、aとbの最大公約数 GCD(a,b) をgとすれば、aもbもgで割り切れるので、rもgで割り切れます(この段階では、gが24であることはわかっていないのですが、r=204(=17×24)はg=24で割り切れます)。すなわち、
GCD(a,b)=GCD(b,r) GCD(924,360)=GCD(360,204)
の関係があります。GCD(a,b) を求めるかわりに、GCD(b,r) を求めればよいし、しかも、r<b<aなので、GCD(b,r) を求めるほうが簡単です。
「b←a、r←a/bの剰余」の置き換えを繰り返すことにより、より小さい数の組み合わせになります。
そして、rが0になるまで繰り返されるのですから、最終的には、
GCD(a,b)=GCD(b,0)
となります(GCD(360,204)=GCD(204,156)=GCD(156,48)=GCD(48,12)=GCD(12,0) )。
ここで、GCD(b,0) とは、bと0との最大公約数であり、0はどんな数でも割り切れるのですから、GCD(b,0) =b であることは明白です。すなわち、
GCD(a,b)=b GCD(924,360)=12
となります。
最大公約数の計算量では、剰余の計算回数が尺度になります。私は理解していないのですが、ラメ(Lame')の定理によれば、最悪の場合でも、剰余の計算回数bの桁数の5倍以下であることが証明されています(b=360は3桁、3×5=15回以下)。
また、k回の剰余計算が必要となるとき、bはk番目のフィボナッチ数以上であることが証明されています(フィボナッチ数とは、F(k)=F(k-1)+F(k-2)、 F(1)=1、F(2)=1で定義される級数のこと。F(13)=233,F(14)=377)。
拡張ユークリッドの互除法は、自然数aとbを与えたとき、不定方程式
ax+by=g (gはaとbの最大公約数)
となる整数(負の場合もある)xとyを求める方法です。
例えば、a=924,b=360とすると、g=12となり、
924×(-7+30t)+360×(18-77t)=12
になります。このような、
x=-7+30t
y=18-77t (tは任意の整数)
を求める方法です。
ax+by=gを
a×xi+b×yi=ri ①
と表記します。
初期値を
x0=1,y0=0、r0=a(=924)
x1=0,y1=1、r1=b(=360)
とすれば、①は
i=0のとき、a×x0+b×y0=r0 → a×1+b×0=a
i=1のとき、a×x1+b×y1=r1 → a×0+b×1=b
となり、成立しています。
qiをri-1/riの整数部分、ri+1をその剰余とし、
ri+1=ri-1-qi×ri ②
xi+1=xi-1-qi×xi ③
yi+1=yi-1-qi×yi ④
と定義します。
i=1のときは、q1=r0/r1の整数部分なので、q1=924/360=2となります。
②③④から、
r2=r0-q1×r1 (924-2×360=204)
x2=x0-q1×x1 (1-2×0=1)
y2=y0-q1×y1 (0-2×1=-2)
となるので、これを①に代入すれば、
a×x2+b×y2=924×1-360×(-2)=204=r2
となり、①が成立しています。
i=kのときは
xk+1=xk-1-qk×xk
yk+1=yk-1-qk×yk
とすれば、
a×xk+1+b×yk+1=(a×xk-1+b×yk-1)+qk×(a×xk+b×yk)
となり、
a×xk-1+b×yk-1=rk-1
a×xk+b×yk=rk
ですから、
a×xk+1+b×yk+1=rk+1
となり、数学的帰納法により、①が常に成立することが証明されます。
ここで、ri-1=0になるときのriがaとbの最大公約数であることは、先のユークリッドの互除法での説明からわかります。
そして、①の関係を保ちながら式の変形をしてきたので、xiとyiは、
ax+by=g
を満足しています。
後述のトレースにより、a=924、b=360のときは、x=-7,y=18になり、
924×(-7)+360×(18)=12
となります。
不定方程式では、xとyの値はいくらでもあります。証明は省略しますが、上の計算で得られたxとyをx*、yをy*とし、最大公約数をgとしたとき、tを任意の整数として
x=x*+(b/g)×t (x=-7+30t)
y=y*-(a/g)×t (x=18-77t)
も、ax+by=gの解になります。
924×(-7+30t)+360×(18-77t)
=-6468+27720t+6480-27720t
=12
function ex_euclid(a, b, gcd, x, y){
r0 = a; r1 = b;
x0 = 1; x1 = 0;
y0 = 0; y1 = 1;
while (r1>0) {
q = Math.floor(r0/r1);
r2 = r0 % r1;
x2 = x0 - q*x1;
y2 = y0 - q*y1;
r0 = r1; r1 = r2;
x0 = x1; x1 = x2;
y0 = y1; y1 = y2;
}
gcd = r0;
x = x0;
y = y0;
}
q r0 r1 r2 x0 x1 x2 y0 y1 y2
2 924 360 204 1 0 1 0 1 -2
1 360 204 156 0 1 -1 1 -2 3
1 204 156 48 1 -1 2 -2 3 -5
3 156 48 12 -1 2 -7 3 -5 18
4 48 12 0 2 -7 30 -5 18-77
4 12 0 0 -7 30 30 18-77-77
↓ ↓ ↓ ↓ ↓ ↓
g 終了 x b/g y -a/g
基本解
a x b y GCD
924×(-7)+360×(18)=12
一般解
b/g=x1= 30 x=-7+30t
-a/g=y1=-77 y=18-77t