luogu P4774 [NOI2018] 屠龙勇士 - 数论 - exgcd - 扩展crt

传送门:https://www.luogu.org/problemnew/show/P4774

因为没a掉此题而cu滚粗的菜鸡来报个到。

首先我们可以显然地发现对于每条龙,攻击它用到的剑是唯一确定的,用一个multiset预处理即可,记作di。

需要强调两点:第一是不要用成set,第二是虽然里面存的元素是int但是你会用ll去查询所以还是要开ll,没错我的25分就是这么丢的……

首先pi=1的点很好处理,因为此时只要打到血量<=0就死定了,所以答案就是max((ai / di)上取整)

然后我们考虑对于最后的答案x,一定满足如下限制:

对于任意i,均有 di * x % pi == ai

也就是n个同余方程的形式。

我们可以对于每个同余方程,用exgcd解出一个x的通解。

注意得到x的最小非负整数解qi后,加上任意整数倍的ki =(pi / gcd(di,pi))也一定是一个合法解

因此我们的同余方程可以化成如下形式:

对于任意i,满足x % pi == qi

然后合并这n个同余方程组就是经典的扩展crt啦

等等……我不会扩展crt!

于是我在考场上推了两个多小时的式子,最后得出了一个玄学做法:

假设我们有两个同余方程:

x % a1 == b1

x % a2 == b2

我们要将其合并为

x % lcm(a1,a2) == b3

b3即为我们的所求。

我们考虑写成不定方程的形式

x == b1 + k1 * a1

x == b2 + k2 * a2

先“钦定”k1 = 0,看一下带入第二个不定方程的结果与a1相差多少,记作q1。

然后考虑k1++,会使得第二个不定方程的结果变化多少,记作q2。

于是我们就有

k1 * q2 % a2 == q1

可以exgcd解出k1,然后就可以还原回b3的值啦。

(并不知道正经的扩展crt是怎么写的,也许跟这个思路类似?)

如果上面任意一个exgcd或扩展crt的步骤无解,则直接返回无解即可。

于是这题就做完了,复杂度O(n log n),是不是很休闲?(逃)

哦对了,由于1e12直接相乘会爆ll,因此需要写快速乘。 

由于是自己一边推一边写的,代码极丑请原谅~

#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<ctime>
#include<set>
using namespace std;
#define li long long
#define gc getchar()
#define pc putchar
li read(){
	li x = 0,c = gc;
	while(!isdigit(c)) c = gc;
	while(isdigit(c)){
		x = x * 10 + (c ^ '0');
		c = gc;
	}
	return x;
}
void print(li x){
	if(x >= 10) print(x / 10);
	pc(x % 10 + '0');
}
int t,n,m;
li a[100010],b[100010],c[100010];
int d[100010];
multiset<li> s;
li gcd(li x,li y){
	return !y ? x : gcd(y,x % y);
}
li lcm(li x,li y){
	return x / gcd(x,y) * y;
}
li exgcd(li a,li b,li &x,li &y){
	if(!b){
		x = 1;y = 0;return a;
	}
	li p = exgcd(b,a % b,x,y);
	li t = x;x = y;y = t - a / b * y;
	return p;
}
li ksc(li q,li w,li p){
	q %= p;w %= p;
	if(q < w) swap(q,w);
	if(p <= 10000000000000ll){
		li as = 0;
		while(w){
			(as += w % 100000 * q) %= p;
			q = q * 100000 % p;
			w /= 100000;
		}
		return as;
	}
	li as = 0;
	while(w){
		if(w & 1) (as += q) %= p;
		(q += q) %= p;
		w >>= 1;
	}
	return as;
}
int main(){
	int i,j,k,l;
	li x,y,p,q,r,a1,a2,b1,b2,c1,c2,tp;
	t = read();
	while(t--){
		s.clear();
		li mx = 0;
		n = read();m = read();
		for(i = 1;i <= n;++i) a[i] = read();
		for(i = 1;i <= n;++i) b[i] = read(),mx = max(mx,b[i]);
		for(i = 1;i <= n;++i) c[i] = read();
		for(i = 1;i <= m;++i){
			j = read();
			s.insert(j);
		}
		for(i = 1;i <= n;++i){
			multiset<li>::iterator it = s.upper_bound(a[i]);
			if(it != s.begin()) --it;
			d[i] = *it;
			s.erase(it);
			s.insert(c[i]); 
		}
		if(mx == 1){
			li as = 0;
			for(i = 1;i <= n;++i) as = max(as,(a[i] - 1) / d[i] + 1);
			print(as);pc('\n');continue;
		}
		p = exgcd(d[1],b[1],x,y);
		if(a[1] % p){
			puts("-1");continue;
		}
		q = b[1] / p;r = a[1] / p % q;
		if(x < 0) x += q;
		b1 = q;b2 = ksc(r,x,q);
		bool fg = 0;
		for(i = 2;i <= n;++i){
			p = exgcd(d[i],b[i],x,y);
			if(a[i] % p){
				puts("-1");fg = 1;break;
			}
			q = b[i] / p;r = a[i] / p % q;
			if(x < 0) x += q;
			a1 = q;a2 = ksc(r,x,q);
			tp = ((b2 - a2) % b1 + b1) % b1;
			p = exgcd(a1 % b1,b1,x,y);
			if(tp % p || a2 % p != b2 % p){
				fg = 1;puts("-1");break;
			}
			c1 = a1 / p * b1;
			q = b1 / p;r = tp / p % q;
			if(x < 0) x += q;
			x = ksc(r,x,q);
			c2 = (a2 + ksc(x,a1,c1)) % c1;
			b1 = c1;b2 = c2;
		}
		if(fg) continue;
		print(!b2 ? b1 : b2);pc('\n');
	}
	return 0;
}

 

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值