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

という問題があります。これを と表し、r₂ － r₁ = cとおくと２元一次不定方程式 となりますが、アールヤバタはこれを割り算でa,bを小さい数に還元して解いています。小さい数に還元するというのはa＞b＞0と仮定してa = bq + r, 0≦r＜bなるq,rを求め、aに代入し というようにaをbより小さい数rに置き換えることです。アールヤバタは簡単に解けるぐらいに係数が小さくなった時点で打ち切っていましたが、これを繰り返すといつかはr=0となり、dx'=c, d=gcd(a,b)という一元一次不定方程式に還元され、d|cの場合に解けることがわかります。この節ではこのように一次不定方程式をユークリッドの互除法を使って解くことについて考えます。

 整数a_i, i=1,⋯,nが生成するイデアル(a₁, ⋯, a_n)は という集合でした。この節ではこの集合の元の形式におけるx_iをn個の変数として、ある整数cに対する という方程式の整数解について考えます。一般に代数方程式の整数解を求めるときに、その方程式を不定方程式ふていほうていしき, indeterminate equationといいます。上の式のように方程式が変数の一次式であるときに一次不定方程式といいます。

 { a₁ x₁ + ⋯ + a_n x_n | x_i∈ℤ } ∩ ℕの最小元をdとすると となることがわかっていますから、この不定方程式の解が存在することとc∈(d)、すなわちcがdの倍数であることは同じことです。まず、c=dのときを考えましょう。２変数の場合を具体例でみてみます。 とします。これはX-Y平面上の直線の方程式と考えることができます。この直線をパラメータ表示してみましょう。

から となりますが、この整数解はちょうどtを整数としたときに得られることがわかります。これを一般化して命題の形にまとめておきましょう。

 上の命題における特殊解はユークリッドの互除法を少し拡張することで最大公約数と同時に計算することができます。互除法とは、0でない整数a,bが与えられたときに、 という数列(a_i),(b_i),(q_i),(r_i)をr_n=0となるまで計算するものでした。これを次のように書き直します。すなわち、0でない整数a,b(a＞b)が与えられたときに、x₀ = a, x₁ = bとし、 という数列(x_i),(k_i)をx_n+2=0となるまで計算するものと考えます。ここで数列(x_i)は狭義の単調減少数列になっているので必ず有限回で止まります。これを行列で表現してみましょう。x_i,x_i+1からx_i+1,x_i+2を計算する部分は次のように書くことができます。 したがって となりますから、行列 を計算していけばd=u_n a + w_n bとなり、特殊解(u_n, w_n)が得られます。しかし行列の計算を行うと計算量が多くなりますので、工夫して計算量を減らします。クヌースによれば以下のアイディアはGordon H. Bradleyによるものです。 であることからdとu_nだけを求めるようにしてみましょう。 において、u_nがu_n－1,v_n－1にしか依存していないことに注目すると、数列(u_i),(v_i)だけを計算していけば最終的にu_nが得られることがわかります。これを算法として整理し、そのままプログラミングしてみます。

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

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

 さて、一般のn元一次不定方程式a₁x₁ + ⋯ + a_nx_n =cの解法ですが、線形代数の立場から見ると方程式がn次元ユークリッド空間ℝⁿの超平面の方程式＜`a,`x＞=cの形をしていることから、この超平面上の整数の格子点を求めることと同じです。ここではn=3の場合だけを考えることにします。 とし、この特殊解(x₁,x₂,x₃)が与えられたとしましょう。a₁(X₁－x₁)+a₂(X₂－x₂)+a₃(X₃－x₃)=0から ここでa₃ = a₃'d, a₁ = a₁'d'd, a₂ = a₂'d'd, d'=gcd(a₁,a₂)とすると、 (a₃',d')=1よりd'|( X₃ － x₃ )ですから、 と書くことができます。これを代入して、 X₁' = X₁ － x₁, X₂' = X₂ － x₂とおくと となります。この一般解はa₁' X₁' + a₂' X₂' = 1の一般解に－ a₃'uを掛けたものになりますのでこの特殊解を(s₁,s₂)とすると、 ここまでの結果を整理すると求める一般解 が得られます。この一般解が実際に求める方程式の解になっていること、つまり整数で3次元空間内の平面a₁ X + a₂ Y + a₃ Z = d上にあることは容易に確かめることができます。これを命題の形にまとめます。

 この命題から３元一次不定方程式の解法はその特殊解を求めることに帰着されます。これはイデアルとして(a₁, a₂, a₃) = ((a₁, a₂), a₃)となることからd'=gcd(a₁, a₂)とおくとd = gcd(a₁, a₂, a₃) = gcd(d', a₃)となり、 の特殊解を(x₀, y₀)とし、 の特殊解を(u₀, z₀)とすれば、(x₀ u₀, y₀ u₀, z₀)が求める特殊解となり、拡張互除法を２回適用することで計算できます。この方法をプログラミングしてみましょう。

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回適用することで求めることができます。一般解の構造については連立一次不定方程式の解法とともに線形代数を使って説明するのが適切ですので、ここでは詳細を述べないことにします。

参考文献