扩展中国剩余定理(EXCRT)
数论基础知识吧
主要还是留个板子,之后就不用动脑子了(
模板:
//在主函数依次调用 //fold(a,b,m,n)和excrt(b,m,n) //求解ax+by=gcd(a,b),返回gcd(a,b) ll exgcd(ll a,ll b,ll &x,ll &y) { if(b==0LL) { x=1,y=0; return a; } ll tmp=exgcd(b,a%b,y,x); y=y-a/b*x; return tmp; } inline ll mod(ll x,ll m) { if(x>=m) x-=m; if(x<0) x+=m; return x; } //防止中间结果爆longlong的x*y%m ll mul(ll x,ll y,ll m) { ll res=0,tmp=x; x=mod(x%m+m,m),y=mod(y%m+m,m); while(y) { if(y&1LL) res=mod(res+tmp,m); tmp=mod(tmp*2,m); y/=2; } return res; } //求解方程组x_i=a_i mod m_i ll excrt(ll *a,ll *m,int n) { ll res=0,lcm=1; for(int i=1;i<=n;i++) { ll A,B,C,x,y,gcd; A=lcm,B=m[i],C=a[i]-res; gcd=exgcd(A,B,x,y); //无解 if(C%gcd!=0) return -1; ll tmp=lcm; lcm=lcm/gcd*m[i]; res=res+mul(C/gcd,mul(tmp,x,lcm),lcm); res=mod(res%lcm+lcm,lcm); } return res; } //将a_i的逆元乘到b_i上 //之后要excrt的是b_i与m_i ll fold(ll *a,ll *b,ll *m,int n) { for(int i=1;i<=n;i++) { ll x,y,gcd; gcd=exgcd(a[i],m[i],x,y); //无逆元 if(b[i]%gcd!=0) return -1; a[i]/=gcd,b[i]/=gcd,m[i]/=gcd; exgcd(a[i],m[i],x,y); b[i]=mul(b[i],x,m[i]); } return 1; }
主要用来求解线性同余方程组
即我们希望对于以下方程组求最小的$x$,使得
\[\begin{equation} \begin{cases} x\equiv a_1 (mod\ \ m_1)\\ x\equiv a_2(mod\ \ m_2)\\ \vdots \\ x\equiv a_n(mod\ \ m_n)\end{cases} \end{equation}\nonumber \]
其中$m_1,m_2,...,m_n$不保证互质
在此之前需要先掌握扩展欧几里德(exgcd)
参考这篇就足够了,很详细(最后一步式子错了,不过结论对的):RiverHamster - exgcd
exgcd解决的问题是:$ax+by=gcd(a,b)$
//求解ax+by=gcd(a,b),返回gcd(a,b) ll exgcd(ll a,ll b,ll &x,ll &y) { if(b==0LL) { x=1,y=0; return a; } ll tmp=exgcd(b,a%b,y,x); y=y-a/b*x; return tmp; }
一些除了本身之外的用法:
1. 解$ax+by=c$的不定方程:将解出的$x,y$同乘$\frac{c}{gcd(a,b)}$即可得到特解;若$c\ mod\ gcd(a,b)\neq 0$,则方程无解
2. 求逆元:若$a,b$互质,那么求出的$x$就是$a$在模$b$意义下的逆元了
3. 解$Ax\equiv B(mod\ \ M)$的同余方程:令$a=A,b=M,c=B$,则$X=x\cdot \frac{c}{gcd(a,b)}$为该方程的解;若$c\ mod\ gcd(a,b)\neq 0$,则该方程无解
有了上面的铺垫,现在可以愉快地求同余方程组了
参考这篇:niiick - 题解 P4777 【【模板】扩展中国剩余定理(EXCRT)】
假设现在已经算完了前$k-1$个方程组的解,设$M=lcm(m_1,m_2,...,m_{k-1})$
则通解为$x+t\cdot M$($t$为任意整数)
对于第$k$个方程,我们希望$x+t\cdot M\equiv a_k(mod\ \ m_k)$
即,希望找到合适的$t,s$使得$x+t\cdot M=a_k+s\cdot m_k$
移项后就是$t\cdot M+s\cdot m_k=a_k-x$,已经是exgcd可以处理的$ax+by=c$的形式了;用$exgcd(M,m_k,t,s)$,再乘以$\frac{a_k-x}{gcd(M,m_k)}$,即可得到$t,s$
于是前$k$个方程的特解为$x'=x+t\cdot M$,如令$M'=lcm(M,m_k)$则通解为$x'+t'M'$($t'$为任意整数)
这样一步一步做下去,直到算出前$n$个方程的特解
于是可以写出下面的代码(题目为:Luogu P4777 (【模板】扩展中国剩余定理(EXCRT)))
#include <cstdio> #include <cstring> #include <algorithm> using namespace std; typedef long long ll; const int N=100005; //求解ax+by=gcd(a,b),返回gcd(a,b) ll exgcd(ll a,ll b,ll &x,ll &y) { if(b==0LL) { x=1,y=0; return a; } ll tmp=exgcd(b,a%b,y,x); y=y-a/b*x; return tmp; } inline ll mod(ll x,ll m) { if(x>=m) x-=m; if(x<0) x+=m; return x; } //防止中间结果爆longlong的x*y%m ll mul(ll x,ll y,ll m) { ll res=0,tmp=x; x=mod(x%m+m,m),y=mod(y%m+m,m); while(y) { if(y&1LL) res=mod(res+tmp,m); tmp=mod(tmp*2,m); y/=2; } return res; } //求解方程组x_i=a_i mod m_i ll excrt(ll *a,ll *m,int n) { ll res=0,lcm=1; for(int i=1;i<=n;i++) { ll A,B,C,x,y,gcd; A=lcm,B=m[i],C=a[i]-res; gcd=exgcd(A,B,x,y); //无解 if(C%gcd!=0) return -1; ll tmp=lcm; lcm=lcm/gcd*m[i]; res=res+mul(C/gcd,mul(tmp,x,lcm),lcm); res=mod(res%lcm+lcm,lcm); } return res; } int n; ll a[N],m[N]; int main() { scanf("%d",&n); for(int i=1;i<=n;i++) scanf("%lld%lld",&m[i],&a[i]); printf("%lld\n",excrt(a,m,n)); return 0; }
也可以有再稍微扩展一点的做法,即解决
\[\begin{equation} \begin{cases} a_1x\equiv b_1 (mod\ \ m_1)\\ a_2x\equiv b_2(mod\ \ m_2)\\ \vdots \\ a_nx\equiv b_n(mod\ \ m_n)\end{cases} \end{equation}\nonumber \]
当$gcd(a_i,m_i)=1$时,可以通过exgcd求逆元
当$gcd(a_i,m_i)\neq 1$,若$b_i\ mod\ gcd(a_i,m_i)\neq 0$,则该方程必然无解;否则将$a_i,b_i,m_i$同除$gcd(a_i,m_i)$,则$\frac{a_i}{gcd(a_i,m_i)},\frac{m_i}{gcd(a_i,m_i)}$互质,那么也可以通过exgcd求逆元
之前看错题意时想的,没找到题目
#include <cstdio> #include <cstring> #include <algorithm> using namespace std; typedef long long ll; const int N=100005; //求解ax+by=gcd(a,b),返回gcd(a,b) ll exgcd(ll a,ll b,ll &x,ll &y) { if(b==0LL) { x=1,y=0; return a; } ll tmp=exgcd(b,a%b,y,x); y=y-a/b*x; return tmp; } inline ll mod(ll x,ll m) { if(x>=m) x-=m; if(x<0) x+=m; return x; } //防止中间结果爆longlong的x*y%m ll mul(ll x,ll y,ll m) { ll res=0,tmp=x; x=mod(x%m+m,m),y=mod(y%m+m,m); while(y) { if(y&1LL) res=mod(res+tmp,m); tmp=mod(tmp*2,m); y/=2; } return res; } //求解方程组x_i=a_i mod m_i ll excrt(ll *a,ll *m,int n) { ll res=0,lcm=1; for(int i=1;i<=n;i++) { ll A,B,C,x,y,gcd; A=lcm,B=m[i],C=a[i]-res; gcd=exgcd(A,B,x,y); //无解 if(C%gcd!=0) return -1; ll tmp=lcm; lcm=lcm/gcd*m[i]; res=res+mul(C/gcd,mul(tmp,x,lcm),lcm); res=mod(res%lcm+lcm,lcm); } return res; } //将a_i的逆元乘到b_i上 //之后要excrt的是b_i与m_i ll fold(ll *a,ll *b,ll *m,int n) { for(int i=1;i<=n;i++) { ll x,y,gcd; gcd=exgcd(a[i],m[i],x,y); //无逆元 if(b[i]%gcd!=0) return -1; a[i]/=gcd,b[i]/=gcd,m[i]/=gcd; exgcd(a[i],m[i],x,y); b[i]=mul(b[i],x,m[i]); } return 1; } int n; ll a[N],b[N],m[N]; int main() { scanf("%d",&n); for(int i=1;i<=n;i++) scanf("%lld%lld",&m[i],&b[i]),a[i]=1; fold(a,b,m,n); printf("%lld\n",excrt(b,m,n)); return 0; }
只是留个模板,毕竟玩不出什么花里胡哨的
(完)