一次不定方程式
著者:梅谷 武
語句:オイラー,アールヤバタ,不定方程式
一次不定方程式のユークリッドの互除法を使った解法について述べる。
作成:2006-04-25
更新:2013-06-16
 ガウスの『整数論』には、一次不定方程式の一般的な解法はオイラーEuler, Leonhard, 1707-1783によって最初に報告されたとありますが、その千年以上前にインドのアールヤバタĀrayabhata, 476-?は当時のインドの数学においてはクッタカ(粉砕法)と呼ばれていたユークリッドの互除法によって一次不定方程式を解いていました。彼の著書『アールヤバティーヤ』の詩節32-33に、

例題2.6.2

aで割ればr1余り、bで割ればr2余るような整数nを求めよ。
という問題があります。これを
n = ax + r1 = by + r2, 0≦r1<a, 0≦r2<b
と表し、r2 - r1 = cとおくと2元一次不定方程式
ax - by = c
となりますが、アールヤバタはこれを割り算でa,bを小さい数に還元して解いています。小さい数に還元するというのはa>b>0と仮定してa = bq + r, 0≦r<bなるq,rを求め、aに代入し
(bq + r)x - by
=
c
rx - b(y - qx)
=
c
rx - by'
=
c, y' = y - qx
というようにabより小さい数rに置き換えることです。アールヤバタは簡単に解けるぐらいに係数が小さくなった時点で打ち切っていましたが、これを繰り返すといつかはr=0となり、dx'=c, d=gcd(a,b)という一元一次不定方程式に還元され、d|cの場合に解けることがわかります。この節ではこのように一次不定方程式をユークリッドの互除法を使って解くことについて考えます。
 整数ai, i=1,⋯,nが生成するイデアル(a1, ⋯, an)
{ a1 x1 + ⋯ + an xn | xi∈ℤ }
という集合でした。この節ではこの集合の元の形式におけるxin個の変数として、ある整数cに対する
a1 x1 + ⋯ + an xn = c
という方程式の整数解について考えます。一般に代数方程式の整数解を求めるときに、その方程式を不定方程式ふていほうていしき, indeterminate equationといいます。上の式のように方程式が変数の一次式であるときに一次不定方程式といいます。
{ a1 x1 + ⋯ + an xn | xi∈ℤ } ∩ ℕの最小元をdとすると
(d) = { a1 x1 + ⋯ + an xn | xi∈ℤ }
となることがわかっていますから、この不定方程式の解が存在することとc∈(d)、すなわちcdの倍数であることは同じことです。まず、c=dのときを考えましょう。2変数の場合を具体例でみてみます。
15X + 12Y = 3, 3 = gcd(15,12)
とします。これはX-Y平面上の直線の方程式と考えることができます。この直線をパラメータ表示してみましょう。
X-1
-4
=
Y+1
5
= t
から
lb48
X
Y
rb48 = lb48
1
-1
rb48 + lb48
-4
5
rb48 t
となりますが、この整数解はちょうどtを整数としたときに得られることがわかります。これを一般化して命題の形にまとめておきましょう。

命題2.6.8

0でない整数a,bに関する一次不定方程式
aX + bY = d, d = gcd(a,b), a=a'd, b=b'd
は無限個の解をもち、その一般解は一つの特殊解を(x0,y0)としたときに
{ (x0 - b't, y0 + a't) | t∈ℤ }
と書くことができる。

証明

 一般解を(x,y)とするとa'(x-x0)=-b'(y-y0)となり(a',b')=1よりx-x0=-b't, y-y0=a't, t∈ℤと書ける。■
 上の命題における特殊解はユークリッドの互除法を少し拡張することで最大公約数と同時に計算することができます。互除法とは、0でない整数a,bが与えられたときに、
a0 = a, b0 = b, a0
=
q0 b0 + r0, 0<r0<|b0|
a1 = b0, b1 = r0, a1
=
q1 b1 + r1, 0<r1<|b1|
an = bn-1, bn = rn-1, an
=
qn bn
という数列(ai),(bi),(qi),(ri)rn=0となるまで計算するものでした。これを次のように書き直します。すなわち、0でない整数a,b(a>b)が与えられたときに、x0 = a, x1 = bとし、
               x0
=
k0 x1 + x2, 0<x2<|x1|
x1
=
k1 x2 + x3, 0<x3<|x2|
xi-1
=
ki-1 xi + xi+1, 0<xi+1<|xi|
xn
=
kn xn+1
という数列(xi),(ki)xn+2=0となるまで計算するものと考えます。ここで数列(xi)は狭義の単調減少数列になっているので必ず有限回で止まります。これを行列で表現してみましょう。xi,xi+1からxi+1,xi+2を計算する部分は次のように書くことができます。
lb48
xi+1
xi+2
rb48 = lb48
0
1
1
-ki
rb48 lb48
xi
xi+1
rb48
したがって
lb48
d
0
rb48 = lb48
0
1
1
-kn
rb48lb48
0
1
1
-k0
rb48 lb48
a
b
rb48
となりますから、行列
lb48
0
1
1
-kn
rb48lb48
0
1
1
-k0
rb48= lb48
un
wn
vn
zn
rb48
を計算していけばd=un a + wn bとなり、特殊解(un, wn)が得られます。しかし行列の計算を行うと計算量が多くなりますので、工夫して計算量を減らします。クヌースによれば以下のアイディアはGordon H. Bradleyによるものです。
wn =
d - un a
b
であることからdunだけを求めるようにしてみましょう。
lb48
un
wn
vn
zn
rb48 = lb48
0
1
1
-kn
rb48 lb48
un-1
wn-1
vn-1
zn-1
rb48
において、unun-1,vn-1にしか依存していないことに注目すると、数列(ui),(vi)だけを計算していけば最終的にunが得られることがわかります。これを算法として整理し、そのままプログラミングしてみます。

算法2.6.11 拡張互除法

0でない二つの整数a,bの最大公約数dax+by=dなるx,yを求める。
Step 1. s = a; t = b; u = 1; v = 0
Step 2. s = tq + r, 0≦r<|b|となるq,rを求める。
Step 3. もしr = 0であればStep.5へ。
Step 4. s = t; t = r; temp = v; u = u - q*v; u = temp
Step 5. d = s; x = u; y = ( s - u * a ) / b

sample2.6.3.py

def gcd_ext( a, b, val = [ 0, 0 ] ):
  if ( a == 0 ) or ( b == 0 ):
    return 0
  s = a; t = b; u = 1; v = 0
  while t != 0:
    q = s / t; r = s % t
    s = t; t = r
    temp = v; v = u - q * v; u = temp
  val[0] = u
  val[1] = ( s - u * a ) / b
  if s < 0:
    val[0] = -val[0]
    val[1] = -val[1]
    return -s
  else:
    return s
 
a = 9749560641517; b = 8770036831691
val = [ 0, 0 ]
d = gcd_ext( a, b, val )
x = val[0]; y = val[1]
print 'gcd( %d, %d ) = %d' % ( a, b, d )
print 'a*(%d) + b*(%d) = %d' % ( x, y, (a*x + b*y) )

sample2.6.3.pyの実行結果

gcd( 9749560641517, 8770036831691 ) = 6972593
a*(251903) + b*(-280038) = 6972593

sample2.6.3.rb

def gcd_ext(a, b, val = [0, 0])
  if ( a == 0 ) or ( b == 0 )
    return 0
  end
  s = a; t = b; u = 1; v = 0
  while t != 0
    q = s / t; r = s % t
    s = t; t = r
    temp = v; v = u - q * v; u = temp
  end
  val[0] = u
  val[1] = (s - u * a) / b
  if s < 0
    val[0] = -val[0]
    val[1] = -val[1]
    return -s
  else
    return s
  end
end
 
a = 9749560641517; b = 8770036831691; val = [0, 0]
d = gcd_ext(a, b, val)
x = val[0]; y = val[1]
printf "gcd(%d, %d) = %d\n", a, b, d
printf "a*(%d) + b*(%d) = %d\n", x, y, (a*x + b*y)

sample2.6.3.rbの実行結果

gcd(9749560641517, 8770036831691) = 6972593
a*(251903) + b*(-280038) = 6972593
 ここまでで2元一次不定方程式aX+bY=d, d=gcd(a,b)を解くことができました。次に任意の整数cに対するaX+bY=cの解法について考えます。cdの倍数でなければ解は存在しませんから、c=c'dと仮定します。図2.6.6の例で考えましょう。
15X + 12Y = 2・3, 3 = gcd(15,12)
とすると、この直線のパラメータ表示は元のパラメータ表示で始点ベクトルをちょうど2倍することで得られます。
lb48
X
Y
rb48 = lb48
2
-2
rb48 + lb48
-4
5
rb48 t
この整数解はtを整数としたときに得られます。これを一般化して命題の形にまとめます。

命題2.6.17

0でない整数a,bに関する一次不定方程式
aX + bY = c, d = gcd(a,b), a=a'd, b=b'd, c=c'd
は無限個の解をもち、その一般解はaX + bY = dの一つの特殊解を(x0,y0)としたときに
{ (c'x0 - b't, c'y0 + a't) | t∈ℤ }
と書くことができる。

証明

 演習とする。■
 さて、一般のn元一次不定方程式a1x1 + ⋯ + anxn =cの解法ですが、線形代数の立場から見ると方程式がn次元ユークリッド空間nの超平面の方程式<`a,`x>=cの形をしていることから、この超平面上の整数の格子点を求めることと同じです。ここではn=3の場合だけを考えることにします。
a1 X1 + a2 X2 + a3 X3 = d, d = gcd(a1,a2,a3)
とし、この特殊解(x1,x2,x3)が与えられたとしましょう。a1(X1-x1)+a2(X2-x2)+a3(X3-x3)=0から
a3 ( X3 - x3 ) = - ( a1 ( X1 - x1 ) + a2 ( X2 - x2 ) )
ここでa3 = a3'd, a1 = a1'd'd, a2 = a2'd'd, d'=gcd(a1,a2)とすると、
a3' ( X3 - x3 ) = -d'( a1' ( X1 - x1 ) + a2' ( X2 - x2 ) )
(a3',d')=1よりd'|( X3 - x3 )ですから、
X3 - x3 = d'u, u ∈ ℤ
と書くことができます。これを代入して、
a3' u = - ( a1' ( X1 - x1 ) + a2' ( X2 - x2 ) )
X1' = X1 - x1, X2' = X2 - x2とおくと
a1' X1' + a2' X2' = - a3'u, u ∈ ℤ
となります。この一般解はa1' X1' + a2' X2' = 1の一般解に- a3'uを掛けたものになりますのでこの特殊解を(s1,s2)とすると、
lc48
X1'
=
-a3's1u - a2't
X2'
=
-a3's2u + a1't
, t,u ∈ ℤ
ここまでの結果を整理すると求める一般解
lc96
X1
=
x1 - a3's1u - a2't
X2
=
x2 - a3's2u + a1't
X3
=
x3 + d'u
, t,u ∈ ℤ
が得られます。この一般解が実際に求める方程式の解になっていること、つまり整数で3次元空間内の平面a1 X + a2 Y + a3 Z = d上にあることは容易に確かめることができます。これを命題の形にまとめます。

命題2.6.20

0でない整数a,b,cに関する一次不定方程式
a1 X1 + a2 X2 + a3 X3 = d, d = gcd(a1,a2,a3)
は無限個の解をもち、その一般解は特殊解を(x1,x2,x3)とし、
a3 = a3'd, a1 = a1'd'd, a2 = a2'd'd, d'=gcd(a1,a2)
とおいたときのa1' X + a2' Y = 1の特殊解を(s1,s2)としたとき、次のように書くことができる。
lc96
X1
=
x1 - a3's1u - a2't
X2
=
x2 - a3's2u + a1't
X3
=
x3 + d'u
, t,u ∈ ℤ
 この命題から3元一次不定方程式の解法はその特殊解を求めることに帰着されます。これはイデアルとして(a1, a2, a3) = ((a1, a2), a3)となることからd'=gcd(a1, a2)とおくとd = gcd(a1, a2, a3) = gcd(d', a3)となり、
a1 X + a2 Y = d'
の特殊解を(x0, y0)とし、
d' U + a3 Z = d
の特殊解を(u0, z0)とすれば、(x0 u0, y0 u0, z0)が求める特殊解となり、拡張互除法を2回適用することで計算できます。この方法をプログラミングしてみましょう。

sample2.6.5.py

def gcd_ext( a, b, val = [ 0, 0 ] ):
  if ( a == 0 ) or ( b == 0 ):
    return 0
  s = a; t = b; u = 1; v = 0
  while t != 0:
    q = s / t; r = s % t
    s = t; t = r
    temp = v; v = u - q * v; u = temp
  val[0] = u
  val[1] = ( s - u * a ) / b
  if s < 0:
    val[0] = -val[0]
    val[1] = -val[1]
    return -s
  else:
    return s
 
def gcd_ext3( a1, a2, a3, val = [ 0, 0, 0 ] ):
  val1 = [ 0, 0 ]; val2 = [ 0, 0 ]
  d1 = gcd_ext( a1, a2, val1 )
  d = gcd_ext( d1, a3, val2 )
  val[0] = val1[0]*val2[0]
  val[1] = val1[1]*val2[0]
  val[2] = val2[1]
  return d
 
a = 2418976464533; b = 2587070861497
c = 69314547013;
val = [ 0, 0, 0 ]
d = gcd_ext3( a, b, c, val )
x = val[0]; y = val[1]; z = val[2]
print 'gcd( %d, %d, %d ) = %d' % ( a, b, c, d )
print 'a*(%d) + b*(%d) + c*(%d) = %d' % ( x, y, z, (a*x + b*y + c*z) )

sample2.6.5.pyの実行結果

gcd( 2418976464533, 2587070861497, 69314547013 ) = 9941
a*(4809212484) + b*(-4496734855) + c*(-1392) = 9941

sample2.6.5.rb

def gcd_ext(a, b, val = [0, 0])
  if ( a == 0 ) or ( b == 0 )
    return 0
  end
  s = a; t = b; u = 1; v = 0
  while t != 0
    q = s / t; r = s % t
    s = t; t = r
    temp = v; v = u - q * v; u = temp
  end
  val[0] = u
  val[1] = (s - u * a) / b
  if s < 0
    val[0] = -val[0]
    val[1] = -val[1]
    return -s
  else
    return s
  end
end
 
def gcd_ext3(a1, a2, a3, val = [0, 0, 0])
  val1 = [0, 0]; val2 = [0, 0]
  d1 = gcd_ext(a1, a2, val1)
  d  = gcd_ext(d1, a3, val2)
  val[0] = val1[0]*val2[0]
  val[1] = val1[1]*val2[0]
  val[2] = val2[1]
  return d
end
 
a = 2418976464533; b = 2587070861497; c = 69314547013
val = [0, 0, 0]
d = gcd_ext3(a, b, c, val)
x = val[0]; y = val[1]; z = val[2]
printf "gcd(%d, %d, %d) = %d\n", a, b, c, d
printf "a*(%d) + b*(%d) + c*(%d) = %d\n",
       x, y, z, (a*x + b*y + c*z)

sample2.6.5.rbの実行結果

gcd(2418976464533, 2587070861497, 69314547013) = 9941
a*(4809212484) + b*(-4496734855) + c*(-1392) = 9941
 一般のn元一次不定方程式の特殊解は同じように拡張ユークリッド算法をn-1回適用することで求めることができます。一般解の構造については連立一次不定方程式の解法とともに線形代数を使って説明するのが適切ですので、ここでは詳細を述べないことにします。

参考文献

[1] 林 隆夫, インドの数学―ゼロの発明(中公新書), 中央公論社, 1993
[2] 木田 祐司, 牧野 潔夫, UBASICによるコンピュータ整数論, 日本評論社, 1994
[3] D.E.Knuth(中川圭介訳), 準数値算法―算術演算, サイエンス社, 1986
img
 
 
[4] D.E.Knuth(有沢誠他訳), The Art of Computer Programming (2) 日本語版 Seminumerical algorithms, アスキー, 2004
img
 
 
人  物
オイラー Euler, Leonhard, 1707-1783
スイスの数学者。ライプニッツの形式的代数計算を継承、発展させ数学史上最大級の業績を残す。その全集の刊行は21世紀となった現代においてもまだ完結していない。
アールヤバタ Ārayabhata, 476-?
インドの天文学者、数学者。
 
数  学
不定方程式 ふていほうていしき, indeterminate equation
 
Published by SANENSYA Co.,Ltd.