中国剰余定理のプログラム

　中国剰余定理のプログラム（連立合同式を解くプログラム）を書いてみた。

中国剰余定理とは？

2で割ると1余り、3で割ると2余り、5で割ると3余る数なーんだ？
というような問題を考える。

この問題を式で表すと、以下のような連立合同式になる。

x = 1 (mod 2)
x = 2 (mod 3)
x = 3 (mod 5)

中国剰余定理によると、
法が互いに素であれば、法の総積を法としてxは一意に定まる。

以下ではそのようなxを効率的に求めることを考える。

解法

法が互いに素な場合の解法は、以下のとおり。

x = a₀ (mod m₀)
x = a₁ (mod m₁)

を解きたいとする。
(m₀, m₁) = 1であるため、

p m₀ + q m₁ = 1

となるような、整数p, qが存在する。

x = a₀ q m₁ + a₁ p m₀

とすれば、xは上の方程式の解である。

これは以下のように確認できる。

x = a₀ (1 - p m₀) + a₁ p m₀ = a₀ - a₀ p m₀ + a₁ p m₀ = a₀ (mod m₀)
x = a₀ q m₁ + a₁ (1 - q m₁) = a₀ q m₁ + a₁ - a₁ q m₁ = a₁ (mod m₁)

次にm₀, m₁が互いに素ではない場合を考える。この場合は、m₀、m₁が互いに素になるように素因数を片方に寄せればよい。

g = (m₀, m₁)

とし、

g = p0^n0 p1^n1 p2^n2 .... pk^nk ...

と素因数分解されるとする。

素因数pkの個数を考える。m₀の方がm₁よりpkを多く持っていれば、pkはすべてm₀に寄せる。そうでない場合は、pkはすべてm₁に寄せる。

これをすべての素因数について行えば、最終的に(m0, m₁) = 1となる。
よってm0, m₁が互いに素な場合に帰着する。

プログラム

def gcd(a, b):
    if b == 0:
        return a
    return gcd(b, a%b)

# return [g, x, y]
# g = gcd(a, b)
# x, y satisfies a x + b y = g
def extgcd(a, b):
    if b == 0:
        return [a, 1, 0]
    g, x, y = extgcd(b, a%b)
    return [g, y, x - a/b * y]

# eq0: x = a0 (mod m0)
# eq1: x = a1 (mod m1)
# returns [xt, mod] such that x = xt + k mod for integer k.
def crt(eq0, eq1):
    a0, m0 = eq0
    a1, m1 = eq1

    g = gcd(m0, m1)

    if a0 % g != a1 % g:
        raise Exception("x doesn't exist.")

    if g > 1:
        m0 /= g
        m1 /= g

        while True:
            gt = gcd(m0, g)
            if gt == 1:
                break
            m0 *= gt
            g /= gt
        
        m1 *= g

        a0 %= m0
        a1 %= m1

    g, p, q = extgcd(m0, m1)
    
    x = a0 * q * m1 + a1 * p * m0
    mod = m0 * m1
    x = x % mod

    return [x, mod]

テスト

テストコード。

# for debug use only
def test():
    for x in range(1000):
        for m1 in range(1, 1000):
            for m2 in range(1, 1000):
                a1 = x % m1
                a2 = x % m2

                y, mod = crt((a1, m1), (a2, m2))
                
                assert x % mod == y % mod

    print "test success"

サンプル

標準入力から連立合同方程式を受け取って、解くプログラム。
解を10個出力するようにしている。

# read system of linear congruence
n = int(raw_input('input # of system:'))
eqs = []
for _ in xrange(n):
    a, m = map(int, raw_input().split())
    eqs.append((a, m))

# solve the system
x, mod = reduce(crt, eqs, (0, 1))

# print the answer
for _ in range(10):
    print x
    x += mod

まず、 連立合同式
x = 1 (mod 2)
x = 2 (mod 3)
x = 3 (mod 5)
を解いてみる。

input # of system:3
1 2
2 3
3 5
x = 
23
53
83
113
143
173
203
233
263
293

次に、法が互いに素ではない連立式を解いてみる。

input # of system:3
15471 102102
357813 504504
24691357 98765432
x = 
123456789
105883678906461
211767234356133
317650789805805
423534345255477
529417900705149
635301456154821
741185011604493
847068567054165
952952122503837