前言
我們熟知的中國剩余定理,在使用條件上其實是很苛刻的,要求模線性方程組\(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\)兩兩互質,就用不了中國剩余定理了。
擴展中國剩余定理是這樣做的。
我們先把問題簡化到一個方程組只有兩個方程的情況,形如
把它寫成不定方程的形式
這樣就可以合並了,解一下\(x_1\),所以\(x_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\equiv 0(\mod1)\),每次合並兩個方程得到新的方程。當然中途如果有一次出現\(\frac{c_1-c_2}g\)不為整數則整個方程組無解。
再擴展
那么模線性方程組,每個形如\(kx\equiv a(\mod p)\)該怎么解好呢?
還是要合並方程。仍然設一個總方程\(x\equiv c(\mod m)\),將它與當前方程合並。
下面是同步賽上手推的式子,請直接跳過這一段,因為式子是錯的,只有45分。說不定有些Dalao和我的想法一樣?
直接合並
設\(g=gcd(km,p)\),解不定方程得到一個解\(x'\),有
為什么不能直接合並呢?因為連當前\(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;
}
