（伪）再扩展中国剩余定理（洛谷P4774 [NOI2018]屠龙勇士）（中国剩余定理，扩展欧几里德，multiset）

前言

我们熟知的中国剩余定理，在使用条件上其实是很苛刻的，要求模线性方程组$x\equiv c(\mod m)$的模数两两互质。

于是就有了扩展中国剩余定理，其实现方法大概是通过扩展欧几里德把两个同余方程合并，具体会在下面提到。

但是，使用仍有限制，那就是$x$的系数必须为$1$。

没关系，把它再扩展一下

题目及实现

洛谷题目传送门

题意分析

显然，如果我们能干掉所有龙，那么每一次使用的剑的攻击力是已知的，设为$k$。那么对于每一条龙，攻击次数$x$必须满足$kx\equiv a(\mod p)$；当然别忘记要先把生命值降到非正数，也就是说还要满足$kx\geq a$即$x\geq\lceil\frac a k\rceil$。

可以看出不能直接用扩展中国剩余定理，但不妨碍先介绍一下它。

扩展中国剩余定理

如果有一个模线性方程组，每个形如$x\equiv c(\mod m)$，而且不保证$m$两两互质，就用不了中国剩余定理了。

扩展中国剩余定理是这样做的。

我们先把问题简化到一个方程组只有两个方程的情况，形如

$$\begin{cases}x\equiv c_1(\mod m_1)\\x\equiv c_2(\mod m_2)\end{cases}$$

把它写成不定方程的形式

$$\begin{cases}x=c_1+m_1x_1\\x=c_2+m_2x_2\end{cases}$$

这样就可以合并了，解一下$x_1$，所以$x_2$可变号

$$c_1+m_1x_1=c_2+m_2x_2\\m_1x_1+m_2x_2=c_1-c_2$$

发现变成了一个二元一次不定方程的样子，设$g=gcd(m_1,m_2)$，用扩展欧几里德求$m_1x_1+m_2x_2=g$中$x_1$的一个解$x'$，于是用$\frac{c_1-c_2}gx'+\frac{m_2}gt(t\in\mathbb Z)$可以表示$x_1$的解集（关于一般二元一次不定方程的解法和解的周期性证明可以看看蒟蒻之前写的一篇题解）

把解集带回$x=c_1+m_1x_1$得到

$$x=c_1+\frac{m_1(c_1-c_2)}gx'+\frac{m_1m_2}gt$$

$$x\equiv c_1+\frac{m_1(c_1-c_2)}gx'(\mod {\frac{m_1m_2}g})$$

这样，我们设初始方程为$x\equiv 0(\mod1)$，每次合并两个方程得到新的方程。当然中途如果有一次出现$\frac{c_1-c_2}g$不为整数则整个方程组无解。

再扩展

那么模线性方程组，每个形如$kx\equiv a(\mod p)$该怎么解好呢？

还是要合并方程。仍然设一个总方程$x\equiv c(\mod m)$，将它与当前方程合并。

下面是同步赛上手推的式子，请直接跳过这一段，因为式子是错的，只有45分。说不定有些Dalao和我的想法一样？

直接合并

$$\begin{cases}x=c+mx_0\\kx+py_0=a\end{cases}\\k(mx_0+c)+py_0=a\\kmx_0+py_0=a-kc$$

设$g=gcd(km,p)$，解不定方程得到一个解$x'$，有

$$x_0=\frac{a-kc}gx'+\frac p gt\\x=c+\frac{m(a-kc)}gx'+\frac{mp}gt$$

$$x=c+\frac{m(a-kc)}gx'(\mod\frac{mp}g)$$

为什么不能直接合并呢？因为连当前$kx\equiv a(\mod p)$能不能解都没有考虑。

所以正确的方法应该是先解当前方程$kx+py=a$。

设$g=gcd(k,p)$，解$x'$，得$x=\frac a gx'+\frac p gt$即$x\equiv\frac a gx'(\mod\frac p g)$

方程里$x$的系数被去掉了！可以从求逆元的角度来理解这一个过程。

当然，会有一些特殊情况。如果$p\mid k$且$p\mid a$，方程恒成立，我们不把它与总方程合并。如果$p\mid k$且$p\nmid a$，显然无解。

最后就是用扩展中国剩余定理合并啦。

具体实现

确定每次攻击使用的剑，直接用multiset解决，二分找到满足要求的剑（比较舒服的是upper_bound找第一个比要求大的剑，如果等于begin-iterator的话就说明没有不大于要求值的，直接选它，否则--iterator就是满足要求的），把它删掉，再加入当前奖励的剑。注意删的时候删的是iterator而不是数字。

再次提醒要注意$x\geq\lceil\frac a k\rceil$。可以求出$\max\{\lceil\frac a k\rceil\}$，如果最后总方程的$c$小于它，则要补至满足条件的最小值，用式子写一下大概是$c+m\lceil\frac{\max-c}{m} \rceil$（可以理解成，现在c和max还有差距，但是为了保证c模m的值不变，所以一次次给c加上m直到大于等于max）。

注意数据范围，会爆longlong的地方用快速乘。注意处理负数。

多组数据，注意清空和初始化。

#include<cstdio>
#include<set>
#define RG register
#define R RG LL
#define GC c=getchar()
using namespace std;
typedef long long LL;
const int N=1e5+9;
LL a[N],p[N],b[N],X,Y,G;
multiset<LL>s;
inline LL in(){
	RG char GC;
	while(c<'-')GC;
	R x=c&15;GC;
	while(c>'-')x*=10,x+=c&15,GC;
	return x;
}
void exgcd(R a,R b){
	if(!b){X=1;Y=0;G=a;return;}
	exgcd(b,a%b);
	(Y^=X^=Y^=X)-=a/b*X;
}
inline LL mul(R b,R k,R m){//快速乘
	R a=0;
	for(;k;k>>=1,b=(b<<1)%m)
		if(k&1)a=(a+b)%m;
	return a;
}
int main(){
	freopen("dragon.in","r",stdin);
	freopen("dragon.out","w",stdout);
	R T=in(),n,m,i,c,k,mx;
	RG multiset<LL>::iterator it;
	E:while(T--){
		n=in();m=in();
		for(i=1;i<=n;++i)a[i]=in();
		for(i=1;i<=n;++i)p[i]=in();
		for(i=1;i<=n;++i)b[i]=in();
		s.clear();//注意清空
		for(i=1;i<=m;++i)s.insert(in());
		mx=c=0;m=1;//初始总方程
		for(i=1;i<=n;++i){
			it=s.upper_bound(a[i]);//谨慎选择lower_bound和upper_bound
			if(it!=s.begin())--it;
			k=*it;s.erase(it);s.insert(b[i]);//小心手一滑erase(*it)(居然还有90分)
			mx=max(mx,(a[i]-1)/k+1);//更新限制
			k%=p[i];a[i]%=p[i];//开始解方程，去掉系数
			if(!k&&a[i]){puts("-1");goto E;}
			if(!k&&!a[i])continue;//这两个要特判
			exgcd(k,p[i]);
			if(a[i]%G){puts("-1");goto E;}
			p[i]/=G;
			a[i]=mul(a[i]/G,(X%p[i]+p[i])%p[i],p[i]);
			exgcd(m,p[i]);//开始合并，X和a-c都可能是负数
			if((a[i]-c)%G){puts("-1");goto E;}
			m=m/G*p[i];
			c=(c+mul(mul(m/p[i],((a[i]-c)%m+m)%m,m),(X%m+m)%m,m))%m;
		}
		printf("%lld\n",c>=mx?c:c+m*((mx-c-1)/m+1));//满足限制
	}
	return 0;
}