晴耕雨読

working in the fields on fine days and reading books on rainy days

中国人剰余定理

中国人剰余定理とは整数の剰余に関する定理の一つです。 ここでは、連立合同方程式を中国人剰余定理を使ったPythonプログラムで解く方法について説明します。

  • 中国人剰余定理

    $k$ 個の整数 $n_1, n_2, …, n_k$ が互いに素ならば、任意の整数 $a_1, a_2, …, a_k$ に対して、次の連立合同方程式を満たす整数 $x$ が $n_1 n_2 \cdots{} n_k$ を法として一意的に存在する。

    唯一の解 $x$ は、$1 \le i \le k$ について とおくと、次の式で求めることができる。

中国人剰余定理で連立合同方程式 $x \equiv a_i \pmod{n_i} \;\;\text{for}\; i = 1,…,k$ を解くためのアルゴリズムは以下の通りです。

  • アルゴリズム

    1. $N = n_1 n_2 \cdots{} n_k$ を求める。
    2. に対して、整数 $n_i$ と $b_i = N / n_i$ を求める。この2つの整数は互いに素となる。
    3. $n_i$ を法とする $b_i$ の逆元 $b'_i = b_i^{-1} \pmod{n_i}$ を求める。
    4. 次の式を計算して解 $x$ を求める。

    5. 最小の解 $x$ を求める。

このアルゴリズムをPythonで書くと次のようになります。

from functools import reduce

def xgcd(a, b):
    x0, y0, x1, y1 = 1, 0, 0, 1
    while b != 0:
        q, a, b = a // b, b, a % b
        x0, x1 = x1, x0 - q * x1
        y0, y1 = y1, y0 - q * y1
    return a, x0, y0

def modinv(a, m):
    g, x, y = xgcd(a, m)
    if g != 1:
        raise Exception('modular inverse does not exist')
    else:
        return x % m

def chinese_remainder(a, n):
    # a := [a1, a2, ..., ak]
    # n := [n1, n2, ..., nk]
    total = 0
    prod = reduce(lambda x, y: x*y, n)
    for n_i, a_i in zip(n, a):
        b_i = prod // n_i
        total += a_i * b_i * modinv(b_i, n_i)
    return total % prod

連立合同方程式を中国人剰余定理で解く問題として『孫子算経』の問題を使います。 つまり、次の連立合同方程式を解いてみます。

先ほどのプログラムを使って計算してみると、解は 23 となります。

# 連立合同方程式を中国人剰余定理で解く
a = [2, 3, 2]
n = [3, 5, 7]
print(chinese_remainder(a, n))
# => 23

『孫子算経』の問題の答えと一致しました。

以上です。

参考文献