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<cstring>
#include<cstdio>
#include<cstdlib>
#include<map>
#include<cmath>
using namespace std;
#define LL long long int
LL a,b,p;
int ans;
map<LL,int> myp;
LL gcd(LL x,LL y)
{
    return (!y)?x:gcd(y,x%y);
}
LL ksm(LL x,LL y,LL mod)
{
    LL res=1;
    while(y)
    {
        if(y&1)
        res=res*x%mod;
        x=x*x%mod;
        y>>=1;
    }
    return res;
}
LL bsgs()
{
 
    if(b==1)return 0;
    int cnt=0;
    LL tmp=1,g;
    while((g=gcd(a,p))!=1)
    {
        if(b%g)return -1;//no solution
        cnt++;
        b=b/g;
        p=p/g;
        if(p==1&&b!=0)return -1;
        if(p==1&&b==0)return cnt;
        if(b==1)return cnt;
        tmp=tmp*(a/g)%p;
    }
    int m=ceil(sqrt(p));
    LL aa=1;
    myp.clear();
    for(int i=1;i<=m;i++)
    {
        aa=aa*a%p;
        myp[b*aa%p]=i;//保存最大的i
 
    }
    for(int i=1;i<=m;i++)
    {
        LL tt1=aa*tmp%p;
        if(myp.count(tt1))
            return i*m-myp[tt1]+cnt;
        tmp*=aa;
        tmp%=p;
    }
    return -1;
}
int main()
{
   while(1)
   {cin>>a>>p>>b;
    if(!(a||p||b))break;
    ans=bsgs();
    if(ans==-1)printf("No Solution\n");
    else
        printf("%d\n",ans);
   }
   return 0;
}