学(EXCRT)的时间挺久了,有点忘了.
写一篇博客记录一下.
CRT
首先,我们要知道中国剩余定理是什么
它是用来求解这样一个同余方程的
其中,(p_1,p_2...p_n)两两互质.
我们先规定一些数.(M=prod_{i=1}^np_i,L=lcm(p_1,p_2...p_n))
(CRT)的流程是这样的.
假设现在执行到(i)
令(m=frac{M}{p_i})
找到一个(x),使得(mxequiv 1(mod p_i))
找(x)可以用逆元的方法实现.
那么,显然地,(a_imxequiv a_i(mod p_i))
我们令(t_i=a_imx)
那么,答案就是((sum_{i=1}^nt_i)\% L)
这个证明显然吧
EXCRT
有的时候,出题人并不保证所有的(p_i)互质.
那么我们该怎么办呢?
其实也不难.
假设有(2)个方程(xequiv c_1(mod p_1),xequiv c_2(mod p_2))
我们可以列出相应的方程.
(x=c_1+k_1p_1,x=c_2+k_2p_2)
那么(c_1+k_1p_1=c_2+k_2p_2)
(k_1p_1=c_2-c_1+k_2p_2)
令(g=gcd(p_1,p_2)),那么(frac{p_1}{g}k_1=frac{c_2-c_1}{g}+frac{p_2}{g}k_2).
因此只有(g)能整除(c_2-c_1)时才可以合并,否则无解.
此时,可以化简方程为(frac{p_1}{g}k_1equivfrac{c_2-c_1}{g}(mod frac{p_2}{g})).
由于((frac{p_1}{g},frac{p_2}{g})=1),我们可以用(exgcd)求出(frac{p_1}{g})在模(frac{p_2}{g})意义下的逆元,假设为(inv).
那么(k_1equivfrac{c_2-c_1}{g}*inv(mod frac{p_2}{g})).
我们把这个方程回代到第一个方程中,得到
(xequiv c_1+(frac{c_2-c_1}{g}*inv\% frac{p_2}{g})(mod frac{p_1*p_2}{g}))
但是,这个方程中所有的量都是已知的,因此我们就一直合并下去即可.
代码如下
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<vector>
#define N (100010)
#define inf (0x7f7f7f7f)
#define rg register int
#define Label puts("NAIVE")
#define spa print(' ')
#define ent print('
')
#define rand() (((rand())<<(15))^(rand()))
typedef long double ld;
typedef long long LL;
typedef unsigned long long ull;
using namespace std;
inline char read(){
static const int IN_LEN=1000000;
static char buf[IN_LEN],*s,*t;
return (s==t?t=(s=buf)+fread(buf,1,IN_LEN,stdin),(s==t?-1:*s++):*s++);
}
template<class T>
inline void read(T &x){
static bool iosig;
static char c;
for(iosig=false,c=read();!isdigit(c);c=read()){
if(c=='-')iosig=true;
if(c==-1)return;
}
for(x=0;isdigit(c);c=read())x=((x+(x<<2))<<1)+(c^'0');
if(iosig)x=-x;
}
inline char readchar(){
static char c;
for(c=read();!isalpha(c);c=read())
if(c==-1)return 0;
return c;
}
const int OUT_LEN = 10000000;
char obuf[OUT_LEN],*ooh=obuf;
inline void print(char c) {
if(ooh==obuf+OUT_LEN)fwrite(obuf,1,OUT_LEN,stdout),ooh=obuf;
*ooh++=c;
}
template<class T>
inline void print(T x){
static int buf[30],cnt;
if(x==0)print('0');
else{
if(x<0)print('-'),x=-x;
for(cnt=0;x;x/=10)buf[++cnt]=x%10+48;
while(cnt)print((char)buf[cnt--]);
}
}
inline void flush(){fwrite(obuf,1,ooh-obuf,stdout);}
int n; LL c,m,c2,mo;
LL mult(LL a,LL b,LL mod){//a*b
LL res=0,fu=1;
if(a<0)fu=-fu,a=-a;
if(b<0)fu=-fu,b=-b;
while(b){
if(b&1)res=(res+a)%mod;
a=(a+a)%mod,b>>=1;
}
res*=fu;
if(res<0)(res+=((-res-1)/mod+1)*mod);
return res;
}
void exgcd(LL a,LL b,LL &x,LL &y){
if(!b){x=1,y=0;return;}
exgcd(b,a%b,y,x),y-=(a/b)*x;
}
LL gcd(LL a,LL b){return b?gcd(b,a%b):a;}
LL inv(LL a,LL b){
LL x,y;
exgcd(a,b,x,y);
return (x<=0)?(x+b):x;
}
int main(){
read(n),read(m),read(c);
for(int i=1;i<n;i++){
read(mo),read(c2);
LL t=gcd(m,mo),s=inv(m/t,mo/t),tc=c,tm=m;
m=(mo/t*tm),c=(tc+mult(tm,mult(s,(c2-c)/t,(mo/t)),m))%m;
}
if(c<0)c+=((-c-1)/m+1)*m;
printf("%lld
",c%m);
}