Q:今有物不知其数,三三数之剩二,五五数之剩三,七七数之剩二,问物几何?
A:三人同行七十稀,五树梅花廿一枝,七子团圆正半月,除百零五便得知。
可能高二时候我也没想到,这个当时我特别喜欢的算法书上的定理,再次见到是在数学专业课高代上,而不是计算机算法课上。
中国剩余定理又称孙子定理,中国古代求解一次同余式组(见韩信点兵和曹冲养猪)的方法,最早可见于南北朝的数学著作《孙子算经》。
本文给出中国剩余定理的证明及程序实现。
需要注意的是,由于整数环和多项式环性质极其相似,该定理在多项式环上也有相似结论。
Bezout定理
[latex]\large \forall p,q\in \mathbb{N} ,\exist u,v\in \mathbb{N} \; s.t. \; up+vq=(p,q)[/latex]
其中[latex](p,q)[/latex]为p,q的最大公约数。
Bezout定理通过Euclidean算法易证。
Bezout定理推论
若[latex]p,q[/latex]互质,则有[latex]up\equiv 1\pmod{q}[/latex]
中国剩余定理(Chinese Remainder Theorem)
已知[latex]m_{i}(i\in \mathbb{N} _{+})[/latex]满足[latex]\forall i,j\in \mathbb{N} _{+} ,i\not= j, (m_{i}, m_{j} )=1[/latex],
则[latex]\exist x \in \mathbb{N}, x \equiv bi \pmod{m_{i}}[/latex]。
具体地说,一定有[latex]x[/latex]:
[latex]\large x \equiv \sum b_{i} M_{i} N_{i} \pmod{M}[/latex],
其中[latex]\large M = \prod m_{i} , M_{i} = \prod _{k \not= i} m_{k}, N_{i} M_{i} \equiv 1 \pmod {m_{i}}[/latex]。
[latex]N_{i}[/latex]的存在性
由于[latex]m_{i}[/latex]两两互质,而[latex]M_{i} = \prod _{k \not= i} m_{k}[/latex],总是有[latex](m_{i} , M_{i} )=1[/latex],从而由Bezout定理:
[latex]\large \exist u_{i},N_{i}\in \mathbb{N} \; s.t. \; u_{i} m_{i} + N_{i} M_{i} = 1 \implies N_{i} M_{i} \equiv 1 \pmod{m_{i}}[/latex]
解的合理性
注意到[latex]N_{i} M_{i} \equiv 1 \pmod{m_{i}}[/latex],
而[latex]N_{i} M_{i} \equiv 0 \pmod{m_{j}} (i \not= j)[/latex],
从而满足题设方程,而[latex]\forall k \in \mathbb{Z} , x + k M_{i} \equiv x \pmod{m_{i}}[/latex],故模[latex]M[/latex]意义下解合理。
思考
中国剩余定理的构造和Langrange插值公式极其相似,按罗栗老师的话:都有类似开关的东西在控制整个式子。精妙所在。
程序实现
洛谷 P1495 【模板】中国剩余定理(CRT)/曹冲养猪
其实整个算法关键在于求[latex]M_{i}[/latex]的乘法逆元,用扩展Euclidean算法就可以了。然后再按解的构造累加起来,模一下M得到最小正整数解,为了防止溢出我每一步都模了,但是题解上看似乎最终再模也没有问题。
第一次提交,听取[latex]\color{red}{\rm WA}[/latex]声一片
原因写在代码注释里了,修改后完美[latex]\color{green}{\rm AC}[/latex]
#include <cstdio>
#include <cctype>
typedef long long int int64;
int64 n, a[20], b[20], Mi[20], M = 1, X;
inline int getint()
{
int x = 0;
char ch = 0;
while (!isdigit(ch)) ch = getchar();
while (isdigit(ch)) x = (x << 1) + (x << 3) + (ch ^ 0x30), ch = getchar();
return x;
}
int64 exgcd(int64 a, int64 b, int64 &x, int64 &y)
//扩展Euclidean算法求解Bezout方程
{
if (b == 0) {
x = 1; y = 0;
return a;
}
int64 r = exgcd(b, a % b, y, x);
y -= a / b * x;
return r;
}
int main()
{
n = getint();
for (int i = 1; i <= n; i++) {
a[i] = getint();
b[i] = getint();
Mi[i] = 1;
M *= a[i];
}
for (int i = 1; i <= n; i++)
Mi[i] = M / a[i];
for (int i = 1; i <= n; i++) {
int64 Ni, y;
exgcd(Mi[i], a[i], Ni, y); //求解同余方程Ni * M[i] = 1(mod a[i])
X = (X + b[i] * Mi[i] * (Ni % a[i] + a[i])) % M;
//Ni % a[i] + a[i] 防止Ni为负数,注意这里Ni为模a[i]意义下的乘法逆元,不能将a[i]写为M
//写成X = ((X + b[i] * Mi[i] * Ni) % M + M) % M应该也行;
}
printf("%lld\n", X);
return 0;
}