CRT&EXCRT学习笔记

学(EXCRT)的时间挺久了,有点忘了.
写一篇博客记录一下.

---

CRT

首先,我们要知道中国剩余定理是什么
它是用来求解这样一个同余方程的

[xequiv a_1(mod p_1)\ xequiv a_2(mod p_2)\ ...\ xequiv a_n(mod p_n) ]

其中,(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);
}