【演算法筆記】中國剩餘定理（CRT）與擴充套件中國剩餘定理（exCRT）
- 本文總計約 8000 字，閱讀大約需要 40 分鐘。
- 警告！警告！警告！本文有大量的 \(\LaTeX\) 公式渲染，可能會導致載入異常緩慢！
三人同行七十稀，五樹梅花廿一支。
七子團圓月正半，除百零五便可知。— —明代數學家程大位對『物不知數』問題解法給出的歌訣
前言
我覺得可以不要前言 QwQ。
問題引入
在我國古代的數學鉅著《孫子算經》中，有這樣一個名叫『物不知數』的問題：
今有物不知數，三三之數剩二，五五之數剩三，七七之數剩二，問物幾何？
用現代的話翻譯過來，就是求下面的關於 \(x\) 的同餘方程組的解：
\[\begin{cases} x \equiv 2\pmod 3,\\ x \equiv 3\pmod 5,\\ x \equiv 2\pmod 7. \end{cases} \]當然，我們可以計算出來最小的 \(x=23\)
它一定有解嗎？如果有解，能不能給出一個高效的演算法，求出來滿足該方程的最小整數解呢？
中國剩餘定理（CRT）
像上面這樣，解決多個線性同餘方程組的問題，被稱為中國剩餘問題。
在先秦時期，就有了類如『隔牆算』，『剪管術』，『韓信點兵』一類的問題，具體內容都是關於求上面這樣線性同餘方程組的數學問題。
然而，這類問題的文獻記載，最早出現在南北朝時期的數學著作《孫子算經》中，作者將其稱為『物不知數』問題，也就是題目引入中的那個問題。
後來，宋朝數學家秦九韶在他的數學著作《大衍章》中給出了比較系統的解答，並且給出了一個重要的定理，即所謂的中國剩餘定理（Chinese Remainder Theorem，CRT，又稱孫子定理）：
設線性同餘方程組：
\[\begin{cases} x\equiv a_1 \pmod {m_1}, & (1)\\ x\equiv a_2 \pmod {m_2}, & (2)\\ \cdots\\ x\equiv a_n \pmod {m_n}. & (n)\\ \end{cases} \]如果 \(m_1,m_2\cdots,m_n\)
剩餘定理對解的構造
『好的，我們已經知道了這個問題一定有解了，但是我們怎麼樣才能得到這個解呢？』
這個時候，我們就要請出我國古代的數學家，程大位了。他將『物不知數』問題的解法變成了一首四句歌訣：
三人同行七十稀，五樹梅花廿一支。七子團圓月正半，除百零五便可知。
用現在的話來翻譯一下，就是用 \(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\) 是怎麼得出來的呢？難道是通過瞎濛濛出來的？』
我們就要考察一下這些數字的含義了。
對於一個線性同餘方程組：
\[\begin{cases} x \equiv a_1 \pmod {m_1},\\ x \equiv a_2 \pmod {m_2},\\ \cdots\\ x \equiv a_n \pmod {m_n}.\\ \end{cases} \]既然 \(x\) 在模 \(M\) 意義下有唯一解，我們就構造一個和 \(M\) 相關的解吧！
定義 \(M=\prod_{i=1}^n{m_i}\)，\(M_i=\dfrac{M}{m_i}\)，那麼我們構造：
\[\begin{aligned} x &= (a_1M_1t_1+a_2M_2t_2+\cdots+a_nM_nt_n)\bmod M\\ &= \left(\sum_{i=1}^n{a_iM_it_i}\right)\bmod M. \end{aligned} \]其中，\(t_i\) 為 \(M_i\) 在模 \(m_i\) 意義下的逆元。
我們下面可以證明，\(x\) 為該同餘方程的一個解：
對於第 \(i\) 個方程 \(x\equiv a_i\pmod {m_i}\)，我們代入上式，得
\[a_1M_1t_1+a_2M_2t_2+\cdots+a_nM_nt_n \equiv a_i \pmod{m_i}. \]因為 \(M_1,M_2,M_{i-1}\) 和 \(M_{i+1},M_{i+2},\cdots,M_n\) 裡都含有因子 \(m_i\)，所以：
\[\sum_{j\neq i}a_jM_jt_j\equiv 0\pmod{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 應運而生，它的思想不再是簡單地構造一個解，而是通過合併兩個同餘方程的方式解答剩餘問題的。
怎麼合併呢
我們考慮下面的兩個同餘方程：
\[\begin{cases} x\equiv a_1 \pmod{m_1},\\ x\equiv a_2 \pmod{m_2}. \end{cases} \]我們可以將其改寫一下，就變成了：\(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)\)。
這樣，我們就有了：
\[\dfrac{a_1k_1}{\gcd(a_1,a_2)}=\dfrac{m_2-m_1}{\gcd(a_1,a_2)}+\dfrac{a_2k_2}{\gcd(a_1,a_2)}. \]它也可以表示為下式：
\[\dfrac{a_1k_1}{\gcd(a_1,a_2)}\equiv \dfrac{m_1-m_2}{\gcd(a_1,a_2)} \quad\left(\bmod{\dfrac{a_2}{\gcd(a_1,a_2)}}\right). \]我們設 \(\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)\)，則上式可轉化為：
\[k_1\equiv K\dfrac{m_1-m_2}{\gcd(a_1,a_2)}\quad\left(\bmod{\dfrac{a_2}{\gcd(a_1,a_2)}}\right). \]即：
\[k_1=K\dfrac{m_1-m_2}{\gcd(a_1,a_2)}+p\dfrac{a_2}{\gcd(a_1,a_2)}. \]由 \(x=k_1a_1+m_1\)，得：
\[x_1=K\dfrac{m_1-m_2}{\gcd(a_1,a_2)}+Kp\dfrac{a_1a_2}{\gcd(a_1,a_2)}+m_1. \]注意到 \(\dfrac{a_1a_2}{\gcd(a_1,a_2)}=\operatorname{lcm}(a_1,a_2)\)，上式寫成同餘式即：
\[x\equiv K\dfrac{m_1-m_2}{\gcd(a_1,a_2)}+m_1\pmod{\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，是一種非常有用的多項式演算法。我們會在後續的學習中講到。
例題
本題目列表會持續更新。