Search on the blog

2015年2月15日日曜日

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

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

中国剰余定理とは?
2で割ると1余り、3で割ると2余り、5で割ると3余る数なーんだ?
というような問題を考える。

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

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

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

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

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

x = a0 (mod m0)
x = a1 (mod m1)

を解きたいとする。
(m0, m1) = 1であるため、

p m0 + q m1 = 1

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

x = a0 q m1 + a1 p m0

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

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

x = a0 (1 - p m0) + a1 p m0 = a0 - a0 p m0 + a1 p m0 = a0 (mod m0)
x = a0 q m1 + a1 (1 - q m1) = a0 q m1 + a1 - a1 q m1 = a1 (mod m1)

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

g = (m0, m1)

とし、

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

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

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

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

プログラム
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

0 件のコメント:

コメントを投稿