数论三连3·CRT和EXCRT

蒟蒻的我终于调出了屠龙勇士。。数论三连终于可以集齐了，留念一下。。。

CRT：构造啊构造。。害怕。。

有一堆形如ai*x=ci mod bi的东西。（ai,bi两两互质）

设M=b1*b2*..*bn;Mi=M/bi； 那么其实就是求解ai*(Mi*x')=ci mod bi->x'=ci*inv(ai*Mi)， 

于是通解x=sum(ci*inv(ai*Mi)*Mi)+k*M。

EXCRT：其实我觉得CRT还是EXCRT好写。害怕。。

首先可以将ai*x mod ci=bi转化成x mod ..=..的形式。

怎么转化呢？ ai*x+bi*y=ci; ->d=gcd(ai,bi,ci); 将aibici都除以d，

则ai'*x+bi'*y=ci'；此时ai'和bi'一定互质，否则就无解。

为什么？如果ai'和bi'仍不互质，则gcd(gcd(ai',bi'),ci')=1，而gcd(ai',bi')!=1，那么gcd(ai',bi')必然不能整除ci，方程无解。

化到这里以后，就可以求逆元转化成x =ci*inv(ai') mod bi。

然后使用数学归纳法。前i个方程的通解形式可以表现为x0+k*m，其中x0是一个已求好的特解，m是lcm(bi)。

然后对于下一个方程j，必然有(x0+k*m)=cj mod bj -> mk=cj-x0 mod bj，exgcd求出随便某个k，x0'=x0+k*m，就得到了新的特解。

最后将特解模m就可以了。

为何可以表现成x0+k*m这种形式？首先对于第一个方程当然是成立的。对于后面的方程，假设前i个方程可以这样表示，

mk=cj-x0 mod bj->mk+by=c-x0; 这个方程的通解形式是k=k0+t*(lcm(m,b)/m)；于是对于所有x0+k*m都可以表示成

x0+k0*m+t*lcm(m,b)=x0'+t*lcm(m,b)，证完了。

最后吐槽一下调不出屠龙勇士的原因->在使用multiset极其傻叉地erase了元素值。。。是菜的啊。。。

luogu4777的EXCRT模板：

#include<bits/stdc++.h>
using namespace std;
#define rep(x,y,z) for (int x=y; x<=z; x++)
#define downrep(x,y,z) for (int x=y; x>=z; x--)
#define ms(x,y,z) memset(x,y,sizeof(z))
#define LL long long
#define repedge(x,y) for (int x=hed[y]; ~x; x=edge[x].nex)
inline LL read(){
	LL x=0; int w=0; char ch=0;
	while (ch<'0' || ch>'9') w|=ch=='-',ch=getchar();
	while (ch>='0' && ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
	return w? -x:x;
}
const int N=100005;
int n;
struct node{ LL a,b; }o[N];
LL multi(LL a,LL b,LL c){
	a=(a%c+c)%c; b=(b%c+c)%c; 
	LL res=0; while (b){ if (b&1) res=(res+a)%c; a=(a+a)%c; b>>=1; }
	return (res%c+c)%c;
}
void exgcd(LL a,LL b,LL c,LL &x,LL &y){
	if (!a){ x=0; y=c/b; return; }
	if (!b){ y=0; x=c/a; return; }
	exgcd(b,a%b,c,y,x);
	y-=(a/b)*x;
}
LL gcd(LL a,LL b){ return b? gcd(b,a%b):a; }
LL lcm(LL a,LL b){ return (a/gcd(a,b))*b; }
LL excrt(node *o,int n){
	//x'=x+k*m=a (mod b) mod lcm(m,b); m*k=a-x (mod b);
	LL m=o[1].b; LL x=(o[1].a%m+m)%m;
	rep(i,2,n){
		o[i].a%=o[i].b; LL newm=lcm(m,o[i].b); 
		LL a=m%o[i].b; LL b=o[i].b; LL c=((o[i].a-x)%o[i].b+o[i].b)%o[i].b; 
	    LL d=gcd(a,b); if (c%d) return -1; a/=d; b/=d; c/=d;
		LL X,Y; exgcd(a,b,1,X,Y); X=multi(X,c,newm); 
		x=((x+multi(X,m,newm))%newm+newm)%newm; m=newm;
	}   
	return x;
}
int main(){
	n=read(); rep(i,1,n) o[i].b=read(),o[i].a=read();
	printf("%lld\n",excrt(o,n)); return 0;
}

NOI2018屠龙勇士：（吐槽：网络赛的时候一看此题突然脑子一抽想到了同余方程。。然后当时不知道EXCRT。。手推了

鬼畜的exgcd反复带入解同余方程的奇怪东西。。。怒调1小时过了样例就很开心地交了，连快速乘都没加，然后成功被卡成0分，然而EXCRT是多么地简洁，真是一个悲伤的故事）

#include<bits/stdc++.h>
using namespace std;
#define rep(x,y,z) for (int x=y; x<=z; x++)
#define downrep(x,y,z) for (int x=y; x>=z; x--)
#define ms(x,y,z) memset(x,y,sizeof(z))
#define LL long long
#define repedge(x,y) for (int x=hed[y]; ~x; x=edge[x].nex)
inline LL read(){
	LL x=0; int w=0; char ch=0;
	while (ch<'0' || ch>'9') w|=ch=='-',ch=getchar();
	while (ch>='0' && ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
	return w? -x:x;
}
const int N=200005;
int n,m;
multiset<LL> s;
LL ATK[N],a[N],b[N],c[N];
struct node{ LL a,b; }o[N];
LL multi(LL a,LL b,LL c){
	LL res=0; LL w=(b<0)? (-1):1; if (b<0) b=-b;
	while (b){ if (b&1) res=(res+a)%c; a=(a+a)%c; b>>=1; } 
	return ((res*w)%c+c)%c;
}
LL gcd(LL a,LL b){ return b? gcd(b,a%b):a; }
LL lcm(LL a,LL b){ return (a/gcd(a,b))*b; }
void exgcd(LL a,LL b,LL c,LL &x,LL &y,LL d){
	if (!a){ x=0; y=c/b; return; }
	if (!b){ y=0; x=c/a; return; }
	exgcd(b,a%b,c,y,x,d); y=(y-multi(a/b,x,d))%d;
}
//x=a mod b
LL excrt(node *o,int n){
	LL x=0; LL m=1; 
	rep(i,1,n){ 
		LL a=m%o[i].b; LL b=o[i].b; LL c=(o[i].a-x)%o[i].b;
		if (!a){ m=lcm(m,o[i].b); if (c) return -1; continue; }
		LL d=gcd(a,b); if (c%d) return -1; a/=d; b/=d; c/=d;
		LL newm=lcm(m,o[i].b); LL X,Y; exgcd(a,b,1,X,Y,newm); X=multi(X,c,newm);
		x=((x+multi(X,m,newm))%newm+newm)%newm; m=newm;  
	}
	x=(x%m+m)%m; LL tmp=0;
	double res=0; rep(i,1,n) res=max(res,(double)c[i]/a[i]); 
	rep(i,1,n){
		double t1=c[i]-a[i]*x; double t2=a[i]*m; 
		tmp=max(tmp,(LL)ceil(t1/t2));
	}
	return x+tmp*m;
}
LL inv(LL u,LL c){
	LL x,y; exgcd(u,c,1,x,y,c); return (x%c+c)%c; 
}
LL solve(int t){
	LL x=a[t]; LL y=b[t]; LL z=c[t];
	LL d=gcd(gcd(x,y),z); x/=d; y/=d; z/=d; 
	if (abs(gcd(x,y))!=1) return -1; 
	o[t].a=multi(z,inv(x,y),y); o[t].b=y;
}
int main(){
	int cas=read();
	rep(cl,1,cas){
		n=read(); m=read(); s.clear();
		rep(i,1,n) c[i]=read(); rep(i,1,n) b[i]=read(); 
		rep(i,1,n) ATK[i]=read(); rep(i,1,m) { LL x=read(); s.insert(x); }
		rep(i,1,n){
			multiset<LL> :: iterator nw=s.upper_bound(c[i]);
			if (nw!=s.begin()) --nw; a[i]=(*nw); s.erase(nw); s.insert(ATK[i]);
		}
		int pd=1; rep(i,1,n) if (solve(i)==-1) { puts("-1"); pd=0; break; }; if (!pd) continue;
		printf("%lld\n",excrt(o,n));
	}
	return 0;
}