BSGS

· · 算法·理论

Baby-Step Giant-Step，简称 BSGS。是一种在 O(\sqrt{\text{ord}_{p}(a)}) 的时间复杂度内求解离散对数问题的算法。本算法属于数论板块，若无特殊说明均在正整数范围内讨论。

处理对象

求解该同余方程的最小解，或者报告无解。保证 a\perp p。

a^l\equiv b\pmod p

算法内容

前置内容

定义离散函数 f(x)=a^x，容易发现该函数在模 p 意义下存在最小正周期 N。对于 a，该最小正周期可以写作 \text{ord}_{p}(a)，也称作 a 在模 p 意义下的阶。若 a^l\equiv b\pmod p，那么 a^{l+N}\equiv a^l\equiv b\pmod p，所以最小 l 一定属于 [0,N)。

::::info[如何求解 N]{open}

考虑如何求这个 a 在模 p 意义下的阶 N。根据欧拉定理 a^{\varphi(p)}\equiv a^0\pmod p，所以 \varphi(p) 是 f(x) 的一个函数正周期。对于 N，一定有 N\mid\varphi(p)。因为实际的 N 是 \varphi(p) 的约数，所以考虑从 \varphi(p) 中一步一步剥离不必要的因数。

记 N=\varphi(p)，然后枚举 \varphi(p) 的质因数 P_i，若 a^{\frac{N}{P_i}}\equiv a^0\pmod p。说明去掉一个 P_i，N 依旧是 f(x) 的函数正周期；对于每个 P_i 重复尝试去除直到无法去除为止。用试除法对 \varphi(p) 进行质因数分解然后暴力考虑剥离每个质因数可以做到 O(\frac{\sqrt{p}}{\log p}+\log^2p)。

如果上述同余方程存在解，那么一定属于 [0,N)。如果 a 为原根，那么 N=p-1，暴力枚举复杂度直接就炸了。考虑一种快速的算法。

BSGS 算法

对同余方程进行变换。令 l=xB-y

a^{xB-y}\equiv b\pmod p

::::info[对于上式的一些思考]{open}

该内容对理解 BSGS 很有帮助，建议阅读。

考虑一种跳步的思想。

将 l 的选取放到数轴上，假如最小答案的位置是绿色箭头所示。暴力做法会从左端点一步一步枚举到答案。将 l 拆解为 x 步向右的大步，y 步向左的小步。若我们在某个大步跳后知道已经越过答案，即使我们不知道具体是什么，但是我们暴力往回走一步一步找到答案。在 B=\sqrt{N} 的时候最优。

但有一个问题，我们怎么知道是否已经越过答案了。不好实现，于是考虑从答案下手进行处理。那么我们再对式子进行变换。

::::

由于 a\perp p，所以 a^{y}\perp p。两边同乘 a^{y} 得到下式。

a^{xB}\equiv ba^y\pmod p

发现将原先拓展大步 x 和回退小步 y 分开了。预处理大步可能到的状态，然后再考虑小步可能回退的状态（就是暴力枚举 x,y）。在枚举 y 的时候，如果在之前预处理的时候存在一个 x 满足 a^{xB} 与当前的 ba^y 同余。说明就找到了一组解。

枚举 x 预处理的时候用哈希或map维护 a^{xB} 的值和对应的 xB。当 y 找到满足 a^{xB}\equiv ba^y\pmod p 的 a^{xB} 时。直接找到对应的 xB，xB-y 就是 l 的一种解。

左边我们预处理 a^{xB}，枚举 x 并哈希的复杂度是 O(\frac{N}{B}\log V) 的。右边枚举 y 是 O(B\log V) 的。整体复杂度 O(\log V(B+\frac{N}{B})) 在 B=\sqrt{N} 的时候最优为 O(\log V\sqrt{N})。最坏情况 N 和 p 同阶复杂度 O(\log V\sqrt{p})。

实现

::::success[核心代码]

int T=ceil(sqrt(mod));//直接取模数作为上界
for(int i=0;i<=T;i++)
{
    if(mp[mypow(A,T*i)]!=0)break;//出现相同说明进入新的函数周期 break 掉即可，某些时候可以优化常数。
    mp[mypow(A,T*i)]=T*i;
}
for(int i=0;i<=T;i++)
{
    int cnt=(B*mypow(A,i))%mod;
    if(mp[cnt])ans=min(ans,mp[cnt]-i);
}

::::

例题

P3846 [TJOI2007] 可爱的质数

这是一道 BSGS 的裸题直接实现上述算法即可。

在实现 BSGS 的时候，求出 a 在模 p 意义下的阶 N 其实是不必要的，但会出现一些问题。这篇帖子指出，如果不求解 N。直接在 [0,p) 的范围内求，在正确性上是没有任何问题的。但要注意的是因为 p 并不是最小正周期，甚至连周期都不是。所以可能会出现在模 p 意义下 a^{xB} 相同的 x 多次出现。因为属于不同的周期，但因为要求的是最小 l。所以当重复出现 a^{xB} 时，说明进入了一个新的函数周期。直接break掉即可。

::::success[Accepted code]

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int V=1LL<<31;
int A,B,mod;
unordered_map<int,int>mp;
int mypow(int a,int k)
{
    if(k==0)return 1;
    int res=a;k--;
    while(k)
    {
        if(k&1)res=(res*a)%mod;
        a=(a*a)%mod;
        k>>=1;
    }
    return res;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    cin>>mod>>A>>B;
    int T=ceil(sqrt(mod));//直接取模数作为上界
    for(int i=0;i<=T;i++)
    {
        if(mp[mypow(A,T*i)]!=0)break;//出现相同说明进入新的函数周期 break 掉即可，某些时候可以优化常数。
        mp[mypow(A,T*i)]=T*i;
    }
    int ans=1e18;
    for(int i=0;i<=T;i++)
    {
        int cnt=(B*mypow(A,i))%mod;
        if(mp[cnt])ans=min(ans,mp[cnt]-i);
    }
    if(ans==1e18)cout<<"no solution";
    else cout<<ans;
    return 0;
}
/*
3 2 2
*/

::::

P2485 [SDOI2011] 计算器

还是模板题，操作 2 直接逆元算就可以了。

注意操作 3 中 y 和 p 可能不满足互质，之前的 BSGS 代码需要特判。

::::success[Accepted code]

#include<bits/stdc++.h>
using namespace std;
#define int long long
int mod;

int mypow(int a,int k)
{
    if(k==0)return 1;
    int res=a;k--;
    while(k)
    {
        if(k&1)res=(res*a)%mod;
        a=(a*a)%mod;
        k>>=1;
    }
    return res;
}
unordered_map<int,int>mp;
int BSGS(int A,int B)
{
    mp.clear();
    if(B==1)
    {
        return 0;
    }
    int T=ceil(sqrt(mod));
    for(int i=1;i<=T;i++)
    {
        if(mp[mypow(A,T*i)]!=0)mp[mypow(A,T*i)]=min(mp[mypow(A,T*i)],T*i);
        else mp[mypow(A,T*i)]=T*i;
    }
    int ans=1e18;
    for(int i=0;i<=T;i++)
    {
        int cnt=(B*mypow(A,i))%mod;
        if(mp[cnt])ans=min(ans,mp[cnt]-i);
    }
    if(ans==1e18)return -1;
    return ans;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    int T,K;cin>>T>>K;
    while(T--)
    {
        int y,z;cin>>y>>z>>mod;
        if(K==1)
        {
            cout<<mypow(y,z)<<"\n";
            continue;
        }
        if(K==2)
        {
            if(z%__gcd(y,mod))cout<<"Orz, I cannot find x!\n";
            else cout<<(z*mypow(y,mod-2))%mod<<"\n";
            continue;
        }
        if(K==3)
        {
            if(y%mod==1&&z%mod==1)cout<<"0\n";
            else if(y%mod==z%mod)cout<<"1\n";
            else if(y%mod==0)
            {
                if(z%mod==1)cout<<"0\n";
                else cout<<"Orz, I cannot find x!\n";
            }
            else
            {
                int res=BSGS(y,z);
                if(res==-1)cout<<"Orz, I cannot find x!\n";
                else cout<<res<<"\n";
            }
        }
    }
    return 0;
}

::::

P3306 [SDOI2013] 随机数生成器

结合了一定拆式子和特判的 BSGS 题。

给出 x_1,a,b,p,t 和 x_i 的递推式 x_{i+1}\equiv a\times x_i+b\pmod p。求最小 l 使 x_l=t。

考虑拆出通项公式。

::::info[拆式子]

x_{i+1}\equiv a\times x_i+b\pmod p

手动推一下前 4 项。

\begin{aligned} x_1&\equiv x_1\ &(\bmod p)\\ x_2&\equiv a^1\times x_1+a^0\times b\ &(\bmod p)\\ x_3&\equiv a^2\times x_1+a^0\times b+a^1\times b\ &(\bmod p)\\ x_4&\equiv a^3\times x_1+a^0\times b+a^1\times b+a^2\times b\ &(\bmod p)\\ \end{aligned}

容易发现。

x_i\equiv a^{i-1}x_1+b\times \sum_{j=0}^{i-2}a^j\pmod p

等比数列求和得到通项公式，注意 a>1 所以 a=0,a=1 的情况后文会提到如何特判。

x_i\equiv a^{i-1}\times x_1+b\times \frac{1-a^{i-1}}{1-a}\pmod p

::::

拆完式子得到。

x_i\equiv a^{i-1}\times x_1+b\times \frac{1-a^{i-1}}{1-a}\pmod p

直接代 t 得到需要求解的同余方程。

t\equiv x_l\equiv a^{l-1}\times x_1+b\times \frac{1-a^{l-1}}{1-a}\pmod p

题目 p 是质数故直接当普通方程移项化简得到下式。

a^{l-1}\equiv \frac{a\times t-t+b}{x_1\times a-x_1+b}\pmod p

::::warning[特殊情况]{open} 接下来就是几个特判要处理。

第一天比较特殊，直接特判是否 x_1=t。

由于使用的等比数列，a\le 1 的情况需要特殊处理。

对于 t\equiv x_1+b\times(l-1)\pmod p，化简得到 t-x_1\equiv b\times(l-1)\pmod p。

根据裴蜀定理，该方程存在解。当且仅当 \gcd(p,b)|(t-x_1)，如果存在解直接逆元计算即可。

if(a==1)
{
    t=(t-x+mod)%mod;
    if(t%__gcd(b,mod))cout<<"-1\n";
    else
    {
        if((t*mypow(b,mod-2)+1)%mod==0)cout<<mod<<"\n";//不可能是第 0 天，最先的一天为第 p 天
        else cout<<(t*mypow(b,mod-2)+1)%mod<<"\n";
    }
}

:::: ::::success[Accepted code]

#include<bits/stdc++.h>
using namespace std;
#define int long long
int mod,a,b,x,t;
int mypow(int a,int k)
{
    if(k==0)return 1;
    int res=a;k--;
    while(k)
    {
        if(k&1)res=(res*a)%mod;
        a=(a*a)%mod;
        k>>=1;
    }
    return res;
}
unordered_map<int,int>mp;
int BSGS(int A,int B)
{
    mp.clear();
    int T=ceil(sqrt(mod));
    for(int i=0;i<=T;i++)
    {
        if(mp[mypow(A,T*i)]!=0)break;
        else mp[mypow(A,T*i)]=T*i;
    }
    int ans=1e18;
    for(int i=0;i<=T;i++)
    {
        int cnt=(B*mypow(A,i))%mod;
        if(mp[cnt])ans=min(ans,mp[cnt]-i);
    }
    if(ans==1e18)return -1;
    return ans;
}
void solve()
{
    cin>>mod>>a>>b>>x>>t;
    if(x==t)
    {
        cout<<1<<"\n";
        return ;
    }
    if(a==0)
    {
        if(b%mod!=t)cout<<"-1\n";
        else cout<<"2\n";
        return ;
    }
    if(a==1)
    {
        t=(t-x+mod)%mod;
        if(t%__gcd(b,mod))
        {
            cout<<"-1\n";
        }
        else
        {
            if((t*mypow(b,mod-2)+1)%mod==0)cout<<mod<<"\n";
            else cout<<(t*mypow(b,mod-2)+1)%mod<<"\n";
        }
        return ;
    }
    int B=(((a*t-t+b)%mod)*mypow((x*a-x+b)%mod,mod-2))%mod;
    int res=BSGS(a,B);if(res!=-1)res++;//计算出来的是 l-1
    cout<<res<<"\n";
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    int T;cin>>T;
    while(T--)solve(); 
    return 0;
}

::::