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

前言

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

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

可是,使用仍有限制,那就是\(x\)的系数必须为\(1\)spa

不要紧,把它再扩展一下code

题目及实现

洛谷题目传送门htm

题意分析

显然,若是咱们能干掉全部龙,那么每一次使用的剑的攻击力是已知的,设为\(k\)。那么对于每一条龙,攻击次数\(x\)必须知足\(kx\equiv a(\mod p)\);固然别忘记要先把生命值降到非正数,也就是说还要知足\(kx\geq a\)\(x\geq\lceil\frac a k\rceil\)blog

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

扩展中国剩余定理

若是有一个模线性方程组,每一个形如\(x\equiv c(\mod m)\),并且不保证\(m\)两两互质,就用不了中国剩余定理了。同步

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

咱们先把问题简化到一个方程组只有两个方程的状况,形如it

\[\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;
}
相关文章
相关标签/搜索
本站公众号
   欢迎关注本站公众号,获取更多信息