中国剩余定理
看介绍和解方程:https://baike.baidu.com/item/%E5%AD%99%E5%AD%90%E5%AE%9A%E7%90%86
前提:
Mi互质
例题:POJ1006
代码:
该代码只能处理模数互质的情况,对于模数不互质的情况,可以依此合并一次同余同余方程来求解。详情看最后kuangbin的模板
过了一年后看自己写的代码感觉好丑而且不使用
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=1e4+10;
ll ExGcd(ll a,ll b,ll &x,ll &y)//求 a b 最大公约数,且得到gcd(a,b)=x*a+y*b;
{
if(!b)
{
x=1;
y=0;
return a;
}
ll gcd=ExGcd(b,a%b,x,y);
ll temp,k=a/b;
temp=x;
x=y;
y=temp-k*y;
return gcd;
}
bool IsOk;
ll calc(ll a,ll c,ll m)
{
ll x,y;
ll gcd=ExGcd(a,m,x,y);
if(c%gcd!=0)
{
IsOk=false;
return 0ll;
}
return c/gcd*x%m;
}
/* 前提 mi互质 中国剩余定理 x=a1(mod m1) x=a2(mod m2) x=an(mod mn) 设M为m1...mn大的乘积 Mi=M/mi 设ti为Mi关于mi的逆元 则 ans=ai*Mi*ti 时间复杂度为O(n*gcd) */
ll aa[maxn],mm[maxn];
ll CRT(int n)//n个方程 返回x的的值
{
ll Mpro=1;
for(int i=0; i<n; ++i)
Mpro*=mm[i];
ll ans=0;
for(int i=0; i<n; ++i)
{
aa[i]%=mm[i];
ans+=(aa[i]*calc(Mpro/mm[i],1,mm[i])%Mpro*(Mpro/mm[i]))%Mpro;
}
return ans;
}
int main()
{
ll a,b,c,d;
mm[0]=23,mm[1]=28,mm[2]=33;
int cas=0;
ll ans=0;
while(cin>>a>>b>>c>>d)
{
if(a==-1)
break;
aa[0]=a;
aa[1]=b;
aa[2]=c;
ans=CRT(3);
ll Mpro=23*28*33;
ans=((ans-d)%Mpro+Mpro)%Mpro;
cout<<"Case "<<++cas<<": the next triple peak occurs in ";
cout<<(ans==0?Mpro:ans)<<" days."<<endl;
}
return 0;
}
kuangbin的代码,处理模数不互质的情况。
long long extend_gcd(long long a,long long b,long long &x,long long &y)
{
if(a == 0 && b == 0)return -1;
if(b ==0 )
{
x = 1;
y = 0;
return a;
}
long long d = extend_gcd(b,a%b,y,x);
y -= a/b*x;
return d;
}
int m[10],a[10];//模数为m,余数为a, X % m = a
bool solve(int &m0,int &a0,int m,int a)
{
long long y,x;
int g = extend_gcd(m0,m,x,y);//得到的x是m0/g关于m1/g的逆元
if( abs(a - a0)%g )return false;
x *= (a - a0)/g;
x %= m/g;
a0 = (x*m0 + a0);
m0 *= m/g; //合并方程后的模数是增大的
a0 %= m0;
if( a0 < 0 )a0 += m0;
return true;
}
/* * 无解返回false,有解返回true; * 解的形式最后为 a0 + m0 * t (0<=a0<m0) */
bool MLES(int &m0,int &a0,int n) //解为 X = a0 + m0 * k
{
bool flag = true;
m0 = 1;
a0 = 0;
for(int i = 0; i < n; i++)
if( !solve(m0,a0,m[i],a[i]) )
{
flag = false;
break;
}
return flag;
}