【算法笔记】中国剩余定理(CRT)与扩展中国剩余定理(exCRT)
- 本文总计约 8000 字,阅读大约需要 40 分钟。
- 警告!警告!警告!本文有大量的 \(\LaTeX\) 公式渲染,可能会导致加载异常缓慢!
三人同行七十稀,五树梅花廿一支。
七子团圆月正半,除百零五便可知。— —明代数学家程大位对『物不知数』问题解法给出的歌诀
前言
我觉得可以不要前言 QwQ。
问题引入
在我国古代的数学巨著《孙子算经》中,有这样一个名叫『物不知数』的问题:
今有物不知数,三三之数剩二,五五之数剩三,七七之数剩二,问物几何?
用现代的话翻译过来,就是求下面的关于 \(x\) 的同余方程组的解:
当然,我们可以计算出来最小的 \(x=23\),那么,对于任意一个下面这样的关于 \(x\) 的同余方程:
它一定有解吗?如果有解,能不能给出一个高效的算法,求出来满足该方程的最小整数解呢?
中国剩余定理(CRT)
像上面这样,解决多个线性同余方程组的问题,被称为中国剩余问题。
在先秦时期,就有了类如『隔墙算』,『剪管术』,『韩信点兵』一类的问题,具体内容都是关于求上面这样线性同余方程组的数学问题。
然而,这类问题的文献记载,最早出现在南北朝时期的数学著作《孙子算经》中,作者将其称为『物不知数』问题,也就是题目引入中的那个问题。
后来,宋朝数学家秦九韶在他的数学著作《大衍章》中给出了比较系统的解答,并且给出了一个重要的定理,即所谓的中国剩余定理(Chinese Remainder Theorem,CRT,又称孙子定理):
设线性同余方程组:
如果 \(m_1,m_2\cdots,m_n\) 两两互质,则该方程组一定有解,且在模 \(M=\prod m_i\) 意义下解唯一。
剩余定理对解的构造
『好的,我们已经知道了这个问题一定有解了,但是我们怎么样才能得到这个解呢?』
这个时候,我们就要请出我国古代的数学家,程大位了。他将『物不知数』问题的解法变成了一首四句歌诀:
三人同行七十稀,五树梅花廿一支。七子团圆月正半,除百零五便可知。
用现在的话来翻译一下,就是用 \(x\) 除以 \(3\) 所得的余数 \(2\),乘以一个 \(70\)(三人同行七十稀);
加上 \(x\) 除以 \(5\) 所得的余数 \(3\),乘上 \(21\)(即廿一,五树梅花廿一支);
再加上 \(x\) 除以 \(7\) 所得的余数 \(2\),乘上 \(15\)(半个月即十五天,七子团圆月正半);
得到的结果为 \(2\times 70+3\times 21+2\times 15=233\),最后对 \(105\) 取模,得到 \(308\bmod 105=23\),就是最终的答案了(除百零五便可知)。
这四句歌谣的确把『物不知数』问题算清楚了,但是,里面还有一些问题没有解决。
『对啊,\(105=3\times 5\times 7\) 尚且不论,那些 \(70,21,15\) 是怎么得出来的呢?难道是通过瞎蒙蒙出来的?』
我们就要考察一下这些数字的含义了。
对于一个线性同余方程组:
既然 \(x\) 在模 \(M\) 意义下有唯一解,我们就构造一个和 \(M\) 相关的解吧!
定义 \(M=\prod_{i=1}^n{m_i}\),\(M_i=\dfrac{M}{m_i}\),那么我们构造:
其中,\(t_i\) 为 \(M_i\) 在模 \(m_i\) 意义下的逆元。
我们下面可以证明,\(x\) 为该同余方程的一个解:
对于第 \(i\) 个方程 \(x\equiv a_i\pmod {m_i}\),我们代入上式,得
因为 \(M_1,M_2,M_{i-1}\) 和 \(M_{i+1},M_{i+2},\cdots,M_n\) 里都含有因子 \(m_i\),所以:
那么原式只剩下一个 \(x\equiv a_iM_it_i\pmod{m_i}\) 了,我们又知道, \(t_i\) 是 \(M_i\) 模 \(m_i\) 意义下的逆元,有 \(M_it_i\equiv 1\pmod {m_i}\),所以 \(x\equiv a_i\pmod {m_i}\) 成立。
综上,\(x=(\sum_{i=1}^n a_iM_it_i)\bmod M\)。
注意到 \(m_i\) 不一定是质数,所以我们求 \(M_i\) 模 \(m_i\) 意义下的逆元需要用扩欧来完成。
代码
代码如下:
#include <iostream>
using namespace std;
const int MAXN = 100001;
typedef long long ll;
int n;
ll a[MAXN], m[MAXN], M = 1;
ll ans;
void exgcd(ll a, ll b, ll &x, ll &y) { //扩欧,不解释
if(!b) {
x = 1, y = 0; return ;
}
exgcd(b, a % b, y, x);
y -= a / b * x;
}
int main() {
scanf("%d", &n);
for(int i = 1; i <= n; ++i) {
scanf("%lld%lld", &m[i], &a[i]); //代表方程 x = a[i] (mod m[i])
M *= m[i]; //求出 M
}
for(int i = 1; i <= n; ++i) {
ll xi, yi;
exgcd(M / m[i], m[i], xi, yi); //扩欧求逆元
xi = (xi % m[i] + m[i]) % m[i]; //要是正的
ans = (ans + a[i] * (M / m[i]) * xi) % M; //一边加一边取模
}
printf("%lld", ans);
return 0;
}
//by CaO
扩展中国剩余定理(exCRT)
CRT 可以解决当 \(m_1,m_2,\cdots,m_n\) 两两互质时的情况,但是如果 \(m_1,m_2,\cdots,m_n\) 不能做到两两互质,那么 CRT 的算法就束手无策了— —在求逆元的时候会出问题。
于是,exCRT 应运而生,它的思想不再是简单地构造一个解,而是通过合并两个同余方程的方式解答剩余问题的。
怎么合并呢
我们考虑下面的两个同余方程:
我们可以将其改写一下,就变成了:\(x=k_1a_1+m_1=k_2a_2+m_2\),其中,\(k_1,k_2\in \mathbb{N}\)。
做个变形,我们就得到:\(a_1k_1-a_2k_2=m_2-m_1\),这是一个关于 \(k_1,k_2\) 的不定方程,由 Bézout 定理,\(k_1,k_2\) 存在,当且仅当 \(\gcd(a_1,a_2)\mid (m_2-m_1)\)。
这样,我们就有了:
它也可以表示为下式:
我们设 \(\left(\dfrac{a_1}{\gcd(a_1,a_2)}\right)^{-1}\equiv K\quad\left(\bmod {\dfrac{a_2}{\gcd(a_1,a_2)}}\right)\),则上式可转化为:
即:
由 \(x=k_1a_1+m_1\),得:
注意到 \(\dfrac{a_1a_2}{\gcd(a_1,a_2)}=\operatorname{lcm}(a_1,a_2)\),上式写成同余式即:
像这样,我们就把两个同余方程合并成了一个,如果是 \(n\) 个的话,像这样一个一个合并,就可以解决了!
代码
代码如下:
#include <bits/stdc++.h>
using namespace std;
const int MAXN = 300001;
typedef __int128 ll;
int n;
ll a[MAXN], b[MAXN], xi, yi;
ll mul(ll A, ll B) {
ll res = 0;
while(B) {
if(B & 1) res = res + A;
B >>= 1, A = A + A;
}
return res;
}
ll exgcd(ll a, ll b, ll &x, ll &y) {
if(!b) {
x = 1, y = 0; return a;
}
ll ans = exgcd(b, a % b, y, x);
y -= a / b * x;
return ans;
}
ll inv(ll a, ll b) {
ll x, k, res = exgcd(a, b, x, k);
ll t = b / res;
return (x % t + t) % t;
}
int main(void) {
scanf("%d", &n);
for(int i = 1; i <= n; ++i)
scanf("%lld%lld", &a[i], &b[i]);
ll r = exgcd(a[1], a[2], xi, yi);
ll A = mul(a[1], a[2] / r);
ll B = (inv(a[1] / r, a[2] / r) * (b[2] - b[1]) / r % (a[2] / r) + a[2] / r) % (a[2] / r) * a[1] + b[1]; //合并两个同余方程
for(int i = 3; i <= n; ++i) {
r = exgcd(A, a[i], xi, yi);
B = (inv(A / r, a[i] / r) * (b[i] - B) / r % (a[i] / r) + a[i] / r) % (a[i] / r) * A + B;
A = mul(A, a[i] / r); //继续合并同余方程
}
printf("%lld", B);
return 0;
}
//by CaO
应用
CRT 在多项式领域有一个非常重要的应用,就是辅助完成任意模数 NTT(MTT)。
传统的 NTT 可以解决多项式乘法问题,但是如果给出的模数不能写成 \(2^ka+1\) 的形式,我们就不一定能找到模数的原根,也不一定能保证原根有成为单位根的价值。
我们一般可以选择三个模数 \(998244353\),\(1004535809\),\(469762049\) 三个质数(都可以写成 \(2^ka+1\) 的形式,且 \(3\) 是它们的原根)分别作 NTT,并且使用 CRT 将答案合并起来。
这种做法,我们成为任意模数 NTT,是一种非常有用的多项式算法。我们会在后续的学习中讲到。
例题
本题目列表会持续更新。