BSGS算法及其扩展

bsgs算法:

我们在逆元里曾经讲到过如何用殴几里得求一个同余方程的整数解。而 $ bsgs $ 就是用来求一个指数同余方程的最小整数解的:也就是对于 $ a^xequiv b mod p $ 我们可以用 $ bsgs $ 在 $ O(sqrt n) $ 的复杂度内求出关于 $ x $ 的最小正整数解。(前提是 $ p $ 为质数)

$ a^xequiv b mod p $

我们可以知道如果我们的模数p是一个质数，我们将同余式的右边以逆元的形式乘到左边来，根据殴拉定理(因为p是质数,所以a,p互质)则我们的x一定有一个小于p的正整数解!(当然,我们还要特判一下a是p倍数这一情况)。如果我们用穷举法来求这一个x最小是多少，复杂度显然不能接受，而 $ bsgs $ 就是做到了用空间来换取时间,对这个穷举进行一个改进。

其实这里 $ bsgs $ 的用空间换时间,和我们搜索用的双向搜是一个原理的。它充分利用了同余的性质，我们设m为 $ sqrt a $ 向上取整,然后原式左边就变成了 $ a^x=(a^m)^i imes a^j _{x=i+j} $ 这样我们第一步从小到大枚举i是多少,然后算出 $ (a^m)^i $ ，然后根据同余方程我们就可以知道我们还需要的 $ a^j $ 是多少。于是第二步我们就可以先预处理出所有 $ a^j $ 然后存到hash表里，然后每算一个 $ (a^m)^i $ 我们就在hash表中查看有没有对应的 $ a^j $ 存在(存hash的同时要保留 $ j $ 的值).这样我们整个算法的复杂度就是 $ O(sqrt n) $ 了!

$ code: $

inline ll bsgs(int x,ll y,ll p){
    int z=pow(p,0.5)+1;
    ll sx=1,sy=1;
    for(rg i=0;i<z;++i,sx=ksc(sx,x,p))
        ha1(i,ksc(sx,y,p));//先把所有的a^j存入hash表
    sy=1;
    for(rg i=0;i<=z;++i,sy=ksc(sy,sx,p)){//可能用到快速乘
        for(rg j=tou[sy%376637];j;j=b[j].p)//查看hash表中是否有相应的a^j
            if(sy==b[j].y)return (ll)i*z-b[j].x;
    }return p-1;//这句话其实可以省。
}

例题：

洛谷 P4884 多少个1？

$ solution $ ：

这道题还是比较裸的，我们可以讲题目中的同余方程左右同乘9再加1，得到： $ 10^xequiv k imes 9+1mod m $ 。这样就直接变成了一道板子题了（不过bsgs的题目向来是要用到快速乘的）：

#include<iostream>
#include<cstdio>
#include<iomanip>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<ctime>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>

#define lb long double
#define ll long long
#define db double
#define rg register int

using namespace std;

struct su{
    int p,x;
    ll y;
}b[376637];

int top;
ll k,m;
int tou[376637];

inline ll qr(){//快读
    char ch;
    while((ch=getchar())<'0'||ch>'9');
    ll res=ch^48;
    while((ch=getchar())>='0'&&ch<='9')
        res=res*10+(ch^48);
    return res;
}

inline ll ksc(ll x,ll y,ll p){//快速乘
    return (x*y-(ll)((lb)x/p*y)*p+p)%p;
}

inline void ha1(int x,ll y){//y用链式前向星来hash
    b[++top].x=x; b[top].y=y; y%=376637;
    b[top].p=tou[y]; tou[y]=top;
}

inline ll bsgs(int x,ll y,ll p){//快速乘
    int z=pow(p,0.5)+1;
    ll sx=1,sy=1;
    for(rg i=0;i<z;++i,sx=ksc(sx,x,p))
        ha1(i,ksc(sx,y,p));
    sy=sx;
    for(rg i=1;i<=z;++i,sy=ksc(sy,sx,p)){
        for(rg j=tou[sy%376637];j;j=b[j].p)
            if(sy==b[j].y)return (ll)i*z-b[j].x;
    }return p-1;
}

int main(){
    //freopen(".in","r",stdin);
    //freopen(".out","w",stdout);
    k=qr(); m=qr(); k=(k*9+1)%m;//稍稍处理一下
    printf("%lld
",bsgs(10,k,m));
    return 0;
}

$ BSGS $ 算法的扩展：

上面我们有讲过一般的 $ BSGS $ 算法是有前提条件的（我们的p必须是素数）。但是当我们的p不是素数的时候我们也是可以用 $ BSGS $ 求出来的，只是要先做处理。因为同余具有这样一条神奇的性质（同除性）： $ aequiv b mod m quad ->quad d=gcd(a,b,p) quad-> quad frac{a}{d} equiv frac{b}{d}mod frac{m}{d} $

根据这条性质，我们可以将a,p 取gcd 。而在扩展欧几里得解一次同余方程中我们说过，如果要有一组整数解，冲要条件就是我们的gcd（a,p）可以整除b！所以如果我们在BSGS中给三者用同除时如果b不能为整数了则无解！对这个进行特判后我们就可以一直同除（注意我们的同余式左边可以是多个a的乘积1）直到 $ a^x $ 与 $ p $ 互质。这时我们就可以再次 $ BSGS $ 来求解了（注意我们此时只保证 $ a^x $ 与 $ p $ 互质，而不是p为素数，所以可能无解！）

$ code $ ：

inline int bsgs(int x,int y,int p,int ex){
    int z=sqrt(p)+1; ll sx=1,sy;
    for(rg i=0;i<z;++i)
        ha(i,sx*y%p),sx=sx*x%p;
    sy=sx;
    for(rg i=1;i<=z;++i,sy=sy*sx%p)
        for(rg j=tou[sy%72727];j;j=b[j].p)
            if(b[j].y==sy)return i*z-b[j].x;
    return -1;//可能无解
}

例题：

洛谷 P4195 【模板】exBSGS/Spoj3105 Mod

这确确实实就是道模板，不过是多组数据，所以我们用手写hash防被卡：

#include<iostream>
#include<cstdio>
#include<iomanip>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<ctime>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>

#define lb long double
#define ll long long
#define db double
#define inf 0x7fffffff
#define rg register int

using namespace std;

struct su{
    int x,y,p;
}b[72727];

int n,m,mod,ans;
int top,tou[72727];

inline int qr(){
    char ch;
    while((ch=getchar())<'0'||ch>'9');
    int res=ch^48;
    while((ch=getchar())>='0'&&ch<='9')
        res=res*10+(ch^48);
    return res;
}

inline int gcd(int x,int y){
    int z;
    while(y){z=x,x=y,y=z%y;}
    return x;
}

inline int ksm(ll x,int y,int p){
    int res=1;
    while(y){
        if(y&1)res=res*x%p;
        x=x*x%p; y>>=1;
    }return res;
}

inline void ha(int x,int y){
    b[++top].x=x; b[top].y=y; y%=72727;
    b[top].p=tou[y]; tou[y]=top;
}

inline int bsgs(int x,int y,int p,int ex){
    int z=sqrt(p)+1; ll sx=1,sy;
    for(rg i=0;i<z;++i)
        ha(i,sx*y%p),sx=sx*x%p;
    sy=sx*ex%p;
    for(rg i=1;i<=z;++i,sy=sy*sx%p)
        for(rg j=tou[sy%72727];j;j=b[j].p)
            if(b[j].y==sy)return i*z-b[j].x;
    return -1;
}

inline int ex_bsgs(int x,int y,int p){
    if(y==1)return 0;
    ll c=1,z=x,d; int k=0;
    while(1){
        d=gcd(x,p);
        if(d==1)break;
        if(y%d)return -1;
        ++k; c=c*(x/d)%p;
        y/=d; p/=d;
        if(c==y)return k;
    }z=bsgs(x,y,p,c);
    return z==-1?z:z+k;
}

signed main(){
    //freopen(".in","r",stdin);
    //freopen(".out","w",stdout);
    while(n=qr(),mod=qr(),m=qr(),n){top=0;
        for(rg i=0;i<72727;++i)tou[i]=0;
        ans=ex_bsgs(n%mod,m%mod,mod);
        if(ans!=-1)printf("%d
",ans);
        else puts("No Solution");
    }
    return 0;
}