【学习笔记】中国剩余定理
感觉之前写的太垃圾了,重新写一下。
没有那么多奇怪的介绍,直接开始吧。
设m1,m2,m3,m4....mk两两互素,则同余方程组
x≡a1(mod m1)
x≡a2(mod m2)
x≡a3(mod m3)
x≡a4(mod m4)
x≡ak(mod mk)
一定有解,x≡(a1*M1*M1^(-1)+a2*M2*M2^(-1)+....)(mod M)
其中M=m1*m2*m3*....mk,Mi=M/mi,Mi^(-1)是Mi在模mi意义下的逆元。
互素情况下的模板,是POJ1006代码
代码:
#include<iostream> #include<cstring> #include<cstdio> using namespace std; int a[4],m[4]; int p,e,i,d,t=1; void exgcd(int a,int b,int &x,int &y){ if(b==0){ x=1;y=0; return; } exgcd(b,a%b,x,y); int t; t=x;x=y;y=t-a/b*y; } int CRT(int a[],int m[],int n){ int M=1,ans=0; for(int i=1;i<=n;i++)M*=m[i];//求出M for(int i=1;i<=n;i++){ int x,y; int Mi=M/m[i]; exgcd(Mi,m[i],x,y);//扩展欧几里得求逆元 ans=(ans+Mi*x*a[i])%M;//统计答案 } if(ans<0)ans+=M; return ans; } int main(){ while(cin>>p>>e>>i>>d){ if(p==-1&&e==-1&&i==-1&&d==-1)break; a[1]=p;a[2]=e;a[3]=i; m[1]=23;m[2]=28;m[3]=33; int ans=CRT(a,m,3); if(ans<=d)ans+=21253; } return 0; }
接下来是m1,m2,m3,m4...mk不互素的情况。
主要的做法是将式子合并。
对于
x mod m1=a1
x mod m2=a2
x mod m3=a3
可以将前两个式子写成
x=k1*m1+a1--------- (1)
x=k2*m2+a2----------(2)
那么k1*m1+a1=k2*m2+a2
整理得:
k1*m1-k2*m2=a2-a1,
对于上面这个式子,其实就是ax+by=k。
用扩展欧几里得求出k1之后,代入(1),就可以求出x。
这里求的x是(1)、(2)两式的解,通解是x'=x+k*lcm(m1,m2)。
将通解的式子变形就是x' mod lcm(m1,m2)=x。这个式子就是我们将前两个方程
合并后的式子。
令M=lcm(m1,m2),A=x*m1+a1.
那么新的方程就是 x' mod M =A
n个方程一直合并就可以了,最后剩下一个方程
x mod M=A,其中M=lcm(m1,m2,m3....mk)。
所求x是方程组的一个特解,通解为x'=x+k*M。
注意:没有解的条件是gcd(m1,m2)|a2-a1是不成立的。
如果要求最小正整数解的话:
x%=M;if(x<0)x+=M;
代码:
POJ2891
#include<iostream> #include<cstdio> #include<cstring> #define LL long long #define maxn 100009 using namespace std; int n; LL m[maxn],r[maxn]; LL exgcd(LL a,LL b,LL &x,LL &y){ if(b==0){ x=1;y=0; return a; } LL t,r=exgcd(b,a%b,x,y); t=x;x=y;y=t-a/b*y; return r; } LL CRT(LL m[],LL r[],LL n){ LL x,y,gcd,M=m[1],R=r[1]; for(int i=2;i<=n;i++){ LL gcd=exgcd(M,m[i],x,y); if((r[i]-R)%gcd)return -1; x=(r[i]-R)/gcd*x%(m[i]/gcd); R+=x*M;//新的余数R=r1+x*m[1] M=M/gcd*m[i];//新的M,lcm(m1,m2) R%=M; } return R>0?R:R+M; } int main(){ while(~scanf("%d",&n)){ for(int i=1;i<=n;i++) scanf("%lld%lld",&m[i],&r[i]); printf("%lld\n",CRT(m,r,n)); } return 0; }