もう一度中国剰余定理

数学

こんにちは。またまた数学です。暗号と数学は密接で切り離すことができないのでもう少し続けます。

中国剰余定理とは

以前の記事で取り上げていました。

整数\(m_1 , m_2, \cdots m_k\)が互いに素ならば、任意に与えられる整数\(a_1 , a_2, \cdots a_k\)に対して

\(x \equiv a_1 \bmod m_1 \)

\(x \equiv a_2 \bmod m_2 \)

\( \cdots \)

\(x \equiv a_k \bmod m_k \)

を満たす\(x\)は、\(m_1 m_2 \cdots m_k\)を法として一意に存在するというものでした。

そこで、実装をしようと思ったところモジュラ逆元を求める必要が出てきて、中途半端になっていました。今回はプログラムを用いて\(x\)の値を求めるところまで進めていきます。

ライブラリを用いて解く

毎度のPythonですが、Pythonのライブラリには本当に便利なものが多くあるので、それを使って求めます。

import sympy.ntheory.modular

if __name__ == '__main__':
  print( sympy.ntheory.modular.crt( [ 3, 5, 7 ], [ 2, 3, 3 ] ) )

実行結果は次になります。 

(38, 105)

ライブラリを使ってしまえば、簡単に求められます。しかし、数学的な知識として、このライブラリによってどのような値が求められているかについての理解が必要になります。

Gaussのアルゴリズム

次にGaussのアルゴリズムについて紹介します。

\(M=m_1m_2\cdots m_k\)として、\(M_i=m_1m_2\cdots m_{i-1} m_{i+1} \cdots m_k = M / m_i \)とします。また、\(y_i M_i \equiv 1 (\bmod m_i ) \)となる\(y_i \equiv {M_i} ^{-1} (\bmod m_i) \)が存在し、これを用いて次の式

\[x \equiv a_1 M_1 y_1 + a_2 M_2 y_2 + \cdots + a_k M_k y_k (\bmod M)\]

で求めることができます。

具体的な数値でやってみます。

\(\begin{eqnarray} x & \equiv & 2 (\bmod 3 ) \\ x & \equiv & 3 (\bmod 5) \\ x & \equiv & 3 (\bmod 7) \end{eqnarray} \)

を満たす\(x\)を求めていきます。

\(M=3 \times 5 \times 7 = 105\)です。

\(M_1 = 105/3=35\)となり、\(y_1=2\)とすると\(M_1 y_1 \equiv 2 \times 35 \equiv 70 \equiv 1 (\bmod 3)\)です。

\(M_2 = 105/5=21\)となり、\(y_2=1\)とすると\(M_2 y_2 \equiv 1 \times 21 \equiv 21 \equiv 1 (\bmod 5)\)です。

\(M_3 = 105/7=15\)となり、\(y_3=1\)とすると\(M_3 y_3 \equiv 1 \times 15 \equiv 15 \equiv 1 (\bmod 7)\)です。

よって、\(x \equiv 2 \times 35 \times 2 + 3 \times 21 \times 1 + 3 \times 15 \times 1 \equiv 248 \equiv 38 (\bmod 105)\)が求められます。

一応、証明につながる考え方を書いておきます。

上で求めている\(x\)を\(\bmod 3\)で考えます。

\(x \equiv 2 \times 35 \times 2 + 3 \times 21 \times 1 + 3 \times 15 \times 1 (\bmod 3)\)を考えます。

第2項と第3項はいずれも\(21=3\times7,15=3\times5\)が含まれているので\(3\)の倍数となり、\(\bmod 3\)では余りは\(0\)となります。

また、\(35\times2\equiv1(\bmod 3)\)となるので、全体として\(x\equiv 2(\bmod 3)\)が得られます。

同様に\(x\equiv 3(\bmod 5)\)，\(x\equiv 3(\bmod 7)\)も得られます。

この手順をプログラムで書いていきます。

def inverse( a, p ):
  e = p - 2
  k = a
  result = 1
  while e > 0:
    bit = e & 1
    if bit==1:
      result = result*k % p
    k = k*k % p
    e = e >> 1
  return result

def gauss( mlist, alist ):
  M = 1
  for m in mlist:
    M *= m
  x = 0
  for m, a in zip(mlist, alist):
    tmp = M // m
    x = ( x + a * tmp * inverse( tmp, m ) ) % M
  return x, M

import sympy.ntheory.modular

if __name__ == '__main__':
  print( gauss( [ 3, 5, 7 ], [ 2, 3, 3 ] ) )

実行結果は次になります。 

(38, 105)

Garnerのアルゴリズム

次にGarnerのアルゴリズムを取り上げます。先に具体的な数値で考えます。

\(x \equiv 2 (\bmod 3)\)は、\(x=2+3k_1(k_1は整数)\)と書けます。

次に、\(x \equiv 3(\bmod 5)\)も満たしています。

\[\begin{eqnarray}  2+3k_1 & \equiv & 3(\bmod 5) \\ 3k_1 & \equiv & 3 – 2(\bmod 5) \\ k_1 & \equiv & 3^{-1}(3-2) (\bmod 5) \\ k_1 & \equiv & 2 (\bmod 5) \end{eqnarray} \]

となって、\(x=2+ 3\times 2 + (3 \times 5)k_2=8+15k_2(k_2は整数)\)と書けます。

同様に続きを考えることができます。\(x \equiv 3(\bmod 7)\)も満たしています。

\[\begin{eqnarray}  8+15k_2 & \equiv & 3(\bmod 7) \\ 15k_2 & \equiv & 3 – 8(\bmod 7) \\ k_2 & \equiv & {15}^{-1}(3-8) (\bmod 7) \\ k_2 & \equiv & 1\times 2 (\bmod 7) \\k_2 & \equiv & 2 (\bmod 7) \end{eqnarray} \]

となって、\(x=8+ 15\times 2 + (3 \times 5 \times 7)k_3=38+105k_3(k_3は整数)\)と書けます。

このように合同式を順に２つずつ組み合わせて、ともに満たす合同式を作っていくアルゴリズムがGarnerのアルゴリズムです。

それではプログラムで書いていきます。

def garner( mlist, alist ):
  m, x = mlist[0], alist[0]
  i = 1
  while i<len(mlist):
    x = x + m * inverse( m, mlist[i] ) * ( alist[i] - x )
    m *= mlist[i]
    i += 1
  return x % m, m

if __name__ == '__main__':
  print( garner( [ 3, 5, 7 ], [ 2, 3, 3 ] ) )

結果は同じになることはわかっていますが、念のため実行しておきます。 

(38, 105)

拡張ユークリッドの互除法を用いて不定方程式の解を求める

はじめに、

\(\begin{eqnarray} x & \equiv & 2 (\bmod 3 ) \\ x & \equiv & 3 (\bmod 5) \end{eqnarray} \)

を満たす\(x\)を求めていきます。

これは、不定方程式\(3u+5v=1\)の解\(u=2,v=-1\)を用いて\(x \equiv 3\times 3 \times 2 + 2 \times 5 \times (-1) \equiv 8 (\bmod 15)\)と求められます。

次に、

\(\begin{eqnarray} x & \equiv & 8 (\bmod 15 ) \\ x & \equiv & 3 (\bmod 7) \end{eqnarray} \)

を満たす\(x\)を求めていきます。

これは、不定方程式\(15u+7v=1\)の解\(u=1,v=-2\)を用いて\(x \equiv 3\times 15 \times 1 + 8 \times 7 \times (-2) \equiv 38 (\bmod 105)\)が求められます。

このように\(\bmod\)の値に着目して不定方程式を解き、二つの合同式をともに満たす値\(x\)を順に求めていくことにより、すべての合同式を満たす\(x\)を求めることができます。

プログラムは次になります。

def euclid( a, b ):
  u1, u = 1, 0
  v1, v = 0, 1
  while b:
    q = a // b
    a, b = b, a-q*b
    u1, u = u, u1-q*u
    v1, v = v, v1-q*v
  return u1, v1

def solve_equation( mlist, alist ):
  m = mlist[0]
  x = alist[0]
  i = 1
  while i<len(mlist):
    u, v = euclid( m, mlist[i] )
    x = alist[i] * m * u + x * mlist[i] * v
    m *= mlist[i]
    x %= m
    i += 1
  return x, m

if __name__ == '__main__':
  print( solve_equation( [ 3, 5, 7 ], [ 2, 3, 3 ] ) )

念のため実行結果を載せておきます。 

(38, 105)

今回はこれでおしまいにします。それではまた。

Posted by 春日井 優