[笔记]中国剩余定理(CRT) & 扩展中国剩余定理(exCRT)
中国剩余定理(CRT)
对于线性同余方程组:
其中\(b_i\)两两互质。
下文中,令\(N=\prod\limits_{i=1}^n b_i\)。
考虑如何求解其中一个特解\(x\)。
根据同余运算的可加性,我们可以将原方程拆成\(n\)个:
则\(x=\sum_{i=1}^n x_i\)。
对于第一个方程组,我们令\(x_1=y_1\times a_1\),则有:
上面的转化过程是CRT的核心思想。
由后\(n-1\)个式子可知\(y_1\equiv 0\pmod{\text{lcm}_{i=2}^n\ b_i}\),即\(y_1\equiv 0\pmod{\frac{N}{b_1}}\),令\(c_i=\frac{N}{b_i}\),则\(y_1=c_1\times k_1(k_1\in \mathbb{Z})\)。
由第\(1\)个式子可知\(c_1\times k_1\equiv 1\pmod{b_1}\),于是\(k_1\equiv c_1^{-1}\pmod{b_1}\),即\(k_1\)是\(c_1\)在模\(b_1\)意义下的乘法逆元。
类比可得:
至此完成\(x\)的求解。
再考虑如何用特解\(x\)表示所有的通解。
考虑\(x,y\)同为该方程组的解,则根据同余的传递性,有:
也即\(x\equiv y\pmod{\text{lcm}_{i=1}^n b_i}\),即\(x\equiv y\pmod N\)。
所以通解可以表示为\(kN+x(k\in \mathbb{Z})\)。
这也告诉我们,对于\(k\in \mathbb{Z}\),该方程组的解在\([kN,(k+1)N)\)上是唯一的。
由于\(b_i\)不一定是质数,所以乘法逆元采用exGCD。
点击查看代码
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=15;
int n,fac=1,a[N],b[N];
__int128 ans;
int exgcd(int a,int b,int& x,int& y){
int d;
if(b==0) x=1,y=0,d=a;
else d=exgcd(b,a%b,y,x),y-=a/b*x;
return d;
}
signed main(){
ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
cin>>n;
for(int i=1;i<=n;i++){
cin>>b[i]>>a[i];
fac*=b[i];
}
for(int i=1;i<=n;i++){
int c=fac/b[i],k,t;
exgcd(c,b[i],k,t);
(ans+=(__int128)c*k*a[i])%=fac;
}
cout<<(long long)((ans+fac)%fac)<<"\n";
return 0;
}
扩展中国剩余定理(exCRT)
CRT不能解决\(b_i\)不两两互质的线性同余方程组。
这是因为模\(p\)意义下\(a\)的逆元存在\(\iff \gcd(a,p)=1\),即\(a,p\)互质。
而如果存在\(j\neq i\)使得\(\gcd(b_i,b_j)\neq 1\),则容易知道\(\gcd(b_i,c_i)\neq 1\),自然就不存在\(k_i\)。
这一缺陷来源于CRT的核心思想:构造\(x_i\)使得\(\begin{aligned}\begin{cases} x_i\equiv 1\pmod{b_k} & k=i\\ x_i\equiv 0\pmod{b_k} & k\neq i \end{cases}\end{aligned}\),这就注定了我们要使用乘法逆元。
exCRT则可以处理\(b_i\)不两两互质下的问题。
所谓exCRT,其实与CRT本身没有什么关系,它采用的是另一种思想——将同余式进行合并。
我们考虑\(2\)个同余式构成的方程组:
要将其合并为\(x\equiv a\pmod b\)的形式。
我们将同余式重新表示为:\(x=a_1+k_1\times b_1=a_2+k_2\times b_2(k_1,k_2\in \mathbb{Z})\)。
变形可得:\(k_1\times b_1-k_2\times b_2=a_2-a_1\)。
这是一个典型的不定方程\(ax+by=c\)的形式,未知量是\((k_1,k_2)\),可以用exGCD求解。
如何解$ax+by=c$
具体来说,该方程有解\(\iff \gcd(a,b)\mid c\),必要性显然,充分性由裴蜀定理可得。
然后我们用exGCD求出\(ax+by=\gcd(a,b)\)的解,再将\(x,y\)同时乘\(\frac{c}{\gcd(a,b)}\)即可。
若存在\(ax+by=c\),则可以根据特解\(x,y\)求出任意通解\(x',y'\):
\(\begin{cases}
x'=x+k\times \frac{b}{\gcd(a,b)}\\
y'=y-k\times \frac{a}{\gcd(a,b)}
\end{cases}(k\in \mathbb{Z})\)
若解得特解为\((k_1',k_2')\),通解\((k_1,k_2)\)即为\((k_1'+k\times\frac{b_2}{\gcd(b_1,b_2)},k_2'+k\times\frac{b_1}{\gcd(b_1,b_2)})\),其中\(k\in \mathbb{Z}\)。
选择其一带入:
即:
然后重复上述步骤,直到所有式子合并为\(1\)个为止。
点击查看代码
#include<bits/stdc++.h>
#define int __int128
using namespace std;
const int N=1e5+10;
long long tmp;
int n,a[N],b[N];
int exgcd(int a,int b,int& x,int& y){
int d;
if(b==0) x=1,y=0,d=a;
else d=exgcd(b,a%b,y,x),y-=a/b*x;
return d;
}
signed main(){
ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
cin>>tmp,n=tmp;
for(int i=1;i<=n;i++) cin>>tmp,b[i]=tmp,cin>>tmp,a[i]=tmp;
int pb=b[1],pa=a[1],k1,k2;
for(int i=2;i<=n;i++){
int g=exgcd(pb,b[i],k1,k2);//实际上应该是-b[i],为了避免负数所以这样写,注意此时k2要取相反数
//if((a[i]-pa)%g) cout<<"no solution\n";//保证有解
k1=(a[i]-pa)/g*k1;
pa=pa+pb*k1;
pb=pb/g*b[i];
pa%=pb;
}
tmp=(pa+pb)%pb;
cout<<tmp<<"\n";
return 0;
}
(实际上,也有在CRT的基础上动“小手术”来解决exCRT的方法,见此文,我没有学就不放了)
例题
会不断添加。
P3868 [TJOI2009] 猜数字
变形一下:
然后是CRT的板子了。
点击查看代码
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=15;
int n,fac=1,a[N],b[N];
__int128 ans;
int exgcd(int a,int b,int& x,int& y){
int d;
if(b==0) x=1,y=0,d=a;
else d=exgcd(b,a%b,y,x),y-=a/b*x;
return d;
}
signed main(){
ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
cin>>n;
for(int i=1;i<=n;i++) cin>>a[i];
for(int i=1;i<=n;i++) cin>>b[i],fac*=b[i];
for(int i=1;i<=n;i++){
int c=fac/b[i],k,t;
exgcd(c,b[i],k,t);
(ans+=(__int128)c*k*a[i])%=fac;
}
cout<<(long long)((ans+fac)%fac)<<"\n";
return 0;
}
P4774 [NOI2018] 屠龙勇士
我们发现,用于对付每条龙的剑是固定的,因此我们预处理出对付第\(i\)条龙的剑的攻击力,记为\(b_i\)。
第\(i\)条龙能被击杀,当且仅当\(p_i\mid a_i-b_ix\),也即\(b_ix\equiv a_i\pmod{p_i}\)。
所以我们实际上是要求解下面的线性同余方程组:
我们发现,相比之前的方程,此题在\(x\)前面加了一个系数\(b_i\)。
由于\(p\)不互质,所以采用exCRT。
首先求\(b_i\)在模\(p_i\)意义下的逆元是不可行的,因为\(b_i,p_i\)未必互质。
让我们看看能否沿用之前的思路来分析。
合并成带参的形式
仍然取两个方程,试图将它们合并。
发现\(x\)前系数不同,无法直接合并,因此不妨让系数变成\(\text{lcm}(b_1,b_2)=L\)。
统一系数后,就和exCRT基本相同了,可以类比上面的分析。
将两式联立并整理,得到:
解线性方程组得到特解\((k_1',k_2')\),根据上面的结论,通解\((k_1,k_2)\)可以表示为\((k\in \mathbb{Z})\):
任选其一带入:
变成同余式的形式:
至此完成合并。
最终得到的是一个线性同余方程\(bx\equiv a\pmod p\),理论上应该可以化成不定方程用exGCD来解。不过代码到现在还没写出来,有一些奇怪的问题(^^;
合并成不带参的形式
我们当然更希望最终合并成的同余方程式中,\(x\)系数是\(1\)。
因此假设之前的方程式已经合并成了该形式,现在向其中添加一个带系数的同余方程式:
将\((1)\)式乘\(b_i\),两式联立并整理,得到:
exGCD解不定方程得到\((k_1',k_2')\),通解\((k_1,k_2)\)表示为\((k\in\mathbb{Z})\):
由于我们想要\(x\)的系数为\(1\),所以选择\(k_1\)带入\((1)\)式:
变成同余式的形式:
至此完成合并。
点击查看代码
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=1e5+10;
int t,n,m,a[N],p[N],z[N],b[N];
multiset<int> se;
int exgcd(int a,int b,int& x,int& y){
int d;
if(b==0) x=1,y=0,d=a;
else d=exgcd(b,a%b,y,x),y-=a/b*x;
return d;
}
int solve(){
se.clear();
cin>>n>>m;
for(int i=1;i<=n;i++) cin>>a[i];
for(int i=1;i<=n;i++) cin>>p[i];
for(int i=1;i<=n;i++) cin>>z[i];
for(int i=1,x;i<=m;i++) cin>>x,se.insert(x);
int pp=1,pa=0,C,g,k1,k2,mx=0,tmp;
for(int i=1;i<=n;i++){
auto it=se.upper_bound(a[i]);
if(it!=se.begin()) --it;
b[i]=*it,se.erase(it),se.insert(z[i]);
mx=max(mx,(a[i]-1)/b[i]+1);
}
for(int i=1;i<=n;i++){
g=exgcd((__int128)b[i]*pp%p[i],p[i],k1,k2);
C=(a[i]-b[i]*pa)%p[i];
if(C%g) return -1;
pa+=(__int128)k1*(C/g)%(p[i]/g)*pp%(tmp=p[i]/g*pp);
pa%=(pp=tmp);
}
if(pa<mx) pa+=((mx-pa-1)/pp+1)*pp;
return pa;
}
signed main(){
ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
cin>>t;
while(t--) cout<<solve()<<"\n";
return 0;
}
浙公网安备 33010602011771号