変数の個数が式の本数よりも多い方程式を「不定方程式」といいます。一次不定方程式 ax + by = c (c != 0) は、a, b, c が整数で、c が a と b の最大公約数の倍数であるとき、整数解を持つことが証明されています。これを「ベズーの等式」といいます。ベズーの等式は「拡張ユークリッドの互除法」を使って簡単に解くことができます。今回は拡張ユークリッドの互除法を簡単に説明します。
拡張といっても原理はユークリッドの互除法とまったく同じです。gcd(a, b) を求めるときに一組の解 (x, y) もいっしょに見つける、というものです。ポイントは数を x = xa * a + xb * b の形で表すところです。次の図を見てください。
x = xa * a + xb * b, y = ya * a + yb * b とおき、
x, y に対してユークリッドの互除法を適用する
(初期値は xa = 1, xb = 0, ya = 0, yb = 1)
x と y の商を q, 余りを z とすると
z = x - q * y = xa * a + xb * b - q * (ya * a + yb * b)
= (xa - q * ya) * a + (xb - q * yb) * b
= za * a + zb * b
となる
次は y と z で同様の処理を繰り返す
具体的には変数の値を
x, xa, xb = y, ya, yb
y, ya, yb = z, za, za
のように更新して処理を繰り返す
y が 0 になったとき、x が a と b の最大公約数で (xa, xb) が一つの解となる
具体的な例を示しましょう。4321 * x + 1234 * y = gcd(4321, 1234) を解いてみます。
(x, xa, xb) (y, ya, yb) q (z, za, zb) ------------------------------------------------------ (4321, 1, 0) (1234, 0, 1) 3 (619, 1, -3) (1234, 0, 1) (619, 1, -3) 1 (615, -1, 4) (619, 1, -3) (615, -1, 4) 1 (4, 2, -7) (615, -1, 4) (4, 2, -7) 153 (3, -307, 1075) (4, 2, -7) (3, -307, 1075) 1 (1, 309, -1082) (3, -307, 1075) (1, 309, -1082) 3 (0, -1234, 4321) (1, 309, -1082) (0, -1234, 4321) 4321 * 309 + 1234 * (-1082) = 1 になる
一般解は次のように求めることができます。
4321 * x + 1234 * y = 1 - 4321 * 309 + 1234 * (-1082) = 1 ----------------------------------- 4321 * (x - 309) + 1234 * (y - (-1082)) = 0 となるから 4321 * (x - 309) = -1234 * (y + 1082) が成り立つ -1234 は x - 309 の約数で 4321 は y + 1082 の約数だから、t を整数とすると、 x - 309 = -1234 * t => x = 309 - 1234 * t y + 1082 = 4321 * t => y = -1082 + 4321 * t となる
それでは、実際に値を代入して (x, y) の組を求めてみましょう。
t x y z = 4321 * x + 1234 * y ----------------------------------------- -4 5245 -18366 1 -3 4011 -14045 1 -2 2777 -9724 1 -1 1543 -5403 1 0 309 -1082 1 1 -925 3239 1 2 -2159 7560 1 3 -3393 11881 1 4 -4627 16202 1
プログラムも簡単です。Python でプログラムを作ると次のようになります。
リスト : 拡張ユークリッドの互除法
def ext_euclid(a, b):
xs = a, 1, 0
ys = b, 0, 1
while ys[0] != 0:
q, z = divmod(xs[0], ys[0])
xs, ys = ys, (z, xs[1] - q * ys[1], xs[2] - q * ys[2])
return xs
x, xa, xb をタプル xs に、y, ya, yb をタプル ys に、z, za, zb をタプルにまとめています。関数 divmod(x, y) は x と y の商と余りを求める関数です。あとはアルゴリズムをそのままプログラムしただけなので、難しいところはないと思います。それでは実際に試してみましょう。
>>> ext_euclid(8, 5) (1, 2, -3) >>> ext_euclid(9, 5) (1, -1, 2) >>> ext_euclid(11, 3) (1, -1, 4) >>> ext_euclid(4321, 1234) (1, 309, -1082) >>> ext_euclid(12357, 100102) (1, 30127, -3719)
次は「鶴亀算」を解いてみましょう。
鶏、犬、タコを変数 x, y, z とすると、この問題は次の連立方程式で表すことができます。
x + y + z = 24 2x + 4y + 8z = 102 z = 24 - x - y なので 2x + 4y + 8(24 - x - y) = 102 => - 6x - 4y + 192 = 102 => 6x + 4y = 90 => 3x + 2y = 45
この問題は不定方程式 3x + 2y = 45 を解けばよいことになります。2 と 3 の最大公約数は 1 なので、この式には整数解があることがわかります。このほかに x, y, z が 1 以上であることが条件になります。
ext_euclid(3, 2) を計算すると次のようになります。
>>> ext_euclid(3, 2) (1, 1, -1)
解の一つは (1, -1) を 45 倍した (45, -45) になります。一般解は次のようになります。
3 * x + 2 * y = 45 - 3 * 45 + 2 * (-45) = 45 ----------------------------------- 3 * (x - 45) + 2 * (y + 45) = 0 となるから 3 * (x - 45) = -2 * (y + 45) が成り立つ -2 は x - 45 の約数で 3 は y + 45 の約数だから、t を整数とすると、 x - 45 = -2 * t => x = 45 - 2 * t y + 45 = 3 * t => y = -45 + 3 * t となる
x >= 1, y >= 1 より t の範囲は 16 <= t <= 22 となります。実際に計算すると次のようになります。
>>> for t in range(16, 23):
... x = 45 - 2 * t
... y = -45 + 3 * t
... z = 24 - x - y
... print('x = {:d}, y = {:d}, z = {:d}'.format(x, y, z))
...
x = 13, y = 3, z = 8
x = 11, y = 6, z = 7
x = 9, y = 9, z = 6
x = 7, y = 12, z = 5
x = 5, y = 15, z = 4
x = 3, y = 18, z = 3
x = 1, y = 21, z = 2
鶴亀算は数理最適化のソルバーや制約論理プログラミングでも簡単に解くことができます。興味のある方は拙作のページ PuLP による数理最適化超入門 や 制約論理プログラミング超入門 をお読みください。
拙作のページ Puzzle DE Programming に 水差し問題 がありますが、それを解くのに拡張ユークリッドの互除法を使うことができます。
たとえば、8 リットルの容器 A と 5 リットルの容器 B を使う場合、8 * x + 5 * y = gcd(8, 5) を解くと、8 * 2 - 5 * 3 = 1 という解を見つけることができます。これは A で 2 回汲み B で 3 回捨てると、1 リットルの水が残ることを表します。4 リットルの水を汲みだす場合、両辺を 4 倍して 8 * 8 - 5 * 12 = 4 になります。A で 8 回汲み B で 12 回捨ててもいいのですが、一般解を求めるともっと少ない回数の解を見つけることができます。
8 * x + 5 * y = 4 - 8 * 8 - 5 * 12 = 4 --------------------------------- 8 * (x - 8) + 5 * (y + 12) = 0 x - 8 = -5 * t => x = 8 - 5 * t y + 12 = 8 * t => y = -12 + 8 * t t = 0, x = 8, y = -12 t = 1, x = 3, y = -4 t = 2, x = -2, y = 4
上記のように、B で 4 回汲み A で 2 回捨てる解を見つけることができます。
中国の剰余定理 (Chinese remainder theorem) は整数の剰余に関する定理です。具体的には、連立合同式に解があることを示す定理です。
プログラムは ext_euclid() を使うと簡単です。
リスト : 中国の剰余定理
def crt(a, m, b, n):
g, x, y = ext_euclid(m, n)
return (a * n * y + b * m * x) % (m * n) if g == 1 else None
関数 crt() は m で割った余りが a で、n で割った余りが b となる整数を求めます。ext_euclid() で m * x + n *y = 1 を満たす (x, y) を求めます。あとは、a * n * y + b * m * x を計算し、m * n との剰余を求めて返すだけです。最大公約数 g が 1 でなければ、m と n は互いに素ではないので None を返します。
それでは実際に試してみましょう。3 で割ると 2 余り、5 で割ると 4 余る数を求めます。
>>> crt(2, 3, 4, 5) 14 >>> 14 % 3 2 >>> 14 % 5 4
解は 14 になります。確かに 14 % 3 は 2 に、14 % 5 は 4 になります。解を合同式で表すと x ≡ 14 (mod 15) になります。
中国の剰余定理は合同式の本数を増やしても成立します。
これをそのままプログラムすると、次のようになります。
リスト : 中国剰余定理 (一般化)
def crtn(xs):
a, m = xs[0]
for b, n in xs[1:]:
z = crt(a, m, b, n)
if z is None: return None
a = z
m *= n
return a
引数 xs はリストで、要素は合同式を表すタプル (ai, ni) です。最初に先頭要素を取り出して、変数 a, m にセットします。あとは for ループで、xs の要素を順番に取り出して変数 b, n にセットし、crt() で解 z を求めます。z が None であれば None を返します。そして、a の値を z に更新し、m に n を乗算します。最後に a の値を返します。これで連立合同式の解を求めることができます。
それでは実際に試してみましょう。
>>> crtn([(2, 3), (3, 5), (2, 7)]) 23 >>> 23 % 3 2 >>> 23 % 5 3 >>> 23 % 7 2
>>> crtn([(2, 3), (3, 5), (2, 7), (4, 11)]) 653 >>> 653 % 3 2 >>> 653 % 5 3 >>> 653 % 7 2 >>> 653 % 11 4
次は、中国剰余定理を整数 ni (i = 1, 2, ..., k) が互いに素でない場合に拡張します。
\(lcm(m, n)\) は m と n の最小公倍数を表します。m と n が互いに素であれば \(\gcd(m, n) = 1\) になるので、必要十分条件を満たします。
もちろん、合同式の本数を増やしても成立します。
証明は互いに素の場合と同様 (数学的帰納法) にできるので省略します。
プログラムは次のようになります。
リスト : 中国剰余定理 (互いに素数でない場合)
def crt1(a, m, b, n):
d, x, y = ext_euclid(m, n)
if (b - a) % d != 0: return None
s = (b - a) // d
return (s * m * x + a) % math.lcm(m, n)
def crtn1(xs):
a, m = xs[0]
for b, n in xs[1:]:
z = crt1(a, m, b, n)
if z is None: return None
a = z
m = math.lcm(m, n)
return a
math.lcm() は最小公倍数を求める Python の関数です。crt1() の場合、b - a が最大公約数 d で割り切れないとき、解が無いので None を返します。あとは、どちらの関数も lcm() を使って最小公倍数を求めていることに注意してください。
それでは実際に試してみましょう。
>>> crt(2, 3, 4, 5) 14 >>> crt1(2, 3, 4, 5) 14 >>> crt(20, 30, 40, 50) >>> crt1(20, 30, 40, 50) 140 >>> 140 % 30 20 >>> 140 % 50 40 >>> crtn([(2, 3), (3, 5), (2, 7)]) 23 >>> crtn1([(2, 3), (3, 5), (2, 7)]) 23 >>> crtn([(20, 30), (30, 50), (20, 70)]) >>> crtn1([(20, 30), (30, 50), (20, 70)]) 230 >>> 230 % 30 20 >>> 230 % 50 30 >>> 230 % 70 20 >>> crtn1([(2, 3), (3, 5), (2, 7), (4, 11)]) 653 >>> crtn1([(20, 30), (30, 50), (20, 70), (40, 110)]) 6530 >>> 6530 % 30 20 >>> 6530 % 50 30 >>> 6530 % 70 20 >>> 6530 % 110 40
正常に動作しているようです。興味のある方はいろいろ試してみてください。
数学の世界では、加減乗除ができる数の集合のことを「体」といいます。特に、数が n 個しかないものを「有限体 (finite field)」または「ガロア体 (Galois field)」といいます。ここでは GF(n) と表記することにします。n が素数のとき、その剰余の集合 {0, 1, ..., n - 1} は有限体になります。簡単な例として GF(2) と GF(3) の加法と乗法を示します。
GF(2) + | 0 1 * | 0 1 ---+------- ---+------- 0 | 0 1 0 | 0 0 1 | 1 0 1 | 0 1 GF(3) + | 0 1 2 * | 0 1 2 ---+---------- ---+---------- 0 | 0 1 2 0 | 0 0 0 1 | 1 2 0 1 | 0 1 2 2 | 2 0 1 2 | 0 2 1
mod n の世界では、加法と乗法は簡単に計算できます。加算または乗算したあと mod n を計算するだけです。減算は負数 (加法の逆元) を、徐算は逆数 (乗法の逆元) を考えます。x + y = 0 を満たす y が x の負数になるので、y = n - x であることがわかります。0 以外の数の逆数が存在すれば徐算が可能になります。乗法の表を見れば、GF(2) は 1 * 1 = 1, GF(3) は 1 * 1 = 1, 2 * 2 = 1 となるので、徐算が成り立つことがわかります。
n が素数ではない場合、逆数が存在しない数があります。たとえば、mod 4 の加法と乗法は次のようになります。
mod 4 の加法と乗法 + | 0 1 2 3 * | 0 1 2 3 --+------------- ---+------------- 0 | 0 1 2 3 0 | 0 0 0 0 1 | 1 2 3 0 1 | 0 1 2 3 2 | 2 3 0 1 2 | 0 2 0 2 3 | 3 0 1 2 3 | 0 3 2 1
mod 4 の場合、加算、乗算、減算は成り立ちますが、徐算で 2 の逆数が存在しないことが乗法の表からわかります。つまり、有限体にはならないわけです。
4 の剰余 (0, 1, 2, 3) は GF(4) になりませんが、GF(2) に新しい数を付け加えることで、数が 4 つある有限体を作ることができます。新しい数を付加して体を拡張することを「体の拡大」といい、拡大された体を「拡大体」といいます。たとえば、実数の世界では代数方程式 x2 + 1 = 0 に解はありませんが、虚数 i を導入した複素数の世界では解くことができますね。複素数の世界は実数の拡大体になるわけです。
GF(2) を拡大する場合、GF(2) 上の代数方程式を考えます。2 次方程式を列挙すると次のようになります。
この中で解を持たないのは式 4 です。『x の関数 f(x) において、f(a) = 0 ならば (x - a) で割り切れる』という性質を「因数定理」といいます。式 1, 2, 3 を因数分解すると次のようになります。
式 4 は GF(2) で因数分解することはできません。これを「既約多項式」といいます。既約多項式は多項式の素数みたいなもので、1 次多項式 x, x + 1 で割り切れることはありません。この式の解の一つを a としましょう。a2 + a + 1 = 0 が成り立つので、a2 = a + 1 であることがわかります。次に、a の累乗を計算します。
a1 = a a2 = a + 1 a3 = a * a2 = a * (a + 1) = a + 1 + a = 1 a4 = a1+3 = a a5 = a2+3 = a + 1 a6 = a3+3 = 1
a, a + 1, 1 が順番に現れています。これを多項式で表すと 1, x, x + 1 になり、1 次多項式がすべて現れていることがわかります。このように、累乗を計算すると n 次未満の多項式がすべて現れる既約多項式のことを「原始多項式」といいます。これに 0 を加えた 4 つの要素 {0, 1, x, x + 1} を数として加法と乗法を考えてみましょう。基本的には式を計算して、それを既約多項式 x2 + x + 1 で割った剰余が値になります。下図を見てください。
+ | 0 1 x x+1 * | 0 1 x x+1 ----+----------------- ----+------------------ 0 | 0 1 x x+1 0 | 0 0 0 0 1 | 1 0 x+1 x 1 | 0 1 x x+1 x | x x+1 0 1 x | 0 x x+1 1 x+1 | x+1 x 1 0 x+1 | 0 x+1 1 x
加法の場合は特に簡単で、係数ごとに加算して mod 2 をとる、つまり係数ごとの排他的論理和になります。乗算の場合は累乗の形式に変換して計算すると簡単です。つまり、am * an = an+m = a(n+m) mod 3 になります。a2 = a + 1 の関係を使って式を簡約してもかまいません。
減法の場合、加法表より負数が自分自身であることがわかるので、加法と同様に係数ごとの排他的論理和になります。徐算の場合、乗法表より 1, x, x + 1 の逆数がすべて存在するので、徐算が成り立つことがわかります。つまり、{0, 1, x, x + 1} は有限体になるわけです。多項式の係数をビットで表すと、加法と乗法は次のようになります。
+ | 00 01 10 11 * | 00 01 10 11 ----+----------------- ----+------------------ 00 | 00 01 10 11 00 | 00 00 00 00 01 | 01 00 11 10 01 | 00 01 10 11 10 | 10 11 00 01 10 | 00 10 11 01 11 | 11 10 01 00 11 | 00 11 01 10
p を素数とし、n を自然数とすると、要素数が pn の有限体は必ず存在することが知られています。これを GF(pn) と表記します。GF(2n) はコンピュータとの相性が良いので、エラー検出、乱数の生成、暗号など多くの分野で応用されています。
ご参考までに、GF(23) と GF(32) の加法表と乗法表を示します。
GF(23) 原始多項式 f(x) = x3 + x + 1 f(a) = 0 となる要素 a を考える a1 = a (010) a2 = a2 (100) a3 = a + 1 (011) a4 = a2 + a (110) a5 = a2 + a + 1 (111) a6 = a2 + 1 (101) a7 = 1 (001)
GF(23) の加法表
+ | 000 001 010 011 100 101 110 111
----+---------------------------------
000 | 000 001 010 011 100 101 110 111
001 | 001 000 011 010 101 100 111 110
010 | 010 011 000 001 110 111 100 101
011 | 011 010 001 000 111 110 101 100
100 | 100 101 110 111 000 001 010 011
101 | 101 100 111 110 001 000 011 010
110 | 110 111 100 101 010 011 000 001
111 | 111 110 101 100 011 010 001 000
GF(23) の乗法表, (n) の n は指数
(7) (1) (3) (2) (6) (4) (5)
* | 000 001 010 011 100 101 110 111
--------+----------------------------------
000 | 000 000 000 000 000 000 000 000
001 (7) | 000 001 010 011 100 101 110 111
010 (1) | 000 010 100 110 011 001 111 101
011 (3) | 000 011 110 101 111 100 001 010
100 (2) | 000 100 011 111 110 010 101 001
101 (6) | 000 101 001 100 010 111 011 110
110 (4) | 000 110 111 001 101 011 010 100
111 (5) | 000 111 101 010 001 110 100 011
GF(32) 原始多項式 f(x) = x2 + x + 2 f(a) = 0 となる要素 a を考える a1 = a (1 0) a2 = -a - 2 = 2a + 1 (2 1) a3 = a(2a + 1) = 4a + 2 + a = 2a + 2 (2 2) a4 = a(2a + 2) = 4a + 2 + 2a = 2 (0 2) a5 = 2a (2 0) a6 = 2a * a = 2(2a + 1) = a + 2 (1 2) a7 = a(a + 2) = 2a + 1 + 2a = a + 1 (1 1) a8 = a(a + 1) = 2a + 1 + a = 1 (0 1)
GF(32) の加法表
+ | (0 0) (0 1) (0 2) (1 0) (1 1) (1 2) (2 0) (2 1) (2 2)
------+-------------------------------------------------------
(0 0) | (0 0) (0 1) (0 2) (1 0) (1 1) (1 2) (2 0) (2 1) (2 2)
(0 1) | (0 1) (0 2) (0 0) (1 1) (1 2) (1 0) (2 1) (2 2) (2 0)
(0 2) | (0 2) (0 0) (0 1) (1 2) (1 0) (1 1) (2 2) (2 0) (2 1)
(1 0) | (1 0) (1 1) (1 2) (2 0) (2 1) (2 2) (0 0) (0 1) (0 2)
(1 1) | (1 1) (1 2) (1 0) (2 1) (2 2) (2 0) (0 1) (0 2) (0 0)
(1 2) | (1 2) (1 0) (1 1) (2 2) (2 0) (2 1) (0 2) (0 0) (0 1)
(2 0) | (2 0) (2 1) (2 2) (0 0) (0 1) (0 2) (1 0) (1 1) (1 2)
(2 1) | (2 1) (2 2) (2 0) (0 1) (0 2) (0 0) (1 1) (1 2) (1 0)
(2 2) | (2 2) (2 0) (2 1) (0 2) (0 0) (0 1) (1 2) (1 0) (1 1)
GF(32) の乗法表, (n) の n は指数
(8) (4) (1) (7) (6) (5) (2) (3)
* | (0 0) (0 1) (0 2) (1 0) (1 1) (1 2) (2 0) (2 1) (2 2)
----------+-------------------------------------------------------
(0 0) | (0 0) (0 0) (0 0) (0 0) (0 0) (0 0) (0 0) (0 0) (0 0)
(0 1) (8) | (0 0) (0 1) (0 2) (1 0) (1 1) (1 2) (2 0) (2 1) (2 2)
(0 2) (4) | (0 0) (0 2) (0 1) (2 0) (2 2) (2 1) (1 0) (1 2) (1 1)
(1 0) (1) | (0 0) (1 0) (2 0) (2 1) (0 1) (1 1) (1 2) (2 2) (0 2)
(1 1) (7) | (0 0) (1 1) (2 2) (0 1) (1 2) (2 0) (0 2) (1 0) (2 1)
(1 2) (6) | (0 0) (1 2) (2 1) (1 1) (2 0) (0 2) (2 2) (0 1) (1 0)
(2 0) (5) | (0 0) (2 0) (1 0) (1 2) (0 2) (2 2) (2 1) (1 1) (0 1)
(2 1) (2) | (0 0) (2 1) (1 2) (2 2) (1 0) (0 1) (1 1) (0 2) (2 0)
(2 2) (3) | (0 0) (2 2) (1 1) (0 2) (2 1) (1 0) (0 1) (2 0) (1 2)
拡張ユークリッドの互除法は、合同式の徐算で逆数を求めるときにも役に立ちます。\(a \times b \equiv 1 \pmod n\) において a の逆数 b を求める場合、\(a \times x + n \times y = 1\) を求めます。ここで、\(\gcd(a, n) = 1\) が条件になることに注意してください。両辺に mod n を施すと \(n \times y\) は 0 になるので、\(a \times x \equiv 1 \pmod n\) が成り立ちます。つまり、x が a の逆数になるわけです。簡単な例を示しましょう。
>>> for a in range(1, 11): ... print(ext_euclid(a, 11)) ... (1, 1, 0) (1, -5, 1) (1, 4, -1) (1, 3, -1) (1, -2, 1) (1, 2, -1) (1, -3, 2) (1, -4, 3) (1, 5, -4) (1, -1, 1)
n が素数の場合、n とその剰余は互いに素になります。合同式で x の負数 -x は (n - x) になるので、mod 11 における 1 から 10 までの逆数は次のようになります。
x 1/x ------------------------- 0 --- 1 1 1 mod 11 = 1 2 6 12 mod 11 = 1 3 4 12 mod 11 = 1 4 3 12 mod 11 = 1 5 9 45 mod 11 = 1 6 2 12 mod 11 = 1 7 8 56 mod 11 = 1 8 7 56 mod 11 = 1 9 5 45 mod 11 = 1 10 10 100 mod 11 = 1
n が素数でない場合、n とその剰余は互いに素になるとは限りません。次の例を見てください。
>>> for x in range(1, 8): print(ext_euclid(x, 8)) ... (1, 1, 0) (2, 1, 0) (1, 3, -1) (4, 1, 0) (1, -3, 2) (2, -1, 1) (1, -1, 1)
2, 4, 6 と 8 の最大公約数は 1 にならないので、その逆数を求めることはできません。つまり、mod 8 の世界では 2, 4, 6 で割り算することはできないわけです。
また、フェルマーの小定理から逆数を求めることもできます。
\(a \times a^{n-2} \equiv 1 \pmod n\) が成り立つので、逆数は \(a^{n-2}\) になります。それでは実際に試してみましょう。
>>> for x in range(1, 11): print(x, x**9, (x**9) % 11) ... 1 1 1 2 512 6 3 19683 4 4 262144 3 5 1953125 9 6 10077696 2 7 40353607 8 8 134217728 7 9 387420489 5 10 1000000000 10
拡張ユークリッドの互除法と同じ結果になりましたね。なお、簡単な例ということで累乗を計算してから mod をとりましたが、実用的なプログラムではやらないほうがよいでしょう。mod n の世界で累乗を計算する関数を作ったほうがよいと思います。興味のある方は挑戦してみてください。