省选算法学习-BSGS与exBSGS

Posted 2021-01-29 dedicatus545

前置知识

扩展欧几里得，快速幂

都是很基础的东西

扩展欧几里得

说实话这个东西我学了好几遍都没有懂，最近终于搞明白，可以考场现推了，故放到这里来加深印象

翡蜀定理

方程$ax+by=gcd(a,b)$一定有整数解

证明：

因为$gcd(a,b)=gcd(b,a$ $mod$ $b)$

所以假设我们已经求出来了$bx+(a$ $mod$ $b)y=gcd(b,a$ $mod$ $b)$的一组整数解$(p,q)$

因为$a$ $mod$ $b=a-(lfloor frac{a}{b} floor ast b)$

所以$bp+(a-a/bast b)q=ax+by$

$b(p-(a/b)q)+aq=ax+by$

所以$x=q,y=(p-(a/b)q)$是一组合法的解

所以我们可以递归$gcd$的过程中倒着算每一层的解

当$b=0$时的解为$x=1,y=0$

BSGS

问题提出

给定$a,b,p$，求最小的$x$，使得$a^x≡b(mod$ $p)$

问题求解

显然这个东西不能直接做

考虑分块的思想

定义$m=sqrt(p)$

设$x=iast m - j$

也就是$a^{iast m}≡a^jast b(mod$ $p)$

那么我们首先把$j=0...m-1$时的$a^jast b$插入一个哈希表

然后我们枚举$i$，在哈希表里面查询$a^{iast m}$有没有出现过，如果出现过，它最大的$j$是多少

然后就可以在$O(sqrt(p))$的时间内解决这个问题了

放个板子

namespace hash{
    ll first[1000010],next[1000010],val[1000010],hash[1000010],mod=926081,cnt=0;
    void init(){memset(first,0,sizeof(first));cnt=0;}
    void insert(ll w,ll pos){
        ll p=w%mod,u;
        for(u=first[p];u;u=next[u]){
            if(hash[u]==w){val[u]=pos;return;}
            if(!next[u]) break;
        }
        if(!next[u]){
            cnt++;
            if(!first[p]) first[p]=cnt;
            else next[u]=cnt;
            val[cnt]=pos;hash[cnt]=w;next[cnt]=0;
        }
    }
    ll find(ll w){
        ll p=w%mod,u;
        for(u=first[p];u;u=next[u]){
            if(hash[u]==w) return val[u];
        }
        return -1;
    }
}
ll qpow(ll a,ll b,ll p){
    ll re=1;
    while(b){
        if(b&1) re=re*a%p;
        a=a*a%p;b>>=1;
    }
    return re;
}
ll gcd(ll a,ll b){
    if(b==0) return a;
    return gcd(b,a%b);
}
ll bsgs(ll a,ll b,ll p){
    if(b==1) return 0;
    ll i,j,m=ceil(sqrt((double)p)),tmp=b,cur,base=qpow(a,m,p);
    hash::init();
    for(j=0;j<m;j++){
        hash::insert(tmp,j);
        tmp=tmp*a%p;
    }
    tmp=1;
    for(i=m;i<=p;i+=m){
        tmp=tmp*base%p;
                cur=hash::find(tmp);
        if(~cur) return i-cur;
    }
    return -1;
}

exBSGS

使用BSGS的时候要求$gcd(a,p)=1$，扩展版的exBSGS则不需要

具体操作是这样的：

除掉公约数

假设$tmp=gcd(a,p)$

那么$(yast tmp)^x=zast tmp(mod$ $qast tmp)$

其中$y,z,q$是新设出来的量，$yast tmp=a$，$zast tmp=b$，$qast tmp=p$

这一步可以看出，如果$b$不能整除$gcd(a,p)$，那么一定无解

转化

把等式两边的含$tmp$的东西提取出来，可以得到：

$y^{tmp}=z(mod$ $q)$

然后就可以继续递归下去处理了

代码

namespace hash{
    ll first[1000010],val[1000010],hash[1000010],next[1000010],cnt=0,mod=926081;
    void init(){memset(first,0,sizeof(first));cnt=0;}
    void insert(ll w,ll pos){
        ll p=w%mod,u;
        for(u=first[p];u;u=next[u]){
            if(hash[u]==w){val[u]=pos;return;}
            if(!next[u]) break;
        }
        if(!next[u]){
            cnt++;
            if(!first[p]) first[p]=cnt;
            else next[u]=cnt;
            next[cnt]=0;val[cnt]=pos;hash[cnt]=w;
        }
    }
    ll query(ll w){
        ll p=w%mod,u;
        for(u=first[p];u;u=next[u]){
            if(hash[u]==w) return val[u];
        }
        return -1;
    }
}
ll qpow(ll a,ll b,ll p){
    ll re=1;
    while(b){
        if(b&1) re=re*a%p;
        a=a*a%p;b>>=1;
    }
    return re;
}
ll gcd(ll a,ll b){
    if(b==0) return a;
    else return gcd(b,a%b);
}
ll bsgs(ll a,ll b,ll p){
    if(b==1) return 0;//不要忘了特判
    ll i,j,tmp=1,d=1,cnt=0;
    hash::init();
    while((tmp=gcd(a,p))!=1){
        if(b%tmp) return -1;
        cnt++;b/=tmp;p/=tmp;d=d*(a/tmp)%p;//注意这个d的写法
        if(b==d) return cnt;//记得写这个
    }
    ll m=ceil(sqrt(double(p))),base=qpow(a,m,p);//注意这两个东西一定要写在这里，不要写在while上面
    tmp=b;
    for(j=0;j<m;j++){
        hash::insert(tmp,j);
        tmp=(tmp*a)%p;
    }
    for(i=m;i<=p+m;i+=m){//这里注意p+m，不然的话可能会有少数情况挂掉
        d=(d*base)%p;//同时注意这里的tmp相当于是一开始就是上面的d而不是1，也就是一开始要乘上已经除掉的东西
        tmp=hash::query(d);
        if(tmp!=-1) return i-tmp+cnt;
    }
    return -1;
}