BSGS算法

大步小步算法 BSGS(baby steps giant steps)

BSGS是用来解决离散对数问题的，即$a^x\equiv b(\mathrm{m}\mathrm{o}\mathrm{d}p)$。其中，a,b,p已知，且a和p互质，求x。根据欧拉定理，我们知道$a^{\varphi (p)}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$,而$\varphi (p)<p$,所以，采用枚举法，我们可以在$O(p)$时间复杂度求出x。
BSGS可以在$O(\sqrt{p})$的时间复杂度内求出x，它是这么做的：
首先，令 $$m=\lceil \sqrt{p}\rceil ,r=x\mathrm{m}\mathrm{o}\mathrm{d}m$$则 $$x=k\ast m+r\text{\hspace{1em}(1)}$$
其中 $0\le k<m，0\le r<m$
于是有$$a^{k\ast m+r}\equiv b(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
因为a与p互质，所以两边同时乘以$a^{-r}$,则
$$a^{k\ast m}\equiv b\ast a^{-r}(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
对于$(0\le i\le r)$,求出所有的$b\ast a^{-i}\mathrm{m}\mathrm{o}\mathrm{d}p$，将值和对应的i都保存在map中.
然后再求左边$a^{j\ast m}\mathrm{m}\mathrm{o}\mathrm{d}p(0\le j<k)$ ,并在map中查找是否出现过相同的值。如果有，则已经找出答案了,则答案为$j\ast m+i$。如果没有，表示无解。
其实我们可以稍微改进一下，以避免求逆元。
将公式（1）变换一下：
$$x=k\ast m-r\text{\hspace{1em}(2)}$$
此时，$1\le k\le m,0<r\le m$
则$$a^{k\ast m-r}\equiv b(\mathrm{m}\mathrm{o}\mathrm{d}p)\text{\hspace{1em}(3)}$$
两边同时乘以$a^r$,可得：
$$a^{k\ast m}\equiv b\ast a^r(\mathrm{m}\mathrm{o}\mathrm{d}p)\text{\hspace{1em}(4)}$$
先求出右边所有的 $b\ast a^i(\mathrm{m}\mathrm{o}\mathrm{d}p)(1\le i\le r)$ ，保存在map中，然后再求左边的$a^{j\ast m}\mathrm{m}\mathrm{o}\mathrm{d}p$,并在map中查找是否出现过。如果出现过，左边枚举的是j，右边枚举的是i，则答案为$x=j\ast m-i$.这样就可以不用求逆元了，但仍然用到了逆元。因为我们的推导必须是等价推导，只有当a与p互质时，即$a^r$的逆元存在时，式（3）和式（4）才是等价的。因为如果a不与p互质，式（4）是不能推出式（3）的。

例题1：（poj2417）
题意：给出一个质数P，一个整数B，2<=B<P,和一个整数N，1<=N<=P,求底为B，模为P下的离散对数,如果有多个解，则输出最小的非负整数解。即解方程$B^x\equiv N\mathrm{m}\mathrm{o}\mathrm{d}P$，求最小的x，x>=0

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<map>
#include<cmath>
using namespace std;
#define LL long long int
map<LL,int> mp;
int b,n,p,x,ans;
int bsgs()
{
    mp.clear();
    int m=ceil(sqrt(p));
    LL tmp=1;
    for(int i=1;i<=m;i++)
    {
      tmp=tmp*b%p;
      mp[tmp*n%p]=i; //这里会覆盖掉之前的解，即保存的是最大的i。
    }
    int res=1;
    for(int i=1;i<=m;i++)
    {  res=res*tmp%p;
       if(mp[res])
       {
           ans=m*i-mp[res]; //因为mp[res]保存的是最大的，所以ans是最小的解。
           return ans;
       }
    }
    return -1;
}
int main()
{
    while(scanf("%d%d%d",&p,&b,&n)!=-1)
    {
        n%=p,b%=p;
        if(n==1)
            ans=0;
        else if(b==0)
                {if(n!=0)
                   ans=-1;
                 else ans=1;
                }
        else
        ans=bsgs();
        if(ans>=0)printf("%d\n",ans);
        else printf("no solution\n");
    }
    return 0;
}

扩展BSGS

对于$$a^x\equiv b(\mathrm{m}\mathrm{o}\mathrm{d}p)\text{\hspace{1em}(1)}$$

如果a与p不互质，则不能直接套用BSGS。而扩展BSGS正是解决这种情况的。
可以将（1）变形：
$$a^x+kp=b\text{\hspace{1em}(2)}$$
设a与p的最大公约数为g，即$g=\left(a,p\right)$.若$g\nmid b$,则方程无解，不过这里有一个例外，即当b=1时，x=0。
将式(2)左右两边除以g，得到：
$$a^{x-1}\frac{a}{g}+k\frac{p}{g}=\frac{b}{g}\text{\hspace{1em}(3)}$$
设$a^{\prime }=\frac{a}{g},p^{\prime }=\frac{p}{g}$,$b^{\prime }=\frac{b}{g}$,
$$a^{x-1}{a}^{\prime }+k{p}^{\prime }=b^{\prime }$$
若a与$\frac{p}{g}$仍然不互质，则继续找出a与p’的最大公约数g’,(3)两边继续除以g’,直到最后a与p’互质为止。
这里我们假设a与p’互质，此时，a’与p’也是互质的，于是可以转换为：

$$a^{x-cnt}\equiv b^{\prime }(a^{\prime }{)}^{-1}(\mathrm{m}\mathrm{o}\mathrm{d}{p}^{\prime })$$
其中cnt表示两边除以最大公约数g的次数。
此处右边有逆元。为了避免求逆元，可以将a’保留在左边。在用bsgs枚举左边时，初始值就设为a’。
如果在除以g的过程中，发现$b^{\prime }(a^{\prime }{)}^{-1}=1$，则立马可以得到答案，即x-cnt=0，即x=cnt。
此时，直接套用基础的BSGS即可了。
不要忘了最后的答案要加上cnt。

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cstring>
#include<map>
#include<cmath>
#include<vector>
using namespace std;
#define LL long long int
#define HASHL 10007
int x,z,k,m,ans;
int hashcode(int a)
{
    a=(a*998244353LL+100007)%HASHL;
    return a%HASHL;
}
vector<pair<int,int> > myv[HASHL+5];
int getphi(int x)
{
    int res=x;
    for(int i=2;i*i<=x;i++)
    {
        if(x%i==0)
        {
            while(x%i==0)
            x/=i;
        }
        res=res*(i-1)/i;
    }
    if(x>1)
    res=res*(x-1)/x;
    return res;
}
bool flg[HASHL+5];
int gcd(int a,int b)
{
    if(!b)return a;
    else return gcd(b,a%b);
}
int exbsgs()
{
    memset(flg,0,sizeof flg);
    if(z==0)return -1;
    if(z==1)
    {
        if(k==0)return 0;
        else return -1;
    }
    if(k==1)return 0;
    x%=z,k%=z;
    int g=1,cnt=0,tmpk=1;
    while((g=gcd(x,z))>1)
    {
        if(tmpk%z==k%z)return cnt;
        if(k%g)return -1;
        k=k/g;
        z=z/g;
        tmpk=1ll*tmpk*(x/g)%z;
        cnt++;
    }
    int m=ceil(sqrt(double(z)));
    int tmp=1,tmpid=0;
    for(int i=1;i<=m;i++)
    {
        tmp=1ll*tmp*x%z;
        tmpid=hashcode(1ll*tmp*k%z);
        if(flg[tmpid]==0)
        {
            myv[tmpid].clear();
            flg[tmpid]=1;
        }
        myv[tmpid].push_back(make_pair(1ll*tmp*k%z,i));
    }
    int res=tmpk%z,tmpres=1;
    for(int i=1;i<=m;i++)
    {
        res=1ll*res*tmp%z;
        tmpres=hashcode(res);
        if(flg[tmpres]==1)
        {
            for(int j=myv[tmpres].size()-1;j>=0;j--)
            {
                if(myv[tmpres][j].first==res)
                {
                    return (i*m-myv[tmpres][j].second)+cnt;
                }
                
            }
        }
    }
    return -1;
}
int main()
{
    int ans;
    while(scanf("%d %d %d",&x,&z,&k)!=-1)
    {
        if(!(x||z||k))return 0;
        ans=exbsgs();
        if(ans==-1)printf("No Solution\n");
        else printf("%d\n",ans);
    }
    return 0;
}