文章总览:穷高极天,亢盈难久

中国有句古话:“有物不知其数,三三数之剩二,五五数之剩三,七七数之剩二。问物几何?”上过小学一年级的我们都知道,这是非常经典的同余方程组,而它的解法,就是这篇文章要讲的中国剩余定理 Chinese Remainder Theorem

本文内容限定于下:

  1. 中国剩余定理证明
  2. 中国剩余定理实现
  3. 拓展中国剩余定理证明
  4. 拓展中国剩余定理实现

推导其一:威制八毒,灭却炎烟

从实例出发的推导

那么从开篇提到的“物不知数”问题出发,我们可以得到下面的同余方程组:

{x2(mod3)x3(mod5)x2(mod7)

其实这个题目很早就有了答案,“三人同行七十稀,五树梅花廿一支,七子团圆正半月,除百零五便得知。”这是来自于明朝数学家程大位在《算法统宗》中的解答口诀。其含义即是关于这个问题的通解形式,其中第一个字的表示模数,最后一个词中的数字表示每一个同余方程余数要乘以的系数,最后加起来模 105 就是答案,可以用下面的式子来表示这个同余方程组的解:

x=(2×70+3×21+2×15)mod105=23

关于这个通解的形式,和余数的关系不大,也就是说我们可以再再抽象刚刚的问题,得到下面的问题形式和通解形式,可以观察一下中国剩余定理能处理的基本问题:

ProblemSolution{xa1(mod3)xa2(mod5)xa3(mod7)x=(70a1+21a2+15a3)mod105

实际上模数并不是只能用素数,直接瞪眼法的话也看不出来有什么特殊的数字,我最多只能看出来最后的模数是所有模数的乘积 3×5×7,那么接下来正式介绍一下中国剩余定理的内容。

所谓中国剩余定理,正式的内容就是对于下面的一个同余方程组,可以用 M 表示所有模数的乘积,Mi 表示 M/mi 的值,然后可以表示为下面的通解形式:

ProblemSolution{xa1(modm1)xa2(modm2)xak(modmk)x=(i=1kaiMiMi1)modM

其中 Mi1 表示 Mi 在模 mi 意义下的逆元。我们可以用这个通解的公式和刚刚得到诗中提到的的同余方程的结果来对比一下:

M=3×5×7=105idmiMiMi1(modmi)MiMi1×ai133535×21(mod3)70×a1252121×11(mod5)21×a2371515×11(mod7)15×a3x=(70a1+21a2+15a3)modM=(70×2+21×3+15×2)mod105=23

我们会发现一个惊人的事实,刚刚的通解形式完美匹配上了,接下来就是证明这个通解形式正确性的时候了。

通解正确性的证明

那么想要证明通解形式的正确性,就要从这个同余方程的两个方面下手,第一个是解的存在性,要保证一定有解;第二个是解的唯一性,即所有的解都是通解形式表示出来的这种解,或者只有通解形式标识出来的着一种解。好像是一个意思?(雾

重复一遍同余方程的内容,下面的同余方程中,保证模数两两互素

{xb1(modm1)xb2(modm2)xbk(modmk)

首先进行存在性的证明。我们令 A=i=1kbiMiMi1,那么 aZ,aA(modM) 所有这些整数 a 都得是方程组的解,这样就可以证明存在性。

对于所有的 i[1,k],根据 M 的定义都显然有 miM。并且还有 j[1,k],ij 根据 Mj=Mmj 的定义,可以得到 i,j 满足 miMj 也就是 Mj0(modmi),另外还需用到关于 Mi1 的定义,其实就是 Mi 在模 mi 下的逆元,即 MiMi11(modmi)。有了这些就可以证明关于 a 是方程组的解了。

aAi=1kbiMiMi1(modmi)(miM)biMiMi1+j=1ijkbjMjMj1(modmi)biMiMi1(modmi)(miMj)bi(modmi)(MiMi11)

然后就发现 a 满足于方程组中所有的式子,也就是说 a 就是方程的解,即按照通解得出的解一定是方程组的解的一部分。然后呢接下来就是要确定,这个方程组只有通解表示出来的着一类解即证明解的唯一性。

考虑反证法。设 Y 是与通解表现形式不同的方程组的解,则 i[1,k] 都满足 Ybi(modmi)。而根据刚刚解的存在性的证明,我们也可以得到 Abi(modmi),那么根据同余的性质,可以轻松得到 YA(mod[m1,m2,,mk])。而由于模数 mi 两两互素,有 [m1,m2,,mk]=i=1kmi=M,也就是说最后可以得到 YA(modM) 那么和通解形式表现出来的解形式相同,与假设矛盾,唯一性得证。

实现其一:幽明变化,自在蟠跃

那么我们就从模板题开始,LightOJ 1319 - Monkey Tradition 就是一个非常板的 CRT 题目。

题目大意

有一个名为“ Bamboo Climbing”的传统游戏。游戏规则如下:

  1. N 个猴子玩这个游戏,并且有 N 个等高的竹子。其高度为 L 米。
  2. 每只猴子站在竹子前面,每只猴子被分配一个不同的竹子。
  3. 吹口哨时,猴子开始爬竹子,并且在整个游戏过程中不允许它们跳到另一只竹子上。
  4. 由于它们是猴子,因此通常通过跳跃来爬。并且在每次跳跃中,第 i 个猴子都可以精确地跳跃 pi 米( pi 是素数)。当猴子再跳一次可能会使他离开竹子,他就停止了跳跃,他报告了剩余的长度 ri
  5. 在比赛之前,每只猴子都被分配了一个不同的 pi
  6. ri 最低的猴子获胜。

现在,他们已经找到了之前的比赛信息,但不幸的是,他们还没有找到竹子的高度。更确切地说,他们知道 N,全部 pi 和相应的 ri,但不知道 L。因此,您挺身而出,发现任务很艰巨,因此,您想从给定的信息中找到 L

解题思路

显然的根据题目给出的信息,我们可以列出下面一个关于 L 的方程组:

{Lr1(modp1)Lr2(modp2)Lrn(modpn)

而且出题人非常良心的告诉我们 pi 全都是质数,也就是满足上面介绍中国剩余定理时候的模数两两互素这一要求,所以显然的这就是一个裸的中国剩余定理问题。


我们根据上面的通解形式,就可以轻松计算出最终的 L。不过问题是,我们唯一不会的是计算关于 Mi 在模 mi 意义下的逆元。对于这道题来说,所有 pi 都是质数,所以显然可以用费马小定理搞定,但事实上中国剩余定理的问题形式里面只提到了 mi 两两互素,并不强制要求必须得是素数。所以现在的目的是找到一个快速计算逆元的方式。

那么费马小用不了,其实还有一个利器就是欧拉定理,根据欧拉定理有 aφ(m)1(modm) 其中只要满足 am 即可,泛用性比费马小高了不知道多少倍,这边一看,要计算 Mi 在模 mi 意义下的逆元,显然的 Mimi,这就可以用欧拉定理快速幂解决了。


然后除了欧拉定理,其实还有泛用性更高的做法,就是利用拓展欧几里得解决。下面阐述一下为什么可以用拓展欧几里得来解决这个求逆元的问题。关于求逆元的故事,需要从同余方程讲起,所谓同余方程,其实就是组成上面同余方程组的部分,形如 axb(modn) 的方程可以被称作同余方程。那么其实不难发现逆元就是形如 aa11(modm) 的同余方程。不过暂时不用管,我们的目标是先用拓展欧几里得解决常规形式的同余方程。

那么众所周知,拓展欧几里得解决的是形如 ax+by=c 的问题,并且我们知道这个不定方程有整数解的充分必要条件是 gcd(a,b)c

又众所周知,取模的问题都可以看做一个循环的解,也就是说上面的基础同余方程可以表示为 ax+nk=b 的问题,其中 x,k 为未知数且都是整数,这样就可以表述一个同余方程的结果,根据上面的众所周知,可以知道同样的这个不定方程有整数解的充要条件是 gcd(a,n)b

这时候注意到如果用拓展欧几里得求出 x0,k0 满足 ax0+bk0=gcd(a,b),我们就可以通过等式的形式先同除以 gcd(a,b) 再同乘以 b 即使我们需要的一组解。

那么回到一开始的,问题,现在的同余方程是 MiMi11(modmi),根据上面的结论,用拓展欧几里得得到对应的 x0 以后,同时除以最小公倍数再乘上 b 即可,而在这里发现 gcd(Mi,mi)=1 且最后的 b=1,也就是说直接沿用拓展欧几里得的结果即可。

代码实现

下面给出这道题两种方法的代码:

cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#include<bits/stdc++.h>
using namespace std;

int n;
int p[20], r[20];

void Input() {
scanf("%d", &n);
for(int i = 1; i <= n; i++) {
scanf("%d%d", &p[i], &r[i]);
}
}

void Exgcd(long long a, long long b, long long &x, long long &y) {
if(!b) {
x = 1, y = 0;
return;
}
Exgcd(b, a % b, y, x);
y -= a / b * x;
}

void Work() {
long long M = 1;
for(int i = 1; i <= n; i++) {
M *= p[i];
}
long long ans = 0;
for(int i = 1; i <= n; i++) {
long long inv, y;
long long Mi = M / p[i];
// 这个是欧拉定理做法,由于 p 是质数 phi(p)=p-1 等价于费马小
// inv = Pow(Mi, p[i] - 2, p[i]); // 快速幂不写了
// 然后是拓欧做法,字面意思
Exgcd(Mi, p[i], inv, y);
inv %= p[i];
ans += r[i] * inv * Mi;
ans %= M;
}
printf("%lld\n", (ans + M) % M);
}

int main() {
int T, ca = 0;
scanf("%d", &T);
while(T--) {
Input();
printf("Case %d: ", ++ca);
Work();
}
return 0;
}

推导其二:奋迅三味,如日空居

两个方程时的特解

其实是说中国剩余定理是一个非常强的性质题,也就是指它的模数 m 两两互素,这个性质钛强了,大多数时候都是不太可能有这么好用的性质在题目里的。也就是说当 m 两两不一定互素的时候,中国剩余定理还有用吗?这时候,就要掏出进阶版的中国剩余定理——拓展中国剩余定理 ExCRT 出山了。

显然的如果想要沿用 CRT 是不太可能的,因为当 m 两两不一定互素的时候,Mimi 就不再一定互素,那么我们计算答案时的重要数据逆元就不再存在,这时候 CRT 就被踢掉了。

那么接下来就要考虑另辟蹊径,我们可以从两个方程开始寻找规律,考虑如下的同余方程组:

{xb1(modm1)xb2(modm2)

这时候要注意 m1,m2 是不互质的,这个方程的解应该如何计算呢?根据同余的定义,我们可以把上面的两个同余方程表示成一个等式 x=k1m1+b1=k2m2+b2,其中 k1,k2Z。小学我们就学过,相似的未知数要放在一起,所以把刚刚的式子移项一下,k1m1k2m2=b2b1 这时候 k1,k2 是未知数,典型的不定方程,它有解的条件就是 gcd(m1,m2)(b2b1)

等式两边同时除以 gcd(m1,m2),设未知数 d=gcd(m1,m2),p1=m1d,p2=m2d 则刚刚的不定方程可以化作 k1p1k2p2=b2b1d 的形式。

拓展欧几里得解决出关于方程 λ1p1λ2p2=1 的解以后,就可以在等式两边同乘 b2b1d 从而得到正确的 (k1,k2)。这时候我们就可以表示出同余方程组的一个解 x=b2b1dλ1m1+b1

向多方程方向拓展

那么我们注意到,同余方程的解其实也是同余形式的,也就是说关于上面两个方程的特解我们只需要找到一个模数 m,使得 x[0,m) 只有一个解满足上面的同余方程。即设如下方程组有特解 x,则找到一个 m 使得通解 x 满足 xx(modm) 恒成立。

{xb1(modm1)xb2(modm2)

如果这样的 m 存在的话,我们就可以在多方程的同余方程组中通过这个性质两两合并方程组最终得到结果。事实上这个 m 的值就是两个同余方程的模数的最小公倍数即 m=lcm(m1,m2)

我们可以简单的证明一下,这里直接设在 [0,lcm(m1,m2)) 的范围中有两个解 x,y。一样的套路,只要能证明这两个解相等,那么解的唯一性就一目了然了。那么基于上面的同余方程组代入 x,y 得到:

{xb1(modm1)xb2(modm2),{yb1(modm1)yb2(modm2)

那么联立两个方程组,并假设 xy,可以得到:

{(xy)0(modm1)(xy)0(modm2)lcm(m1,m2)(xy)

这时候由于 0x,y<lcm(m1,m2),那么 (xy)<lcm(m1,m2) 但是这时候 lcm(m1,m2) 可以整除 xy,这只能说明 xy=0 也即是 x=y 那么唯一性得证,就可以通过两两合并方程组得到结果,这就是拓展欧几里得的过程,接下来可以看题目来做一做。

实现其二:一蹄天水,六虚洪流

其实还是模拟思路即可……P4777 【模板】扩展中国剩余定理(烟斗不带烟了

题目大意

给定 n 组非负整数 ai,bi ,求解关于 x 的方程组的最小非负整数解。

{xb1(moda1)xb2(moda2)xbn(modan)

数据保证有解。

解题思路

模拟刚刚的思路即可,详见代码。

代码实现

cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
// typedef __int128 int128;
typedef pair<int , int > pii;
typedef unsigned long long ull;

namespace FastIO
{
// char buf[1 << 20], *p1, *p2;
// #define getchar() (p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 1 << 20, stdin), p1 == p2) ? EOF : *p1++)
template<typename T> inline T read(T& x) {
x = 0;
int f = 1;
char ch;
while (!isdigit(ch = getchar())) if (ch == '-') f = -1;
while (isdigit(ch)) x = (x << 1) + (x << 3) + (ch ^ 48), ch = getchar();
x *= f;
return x;
}
template<typename T, typename... Args> inline void read(T& x, Args &...x_) {
read(x);
read(x_...);
return;
}
inline ll read() {
ll x;
read(x);
return x;
}
};
using namespace FastIO;

const int N = 1e5 + 10;

int n;
ll b[N], m[N];

inline void Input() {
read(n);
for(int i = 1; i <= n; i++) {
read(m[i], b[i]);
}
}

ll exgcd(ll a, ll b, ll&x, ll&y) {
if(b == 0) {
x = 1, y = 0;
return a;
}
ll d = exgcd(b, a % b, y, x);
y -= a / b * x; return d;
}

ll Mul(ll a, ll b, ll P) {
ll ans = 0;
while(b) {
if(b & 1) ans = (ans + a) % P;
a = (a + a) % P;
b >>= 1;
}
return ans;
}

inline void Work() {
ll ans = b[1], lm = m[1];
ll x, y, d, bg, c;
for(int i = 2; i <= n; i++) {
c = (b[i] - ans % m[i] + m[i]) % m[i];
d = exgcd(lm, m[i], x, y);
bg = m[i] / d;
if(c % d) {
printf("-1");
return ;
}
x = Mul(x, c / d, bg);
ans += x * lm, lm *= bg;
ans = (ans % lm + lm) % lm;
}
ans = (ans % lm + lm) % lm;
printf("%lld\n", ans);
}

int main() {
int T = 1;
while(T--) {
Input();
Work();
}
return 0;
}

算法实战:须绳缚身,沉潜勿用

查看相关算法标签喵!

参考资料(以下排名不分先后)