首页 > Python算法 > 数论算法 阅读数:34

中国余数定理

在公元 4 世纪中国数学家孙子在《孙子数经》中提出了以下问题:一个数字被 3 除余 2,被 5 除余 5,被 7 除余 2,这个数字是什么?孙子并没有给出解决这种问题的方法,所以我们将要学习的方法是由印度数学家阿耶波多在公元 6 世纪提出的。

余数问题的正式形式描述如下:

给予一元性同余方程组 S,其中为 m1,m2,m3,…,mn 大于 1 的整数, a1,a2,a3,…,an 为任意整数。
求方程组 S 的有解条件,并在有解的情况下求解。

1. 算法介绍

中国余数定理表示,如果 m1,m2,…,mi,…,mn 为互质的整数则 S 有解,否则无解。互质意味着 m1,m2,…,mi,…,mn 除 1 没有其他的公约数。比如(3,5,7)互质,(6,4,7)不互质,因为 6 与 4 的公约数是 2。

在有解的情况下,中国余数定理给出了解的形态表达公式。中国余数定理算法过程如下:
  1. 令 M=m1,m2,…,mi,…,mn ;
  2. 令 Mi=M/mi
  3. 寻找 c1,c2,…,ci,…,cn ,使 ciMi=1(mod mi);
  4. 方程组的解为:

让我们通过解决示例来更好地理解上述算法,我们要找到所有满足条件的 x。

① M 为除数的乘积。
M=3×5×7=105

② Mi 为除了 mi 所有其他除数的乘积
M1=5 × 7=35
M2 =3 × 7=21
M3 =3 × 5=15

③ 寻找 c1,使 c1Mc1≡1(mod m1),也就是 c135≡1(mod 3),我们得出 c1=2。

④ 寻找 c2,使 c2M2 ≡1(mod m2),也就是 c221≡1(mod 5),我们得出 c2 =1。

⑤ 寻找 c3,使 c3M3≡1(mod m3),也就是 c315≡1(mod 7),我们得出 c3 =1。

⑥ 计算解
x≡a1M1c1+a2M2c2+a3M3c3(modM)
x≡2×35×2+3×21×1+2×15×1(mod105)
x≡233(mod105)
x≡23(mod105)

⑦ 检查。
我们得到的解为 x=23+105k,k 为任意整数。首先检查 x=23:23 除 3 余 2,23 除 5 余 3,23 除 7 余 2。正确。接着令 k=1 检查 x=23+105=128:128 除 3 余 2,23 除 5 余 3,23 除 7 余 2。正确。

2. 算法证明

1) 首先证明有解条件:如果 m1,m2,…,mi,…,mn 为互质的整数,则方程式 S 有解。

证明如下:

我们必须确保算法的第三步 ciMi≡1(mod mi) 成立,必要的条件为 gcd(mi,Mi)=1,也就是说,mi 与 Mi 的最大公约数为 1。如果 gcd(mi,Mi)≠1,比如 15ci≡1(mod 5),任何 c都满足不了这个式子,毕竟任何数都不能“阻止”15 整除 5。

因为我们需要所有的 (mi,Mi) 互质,所以我们需要所有的 m1,m2,…,mi,…,mn 互质。证明完毕。

2) 接下来我们证明解的形态表达公式:

x=a1M1c1+a2M2c2+…+aiMici+…+anMncn(mod M)

证明如下:

考察 a1M1c1+a2M2c2+…+aiMici+…+anMncn(mod M),由于除 Mi 以外所有的 M1,M2,…,Mn 都有 mi 为因数,因此除 aiMici 所有的 a1M1c1,a2M2c2, …,anMncn 都能够被 mi 整除,因此我们可以表达以上式子为:

a1M1c1+a2M2c2+…+aiMici+…+anMncn(mod mi)+0+0+…+aiMici+…+0(mod mi)


3) 另外,因为我们设定 Mici≡1(mod mi),所以可以继续简化以上式子为:

aiMici(mod mi)=a1x1(mod mi)=ai(mod mi)

这就是我们一开始的定义,x 除 mi 余 ai。公式成立。

3. 算法代码

中国余数算法需要利用扩展欧几里得算法,去计算。该算法我们需要寻找,使 ciMi≡1(mod mi)。ciMi≡1(mod mi) 也可以被表达为 ciMi−kmi=1,这不就是在讲到的贝祖等式吗(因为 Mi 与 mi 的最大公约数为1)?因此我们可以通过扩展欧几里得算法找到满足条件的 ci

代码如下:
#扩展欧几里得算法,a与b为自然数
defextendedEuclidean(a,b):
    if (b == 0):             #边界条件
        return 1, 0            #输出1,0因为g = 1*a+1*0
    else:
        m,n = extendedEuclidean(b, a % b)
        quotient = a//b
        return n, m - (n * quotient)    #输出整数m,n,使gcd(a,b) =ma + nb
#中国余数算法,参数为两个二维数组
#reminder = [a1,a2,a3,…,an]
#divisor = [m1,m2,…,mi,…,mn]
def CRT(reminder, divisor):
    X=0               #声明变量X用来储存方程解
    M = 1              #M为除数的乘机
Mi = []            #Mi用来储存M1,M2,…,Mi,…,Mn
ci = []            #ci用来储存c1,c2,…,ci,…,cn
    for x in divisor:
        M = M*x
    for x in divisor:
        Mix = M//x
        Mi.append(Mix)
        ci.append(extendedEuclidean(Mix,x)[0])
    fori in range(0,len(reminder)):
        X+=reminder[i]*Mi[i]*ci[i]   #计算X
    return X%M,M
运行如下测试代码:

reminder = [2,3,2]
divisor = [3,5,7]
print(CRT(reminder, divisor))

解为 x=23(mod105),得到输出:

(23,105)