时隔多年,再次复习 CRT(数论全家桶)
作者:互联网
1.CRT
中国剩余定理,用来求解同余方程组
\begin{cases} x\equiv a_1\pmod {m_1} \\ x\equiv a_2\pmod {m_2} \\ x\equiv a_3\pmod {m_3} \\ …………\\ x\equiv a_n\pmod {m_n}\\ \end{cases}⎩⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎧x≡a1(modm1)x≡a2(modm2)x≡a3(modm3)…………x≡an(modmn)
其中,mm 两两互质。
显然,这个柿子一定有无数个解,我们只需要找到最小解即可。
中国剩余定理告诉我们,这个柿子存在一个通解可以计算。
x=\sum_{i=1}^n a_iM_it_i,(M_i=(\prod_{j=1}^nm_j)/mi,t_i=inv(M_i,m_i))x=i=1∑naiMiti,(Mi=(j=1∏nmj)/mi,ti=inv(Mi,mi))
对于逆元的计算,我们可以使用扩展欧几里德算法来计算。
int exgcd(int a,int b,int &x,int &y){
if(b==0){x=1,y=0;return a;}
int ans=exgcd(b,a%b,x,y);
int t=x;x=y,y=t-a/b*y;
return ans;
}
为什么,因为我们构造一个特殊的方程:
M_it_i+m_iy=\gcd(m_i,M_i)Miti+miy=gcd(mi,Mi)
得到这个方程的解,因为两两互质,显然 \gcd(m_i,M_i)=1gcd(mi,Mi)=1
于是方程化为了:
M_it_i+m_iy=1Miti+miy=1
我们要求的逆元是 inv(M_i,m_i)inv(Mi,mi)
于是我们化简简单移项原方程:
M_it_i=1-m_iyMiti=1−miy
再次化简得:
M_it_i\equiv 1 \pmod {m_i}Miti≡1(modmi)
于是我们就可以得知:
M_it_i+m_iy=\gcd(m_i,M_i)Miti+miy=gcd(mi,Mi)
的特解中的t_iti,就是我们要求的逆元。
那么我们反观我们的定理:
\begin{cases} x\equiv a_1\pmod {m_1} \\ x\equiv a_2\pmod {m_2} \\ x\equiv a_3\pmod {m_3} \\ …………\\ x\equiv a_n\pmod {m_n}\\ \end{cases}⎩⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎧x≡a1(modm1)x≡a2(modm2)x≡a3(modm3)…………x≡an(modmn)
x=\sum_{i=1}^n a_iM_it_i,(M_i=(\prod_{j=1}^nm_j)/mi,t_i=inv(M_i,m_i))x=i=1∑naiMiti,(Mi=(j=1∏nmj)/mi,ti=inv(Mi,mi))
是不是感觉清爽很多了?
那么我们尝试证明一下。
证明:
因为 M_i=(\prod_{j=1}^nm_j)/miMi=(∏j=1nmj)/mi ,即除了 m_imi 以外所有模数的倍数。
所以 \forall j,a_iM_it_i\equiv0\pmod {m_k}∀j,aiMiti≡0(modmk) ,因此不会对其他的方程造成影响,这很好,给我们的定理证明创造了机会。
又因为 t_i=inv(M_i,m_i)ti=inv(Mi,mi),即 M_it_i\equiv1\pmod {m_i}Miti≡1(modmi)
所以 a_iM_it_i\equiv a_i\pmod {m_i}aiMiti≡ai(modmi)
因此累加每一个方程组的答案即可得到方程组的通解。
若答案过大,对 xx 取模 \prod_{j=1}^nm_j∏j=1nmj 即可。
时间复杂度 O(n\log n)O(nlogn) 主要是扩展欧几里得算法和循环带来的复杂度。
附上模板代码:
#include<cstdio>
#include<algorithm>
#define int long long
#define N 1919
using namespace std;
int n,a[N],b[N];
int exgcd(int a,int b,int &x,int &y){
if(b==0){x=1,y=0;return a;}
int ans=exgcd(b,a%b,x,y);
int t=x;x=y,y=t-a/b*y;
return ans;
}
int crt(int n,int *a,int *m){
int M=1,ans=0;
for(int i=1;i<=n;i++)M*=m[i];
for(int i=1;i<=n;i++){
int mi=M/m[i],ti=0,y=0;
exgcd(mi,m[i],ti,y);
ans=(ans+(a[i]*mi*ti)%M)%M;
}
return (ans+M)%M;
}
signed main(){
scanf("%lld",&n);
for(int i=1;i<=n;i++)scanf("%lld%lld",&a[i],&b[i]);
printf("%lld",crt(n,b,a));
return 0;
}
2.扩展CRT
扩展 CRT 和 CRT 的思想完全不一样,扩展 CRT 更像是在合并同余方程打到化简的效果。
\begin{cases} x\equiv a_1\pmod {m_1} \\ x\equiv a_2\pmod {m_2} \\ x\equiv a_3\pmod {m_3} \\ …………\\ x\equiv a_n\pmod {m_n}\\ \end{cases}⎩⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎧x≡a1(modm1)x≡a2(modm2)x≡a3(modm3)…………x≡an(modmn)
当 mm 不再满足互质的关系时,我们应该怎么做?
观察一个二元一次同余方程:
\begin{cases} x\equiv a_1\pmod {m_1} \\ x\equiv a_2\pmod {m_2} \\ \end{cases}{x≡a1(modm1)x≡a2(modm2)
可以写作:
\begin{cases} x=a_1+y_1m_1 \\ x=a_2+y_2m_2 \\ \end{cases}{x=a1+y1m1x=a2+y2m2
整理得:
x=a_1+y_1m_1=a_2+y_2m_2x=a1+y1m1=a2+y2m2
继续化简:
a_1+y_1m_1=a_2+y_2m_2a1+y1m1=a2+y2m2
y_1m_1+zm_2=a_2-a_1 ( z=-y_2 )y1m1+zm2=a2−a1(z=−y2)
这个柿子在 a_2-a_1a2−a1 可以通过扩展欧几里得算法得到特解 (y,z)(y,z)。
很显然,对于原方程:
m_1y \frac{a_2-a_1}{\gcd(m_1,m_2)}+m_2z \frac{a_2-a_1}{\gcd(m_1,m_2)} = \frac{a_2-a_1}{\gcd(m_1,m_2)}\gcd(m_1,m_2)=a_2-a_1m1ygcd(m1,m2)a2−a1+m2zgcd(m1,m2)a2−a1=gcd(m1,m2)a2−a1gcd(m1,m2)=a2−a1
因为 x=a_1+y_1m_1x=a1+y1m1,现在让 xx 尽可能小,其实是让 y\frac{a_2-a_1}{\gcd(m_1,m_2)}ygcd(m1,m2)a2−a1 尽可能小
我们现在要最小化 y\frac{a_2-a_1}{\gcd(m_1,m_2)}ygcd(m1,m2)a2−a1 不妨设为 tt。
t=y_1+k\frac{m_2}{\gcd(m_1,m_2)}t=y1+kgcd(m1,m2)m2
带回发现:
x=a_1+y_1m_1+k\times lcm(m_1,m_2)x=a1+y1m1+k×lcm(m1,m2)
要删掉后面那堆,显然只用对 tt 进行取模即可:
t=y_1+k\frac{m_2}{\gcd(m_1,m_2)}\pmod {\frac{m_2}{\gcd(m_1,m_2)}}t=y1+kgcd(m1,m2)m2(modgcd(m1,m2)m2)
于是我们得到和合并后的方程组:
x\equiv m_1y_1+a_1\pmod {lcm(m_1,m_2)}x≡m1y1+a1(modlcm(m1,m2))
#include<cstdio>
#include<algorithm>
#include<cstdlib>
#include<cmath>
#define int __int128
using namespace std;
long long n;
int exgcd(int a,int b,int &x,int &y){
if(b==0){x=1,y=0;return a;}
int ans=exgcd(b,a%b,x,y);
int t=x;x=y,y=t-a/b*y;
return ans;
}
int lcm(int a,int b){return a*b/__gcd(a,b);}
pair<int,int> merge(int a2,int m2,int a1,int m1){
int x,y,c;
int d=exgcd(m1,m2,x,y);
c=a2-a1;
if(c%d){puts("-1");exit(0);}
x=x*c/d%(m2/d);
if(x<0)x+=(m2/d);
int md=lcm(m1,m2);
a1=(m1*x+a1)%md;
if(a1<0)a1+=md;
m1=md;
return {a1,m1};
}
signed main(){
scanf("%lld",&n);
int a1=0,m1=0;
for(int i=1;i<=n;i++){
long long a,m;
scanf("%lld%lld",&m,&a);
if(i==1)a1=a,m1=m;
else{
pair<int,int> fi=merge(a1,m1,a,m);
a1=fi.first,m1=fi.second;
}
}
printf("%lld",(long long)(a1%m1));
return 0;
}
标签:复习,CRT,int,a1,pmod,m1,m2,时隔多年,equiv 来源: https://www.cnblogs.com/masida-hall-LZ/p/16629495.html