算法学习笔记(9): 中国剩余定理(CRT)以及其扩展(EXCRT)
扩展中国剩余定理
讲解扩展之前，我们先叙述一下普通的中国剩余定理
“China Remain Theory” 也叫做孙子定理
难得是以中国命名的定理
，谁敢不掌握
中国剩余定理
中国剩余定理通过一种非常精巧的构造求出了一个可行解
但是毕竟是构造，所以相对较复杂
对于上述同余方程组，首先需要满足 \(m_1, m_2, m_3, \dots, m_n\) 互质
令 \(m = \prod_{i=1}^{n} m_i, M_i = m / m_i,\ t_i\) 是线性同余方程 \(M_it_i \equiv 1 \pmod {m_i}\) 的一个解
那么上述方程组有整数解，为 \(x = \sum_{i=1}^n a_iM_it_i\)
证明：
由于所有 \(m_i\) 互质，所以 \(M_i = m/m_i\) 是除了 \(m_i\) 之外的所有模数的倍数，所以 \(\forall j \ne i, a_iM_it_i \equiv 0 \pmod {m_k}\)
所以代入 \(x = \sum_{i=1}^n a_iM_it_i\)，原方程组成立
中国剩余定理给出了方程的一个特解。通解可以表示为 \(x + km (k \in \Z)\)。如果需要求最小的非负整数解，只需要让 \(x\) 对于 \(m\) 取模就是。
参考代码
#include <iostream>
typedef long long ll;
ll m[11], a[11], Mm[11], y[11];
// 扩展欧几里得算法，由于求出线性同余方程 Mi * ti = 1 (mod p) 的一个特解
// 上述同余式可以转化为 Mi * ti + k * p = 1
// 所以跑一个 extgcd(Mi, p, ti, k) 即可 （k没有用的）
void extgcd(ll a, ll b, ll &x, ll &y) {
if (b == 0) {
x = 1, y = 0;
return;
}
extgcd(b, a % b, y, x);
y -= (a / b) * x;
}
int main() {
int n; scanf("%d", &n);
ll M(1);
// 这里读入，并求出m
for (int i = 0; i < n; i++) {
scanf("%lld%lld", m + i, a + i);
M *= m[i];
}
// 求出Mi
for (int i = 0; i < n; i++) {
Mm[i] = M / m[i];
}
// 求出每一个不定方程的特解
for (int i = 0; i < n; i++) {
ll x, tmp, mi = m[i];
extgcd(Mm[i], mi, x, tmp);
y[i] = (x + mi) % mi;
}
// 求和获得最小正整数解
ll x(0);
for (int i = 0; i < n; i++) {
x = (x + a[i] * Mm[i] * y[i]) % M;
}
printf("%lld\n", x);
return 0;
}
扩展中国剩余定理
还是考虑上述式子
此时，我们不保证 \(m_i\) 不一定两两互质，所以中国剩余定理不再适用
所以扩展中国剩余定理就开始发挥作用了
虽然说是扩展……但是两者思路很不一样
扩展中国剩余定理采用数学归纳法，或者说两两合并法
假如我们需要合并方程
我们将同余式改写
将两者相减，整理后可得
由于我们已知 \(m_1, m_2, a_1, a_2\)，也就是说我们只需要知道 \(k_1, k_2\) 其中任意一个，就可以求出这时的 \(x\) 是个什么东西了
所以，联想到了……扩展欧几里得算法
可以参考文章：算法学习笔记(1): 欧几里得算法及其扩展
我们令
那么化简一下有：\(k_1 m_1 - k_2 m_2 = c\)
所以我们求出 \(k_1' \frac {m_1}d + k_2' \frac {m_2}d = 1\)
那么换回去之后
由于 \(exgcd(a, b, x, y)\) 求出的是 \(xa + yb = \gcd(a, b)\) 的一组解……所以需要 \(\frac cd\)
令 \(l = lcm(m_1, m_2) = m_1 m_2 / \gcd(m_1, m_2)\)
那么此时 \(x \equiv k_1' m_1\frac cd + a_1 \pmod l\)
也就是可以得到 \(x\) 的一个解集 \(x \in X = \{kl + x|k \in \R\}\)
也就是说，我们将上述两个式子合并成了
那么考虑代码怎么写
typedef long long lint;
struct Equation {
lint a, p;
};
inline lint gcd(lint x, lint y) {
return y ? gcd(y, x % y) : x;
}
inline lint exgcd(lint x, lint y, lint &s, lint &t) {
if (!y) {
s = 1, t = 0;
return x;
}
lint r = exgcd(y, x % y, t, s);
t -= (x / y) * s;
return r;
}
bool merge(const Equation &a, const Equation &b, Equation &res) {
// 第一部分：预先求出需要的值
lint d = gcd(a.p, b.p), s, t;
lint lcm = a.p / d * b.p;
// 第二部分：求出上文中的 k1', k2'
lint c = (b.a - a.a) % lcm;
if (c < 0) c += lcm;
if (c % d) return false; // 判断有无解，返回false则无解
exgcd(a.p / d, b.p / d, s, t);
// 第三部分：求出上文中的 k1
s = s * (c / d) % lcm;
if (s < 0) s += lcm;
// 第四部分：求出一个特殊的解
lint x = (a.a + s * a.p % lcm) % lcm;
if (x < 0) x += lcm;
res = {x, lcm};
return true;
}
?：为什么在前三个部分放在了模 \(lcm\) 的情况下计算
其实最小的正整数解在 \(\mod lcm\) 的情况下是唯一的~
唯一性的证明可以参考： 阮行止 的博客
所以放在\(\mod lcm\) 的条件下不会对合并造成影响，并且还保证了求出的解是最小非负整数解。
?: 为什么在第一部分需要提前算出 \(\gcd(m_1, m_2)\)
其实是不用的……只需要
data s, t; data d = exgcd(m1, m2, s, t);也可以算出正确答案……
?: 判断有解是如何判断的
有无解实际上就是拓展欧几里得算法有无解
而判断拓欧算法有无解，也就通过贝祖定理验证
也就是说，对于不定方程
ax + by = c如果
gcd(x, y) | c则有解，否则无解
NOTICE：在模板题上，由于可能会溢出……至于什么地方需要用龟速乘，请读者自行揣度
关于龟速乘：算法学习笔记(2): 逆元及其应用我提过一点
更详细的参考 快速乘总结 - 一只不咕鸟
后记
后来才发现，原来写的都是在胡扯……所以重新写了 exCRT 部分。好多都被网上的一些文章误导了，他们在 exCRT 前三个部分很多都是在\(\mod p_2\) 意义下进行的。虽然这样无可厚非，但是在判断有无解时会出现问题，所以……
