跳转至

中国剩余定理

引入

「物不知数」问题:有物不知其数,三三数之剩二,五五数之剩三,七七数之剩二.问物几何?

即求满足以下条件的整数:除以 3 余 2,除以 5 余 3,除以 7 余 2.

该问题最早见于《孙子算经》中,并有该问题的具体解法.宋朝数学家秦九韶于 1247 年《数书九章》卷一、二《大衍类》对「物不知数」问题做出了完整系统的解答.上面具体问题的解答口诀由明朝数学家程大位在《算法统宗》中给出:

三人同行七十希,五树梅花廿一支,七子团圆正半月,除百零五便得知.

2×70+3×21+2×15=233=2×105+23,故答案为 23.

定义

中国剩余定理 (Chinese Remainder Theorem, CRT) 可求解如下形式的一元线性同余方程组(其中 n1,n2,⋯,nk 两两互质):

{x≡a1(modn1)x≡a2(modn2)⋮x≡ak(modnk)

上面的「物不知数」问题就是一元线性同余方程组的一个实例.

过程

  1. 计算所有模数的积 n;
  2. 对于第 i 个方程:
    1. 计算 mi=nni;
    2. 计算 mi 在模 ni 意义下的 逆元 mi−1;
    3. 计算 ci=mimi−1(不要对 ni 取模).
  3. 方程组在模 n 意义下的唯一解为:x=∑i=1kaici(modn).

实现

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
LL CRT(int k, LL* a, LL* r) {
  LL n = 1, ans = 0;
  for (int i = 1; i <= k; i++) n = n * r[i];
  for (int i = 1; i <= k; i++) {
    LL m = n / r[i], b, y;
    exgcd(m, r[i], b, y);  // b * m mod r[i] = 1
    ans = (ans + a[i] * m * b % n) % n;
  }
  return (ans % n + n) % n;
}
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
def CRT(k, a, r):
    n = 1
    ans = 0
    for i in range(1, k + 1):
        n = n * r[i]
    for i in range(1, k + 1):
        m = n // r[i]
        b = y = 0
        exgcd(m, r[i], b, y)  # b * m mod r[i] = 1
        ans = (ans + a[i] * m * b % n) % n
    return (ans % n + n) % n

证明

我们需要证明上面算法计算所得的 x 对于任意 i=1,2,⋯,k 满足 x≡ai(modni).

当 i≠j 时,有 mj≡0(modni),故 cj≡mj≡0(modni).又有 ci≡mi⋅(mi−1modni)≡1(modni),所以我们有:

x≡∑j=1kajcj(modni)≡aici(modni)≡ai⋅mi⋅(mi−1modni)(modni)≡ai(modni)

即对于任意 i=1,2,⋯,k,上面算法得到的 x 总是满足 x≡ai(modni),即证明了解同余方程组的算法的正确性.

因为我们没有对输入的 ai 作特殊限制,所以任何一组输入 {ai} 都对应一个解 x.另外,若 x≠y,则总存在 i 使得 x 和 y 在模 ni 下不同余.故系数列表 {ai} 与解 x 之间是一一映射关系,方程组总是有唯一解.

解释

下面演示 CRT 如何解「物不知数」问题.

  1. n=3×5×7=105;
  2. 三人同行 七十 希:n1=3,m1=n/n1=35,m1−1≡2(mod3),故 c1=35×2=70;
  3. 五树梅花 廿一 支:n2=5,m2=n/n2=21,m2−1≡1(mod5),故 c2=21×1=21;
  4. 七子团圆正 半月:n3=7,m3=n/n3=15,m3−1≡1(mod7),故 c3=15×1=15;
  5. 所以方程组的唯一解为 x≡2×70+3×21+2×15≡233≡23(mod105).(除 百零五 便得知)

Garner 算法

CRT 的另一个用途是用一组比较小的质数表示一个大的整数.

例如,若 a 满足如下线性方程组,且 a<∏i=1kpi(其中 pi 为质数):

{a≡a1(modp1)a≡a2(modp2)⋮a≡ak(modpk)

我们可以用以下形式的式子(称作 a 的混合基数表示)表示 a:

a=x1+x2p1+x3p1p2+…+xkp1…pk−1

Garner 算法 将用来计算系数 x1,…,xk.

令 rij 为 pi 在模 pj 意义下的 逆:

pi⋅ri,j≡1(modpj)

把 a 代入我们得到的第一个方程:

a1≡x1(modp1)

代入第二个方程得出:

a2≡x1+x2p1(modp2)

方程两边减 x1,除 p1 后得

a2−x1≡x2p1(modp2)(a2−x1)r1,2≡x2(modp2)x2≡(a2−x1)r1,2(modp2)

类似地,我们可以得到:

xk=(…((ak−x1)r1,k−x2)r2,k)−…)rk−1,kmodpk
实现
1
2
3
4
5
6
7
8
for (int i = 0; i < k; ++i) {
  x[i] = a[i];
  for (int j = 0; j < i; ++j) {
    x[i] = r[j][i] * (x[i] - x[j]);
    x[i] = x[i] % p[i];
    if (x[i] < 0) x[i] += p[i];
  }
}
1
2
3
4
5
6
7
for i in range(0, k):
    x[i] = a[i]
    for j in range(0, i):
        x[i] = r[j][i] * (x[i] - x[j])
        x[i] = x[i] % p[i]
        if x[i] < 0:
            x[i] = x[i] + p[i]

该算法的时间复杂度为 O(k2).实际上 Garner 算法并不要求模数为质数,只要求模数两两互质,我们有如下伪代码:

Chinese Remainder Algorithm cra⁡(v,m):Input: m=(m0,m1,…,mn−1), mi∈Z+∧gcd(mi,mj)=1 for all i≠j,v=(v0,…,vn−1) where vi=xmodmi.Output: xmod∏i=0n−1mi.1for i from 1 to (n−1) do2Ci←(∏j=0i−1mj)−1modmi3x←v04for i from 1 to (n−1) do5u←(vi−x)⋅Cimodmi6x←x+u∏j=0i−1mj7return (x)

可以发现在第六行中的计算过程对应上述混合基数的表示.

应用

某些计数问题或数论问题出于加长代码、增加难度、或者是一些其他原因,给出的模数:不是质数!

但是对其质因数分解会发现它没有平方因子,也就是该模数是由一些不重复的质数相乘得到.

那么我们可以分别对这些模数进行计算,最后用 CRT 合并答案.

下面这道题就是一个不错的例子.

洛谷 P2480 [SDOI2010] 古代猪文

给出 G,n(1≤G,n≤109),求:

G∑k∣n(nk)mod999 911 659

首先,当 G=999 911 659 时,所求显然为 0.

否则,根据 欧拉定理,可知所求为:

G∑k∣n(nk)mod999 911 658mod999 911 659

现在考虑如何计算:

∑k∣n(nk)mod999 911 658

因为 999 911 658 不是质数,无法保证 ∀x∈[1,999 911 657],x 都有逆元存在,上面这个式子我们无法直接计算.

注意到 999 911 658=2×3×4679×35617,其中每个质因子的最高次数均为一,我们可以考虑分别求出 ∑k∣n(nk) 在模 2,3,4679,35617 这几个质数下的结果,最后用中国剩余定理来合并答案.

也就是说,我们实际上要求下面一个线性方程组的解:

{x≡a1(mod2)x≡a2(mod3)x≡a3(mod4679)x≡a4(mod35617)

而计算一个组合数对较小的质数取模后的结果,可以利用 卢卡斯定理.

扩展:模数不互质的情况

两个方程

设两个方程分别是 x≡a1(modm1)、x≡a2(modm2);

将它们转化为不定方程:x=m1p+a1=m2q+a2,其中 p,q 是整数,则有 m1p−m2q=a2−a1.

由 裴蜀定理,当 a2−a1 不能被 gcd(m1,m2) 整除时,无解;

其他情况下,可以通过 扩展欧几里得算法 解出来一组可行解 (p,q);

则原来的两方程组成的模方程组的解为 x≡b(modM),其中 b=m1p+a1,M=lcm(m1,m2).

多个方程

用上面的方法两两合并即可.

习题