【学习笔记】中国剩余定理

感觉之前写的太垃圾了，重新写一下。

没有那么多奇怪的介绍，直接开始吧。

设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
",CRT(m,r,n));
    }
    return 0;
}