淺談中國剩餘定理（CRT）及其擴充套件
中國剩餘定理（CRT）是中國古代求解一次同餘式組的方法，是數論中一個重要的定理。《孫子算經》中首次提到了同餘方程組問題以及解法，因此在數學文獻中也會將中國剩餘定理稱為孫子定理。
——摘自度娘百科
引入
有物不知其數，三三數之剩二，五五數之剩三，七七數之剩二。問物幾何？
此題源自《孫子算經》卷下第二十六題，大體意思是：有一個數，除以三餘數為二，除以五餘數為三，除以七餘數為二，問這個數是多少。對於這個問題，我們可以依據一個歌謠解決：
三人同行七十稀，五樹梅花廿一支，七子團圓正半月，除百零五便得知。
這句歌謠的意思是：把除以三的餘數乘上70，除以五的餘數乘上21，除以七的餘數乘上15，他們的和除以105的餘數就是答案。代入一開始的問題就可以得到，答案是\((2 \times 70 + 3 \times 21 + 2 \times 15) \mod 105\ =\ 23\)
簡單的分析
那麼為什麼要乘上70、21和15這三個奇怪的數呢？我們通過觀察發現：
\[\left\{ \begin{array}{lr} 70 \equiv 1 \pmod 3, &\\ 70 \equiv 0 \pmod 5, &\\ 70 \equiv 0 \pmod 7. \end{array} \right. \left\{ \begin{array}{lr} 21 \equiv 0 \pmod 3, &\\ 21 \equiv 1 \pmod 5, &\\ 21 \equiv 0 \pmod 7. \end{array} \right. \left\{ \begin{array}{lr} 15 \equiv 0 \pmod 3, &\\ 15 \equiv 0 \pmod 5, &\\ 15 \equiv 1 \pmod 7. \end{array} \right. \]
可以發現，每加上一個數，就是為以除這個數餘數為一的數為除數的餘數貢獻1，而對其他的除數沒有影響！沒錯這就是個繞口令 又因為\(\operatorname{lcm}(3,5,7)=105\)，所以除以105的餘數就是這個同餘方程的最小整數解。
CRT入門
從簡單到複雜
CRT能夠解決的問題當然不止模數為3、5、7的問題。實際上，CRT可以求出如下同餘方程的解：
\[\left\{ \begin{array}{lr} x \equiv q_1 \pmod {p_1}, &\\ x \equiv q_2 \pmod {p_2}, &\\ x \equiv q_3 \pmod {p_3}, & \left(\gcd\{p_1,p_2,p_3\dots,p_n\} = 1\right)\\ \dots &\\ x \equiv q_n \pmod {p_n}. \end{array} \right. \]
我們令\(\forall i \in [1,n]\quad m_i=\dfrac{\operatorname{lcm}(p_1,p_2\dots,p_n)}{p_i}=\dfrac{\prod_{j=1}^n p_j}{p_i}\)，則答案即為\(\sum_{i=1}^n q_i\times m_i \times m_i^{-1}\)，其中\(m_i^{-1}\)為\(m_i\)在模\(p_i\)意義下的逆元。
證明？
即得易見平凡，仿照上例顯然。留作習題答案略，讀者自證不難。反之亦然同理，推論自然成立，略去過程QED，由上可知證畢。
可以知道，\(\forall i,j\in [1,n]\)
-
當\(i \ne j\)時有$ q_j \times m_j \times m_j^{-1} \equiv 0 \pmod {p_i}$
-
當\(i = j\)時有$ q_j \times m_j \times m_j^{-1} \equiv q_j \pmod {p_i}$
因此，\(\forall j \in [1,n] \sum_{i=1}^n q_i\times m_i \times m_i^{-1} \equiv q_j \pmod {p_j}\)，因此定理成立。
關於逆元的那些事
乘法逆元，主要用來求模\(p\)意義下\(\dfrac{a}{b}\)的值。（\(p\)一般是質數）
若\(a \cdot x \equiv 1 \pmod b\)，其中\(\gcd(a,b)=1\)，則稱\(x\)為\(a\)的乘法逆元，用\(a^{-1}\)表示。
如何求逆元？
-
擴充套件歐幾里得演算法（
exgcd）乘法逆元的本質是\(a \cdot x \equiv 1 \pmod b\)，我們就可以把這個柿子轉化成這個亞子：\(ax+by=1\)，其中\(x,y \in \mathbb{Z}\)。而這個不定方程就可以用
exgcd求解。關於exgcd的詳細解釋可以看 這道題。程式碼：
int exgcd(int a,int b,int& x,int& y) { // exgcd模板 if(!b) { x = 1,y = 0; return a; } int x1,y1,ans = exgcd(b,a % b,x1,y1); x = y1,y = x1 - a / b * y1; return ans; } int inv(int a,int b) { // 求逆元 int x,y; exgcd(a,b,x,y); //exgcd return (x % b + b) % b; // 負數處理 } -
費馬小定理&快速冪
費馬小定理：\(\forall a \in \mathbb{Z},p \in prime,gcd(a,p) = 1\)有\(a^{p - 1}\equiv 1 \pmod p\)
我們注意到這個柿子的右邊正好是個1，於是就可以把這個柿子進行億些小小的變換：
\[\left\{ \begin{array}{lr} ax \equiv 1 \pmod p & \\ a^{p-1} = a \cdot a^{p-2} \equiv 1 \pmod p \end{array} \right. \]
於是我們可以得到:
\[x \equiv a^{p-2} \pmod p \]
而這個柿子可以用快速冪快速求解。
程式碼:
int quick_pow(int a,int b,int p) { if(!b) return 1; if(b == 1) return a % p; int rep = quick_pow(a,b >> 1,p); rep = rep * rep % p; if(b & 1) return rep * a % p; return rep; } int inv(int a,int p) { return quick_pow(a % p,p - 2,p); }
程式碼實現
別裝了，我知道你們都在等著這個（大霧
// n表示同餘方程組數，p陣列表示模數，q陣列表示餘數
int crt() {
int M = 1,ans = 0; // 乘積M，結果ans
for(int i = 1;i <= n;i ++) M *= p[i];
for(int i = 1;i <= n;i ++) { // 迴圈列舉1~n
int m = M / p[i];
int x = inv(m); // 求逆元
ans += q[i] * m * x; //求和
}
return ans % M;
}
來個栗子
右轉美食街栗子店炒鍋
很模板的一道CRT，建立\(a_i\)個豬圈有\(b_i\)頭豬無家可歸就是指\(x \equiv b_i \pmod{a_i}\)，又因為\(\forall i,j\in[1,n] \gcd(a_i,a_j)=1\)，所以這道題可以用CRT來做。
程式碼：
#include <cstdio>
typedef long long ll; // 一定要記得開long long！
ll n,a[15],b[15];
ll exgcd(ll a,ll b,ll& x,ll& y) { // exgcd模板
if(!b) {
x = 1,y = 0;
return a;
}
ll x1,y1,ans = exgcd(b,a % b,x1,y1);
x = y1,y = x1 - a / b * y1;
return ans;
}
ll inv(ll a,ll b) { // exgcd求逆元
ll x,y;
exgcd(a,b,x,y);
return (x % b + b) % b;
}
ll crt() { // crt 模板
ll M = 1,ans = 0;
for(int i = 1;i <= n;i ++) M *= 1LL * a[i];
for(int i = 1;i <= n;i ++) {
ll m = M / a[i],x = inv(m,a[i]);
ans += 1LL * b[i] * m * x;
}
return ans % M;
}
int main() {
scanf("%d",&n);
for(int i = 1;i <= n;i ++) scanf("%d%d",a + i,b + i);
printf("%lld\n",crt());
return 0;
}
還是一道模板題，不過因為模數已知，我們可以用\(\mathcal{O}(1)\)的方式求解。
根據題意我們可以列出這樣的方程式：
\[\left\{ \begin{array}{lr} x_1 \equiv 1 \pmod {23}, &\\ x_1 \equiv 0 \pmod {28}, &\\ x_1 \equiv 0 \pmod {33}. \end{array} \right. \left\{ \begin{array}{lr} x_2 \equiv 0 \pmod {23}, &\\ x_2 \equiv 1 \pmod {28}, &\\ x_2 \equiv 0 \pmod {33}. \end{array} \right. \left\{ \begin{array}{lr} x_3 \equiv 0 \pmod {23}, &\\ x_3 \equiv 0 \pmod {28}, &\\ x_3 \equiv 1 \pmod {33}. \end{array} \right. \]
解得：
\[\left\{ \begin{array}{lr} x_1 = 5544 & \\ x_2 = 14421 & \\ x_3 = 1288 \end{array} \right. \]
根據之前的求解柿子可以得出:
\[\operatorname{lcm}\{21,28,33\} = 21252 \\ ans \equiv 5544p + 14421e + 1288i \pmod{21252} \]
所以輸出的答案就是\((5544p + 14421e + 1288i - d + 21251) \mod 21252 + 1\)，其中括號內加上21251是為了避免負數，括號外加1是特判0的情況，此時應輸出21252。
CRT拓展
擴充套件中國剩餘定理（excrt）
你學會了中國剩餘定理，非常高興，於是去向\(\color{black}l\color{red}ndjy\)大佬吹牛皮。這時\(\color{black}l\color{red}ndjy\)大佬向你拋來了一個同餘方程組：
\[\left\{ \begin{array}{lr} x \equiv q_1 \pmod {p_1}, &\\ x \equiv q_2 \pmod {p_2}, &\\ x \equiv q_3 \pmod {p_3}, &\\ \dots &\\ x \equiv q_n \pmod {p_n}. \end{array} \right. \]
（\(q_1,q_2,q_3\dots,q_n\)不兩兩互質）
這時中國剩餘定理不管用了，因為當\(i,j\in [1,n],i \ne j\)時$ q_j \times m_j \times m_j^{-1} \mod {p_i}$不一定為0！
怎麼辦呢？多半是廢了
歡迎我們的嘉賓——擴充套件中國剩餘定理（excrt）！別問，問就是作者中二病又發作了
我們考慮分開求解然後合併，假設求解到了第\(i\)個方程，前\(i-1\)個方程求解的結果是\(ans\)，那麼可以得到：
\[\left\{ \begin{array}{lr} x \equiv ans \pmod {M}, &\\ x \equiv q_i \pmod {p_i}. \end{array} \right. \\ \left(M = \operatorname{lcm}\{p_1,p_2,p_3\dots,p_n\}\right) \]
為了簡便我們設\(a_1=ans,b_1=M,a_2=q_i,b_2=p_i\)，那麼：
\[x=a_1+b_1k_1=a_2+b_2k_2 \\ \therefore b_1k_1-b_2k_2=a_2-a_1 \]
我們利用擴充套件歐幾里得求解\(b_1x'-b_2y'=\gcd(b_1,b_2)\)，那麼就可以化簡得到：
\[k_1=x'\times\dfrac{a_2-a_1}{\gcd(b_1,b_2)}\\ \therefore x=b_1k_1+a_1 \]
這部分推柿子可能有些難理解，建議自己手推一遍。
例題：洛谷P4777 【模板】擴充套件中國剩餘定理（EXCRT）
程式碼：
#include <cstdio>
typedef long long ll;
ll n,ans,p[100005],q[100005];
inline ll read() { // 快讀
#define gc c = getchar()
ll d = 0;
int f = 0,gc;
while(c < 48 || c > 57) f |= (c == '-'),gc;
while(c > 47 && c < 58) d = (d << 1LL) + (d << 3LL) + (c ^ 48),gc;
#undef gc
return f ? -d : d;
}
ll exgcd(ll a,ll b,ll& x1,ll& y1) { // exgcd模板
if(!b) {
x1 = 1,y1 = 0;
return a;
}
ll x2,y2,ans = exgcd(b,a % b,x2,y2);
x1 = y2;
y1 = x2 - a / b * y2;
return ans;
}
ll pmul(ll a,ll b,ll mod) { // 快速乘 不用的話會爆long long
ll rep = 0;
while(b) {
if(b & 1) rep = (rep + a) % mod;
a = (a << 1) % mod,b >>= 1;
}
return rep;
}
bool excrt(ll& ans) { // 返回值表示成功與否
ll M = p[1];
ans = q[1];
for(int i = 2;i <= n;i ++) {
ll k1,x,y;
ll g = exgcd(M,p[i],x,y); // 擴歐
ll c = (q[i] - ans % p[i] + p[i]) % p[i];
if(c % g != 0) return false; // 不能整除，無解
k1 = pmul(x,c / g,p[i] / g); // 見上面的公式 這裡加了一些取模操作
ans = M * k1 + ans;
M = M / g * p[i];
ans = (ans % M + M) % M;
}
ans = (ans % M + M) % M;
return true;
}
int main() {
n = read();
for(int i = 1;i <= n;i ++) p[i] = read(),q[i] = read();
excrt(ans);
printf("%lld\n",ans);
return 0;
}
擴充套件盧卡斯定理（exlucas）
你又學會了擴充套件中國剩餘定理，非常高興，於是又去向\(\color{black}l\color{red}ndjy\)大佬吹牛皮。這時\(\color{black}l\color{red}ndjy\)大佬向你拋來了一個柿子：
\[C_n^m \mod p\\ (p \notin prime) \]
求這個柿子的值。
這個問題可以用到擴充套件盧卡斯定理（exlucas）。這個演算法的思路大概是這樣的：先把模數p分解質因數：
\[p=\alpha_1^{\beta_1}\alpha_2^{\beta_2}\alpha_3^{\beta_3}\dots\alpha_n^{\beta_n} \]
然後再列出一個同餘方程組：
\[\left\{ \begin{array}{lr} C_n^m \equiv q_1 \pmod{\alpha_1^{\beta_1}},& \\ C_n^m \equiv q_2 \pmod{\alpha_2^{\beta_2}},& \\ C_n^m \equiv q_3 \pmod{\alpha_3^{\beta_3}},& \\ \dots & \\ C_n^m \equiv q_n \pmod{\alpha_n^{\beta_n}}. \end{array} \right. \]
然後我們就可以對形如\(C_n^m \equiv q \pmod{p^k} \quad (p \in prime)\)的同餘式分別求解，然後用CRT合併即可。求解單個方程式的過程比較繁瑣，這裡就不再贅述。