【算法】BSGS算法的推导与实现

BSGS简介

BSGS算法，全称Baby Step Giant Step算法，用于求解关于$x$的形如$a^x\equiv b(\mathrm{m}\mathrm{o}\mathrm{d}p)$ ，$p$为质数的方程。

求解过程

不妨设$x=im-j$，其中$m=\lceil \sqrt{p}\rceil$，$j\in [0,m)$，$i\in [1,m]$。于是原方程转化为:
$$a^{im-j}\equiv b(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
继续转化:
$$a^{im}\equiv a^{j}b(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
于是我们可以枚举$j$，并将$a^{j}b\mathrm{m}\mathrm{o}\mathrm{d}p$的答案记录在map中（当然如果你是手写哈希表的巨佬，那就当我没说）。map中的key存的是答案，value存的是$j$。

然后我们再枚举$i$，在map中查找$a^{im}\mathrm{m}\mathrm{o}\mathrm{d}p$,如果查询到了，那么查询到的即为$j$，那么$x=im-j$即为答案。

---

一些证明

注意到我们求解时取了$m=\lceil \sqrt{p}\rceil$，这能保证我们枚举的时间复杂度均为$O(\sqrt{p})$，但是却限定了答案$x\in [0,p]$，那么如何保证解一定在这个范围内呢？

证明过程需要费马小定理：$a^{p-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$，其中$p$为质数，$gcd(a,p)=1$。

引理

$$a^{k\mathrm{m}\mathrm{o}\mathrm{d}(p-1)}\equiv a^k(\mathrm{m}\mathrm{o}\mathrm{d}p)$$

条件与费马小定理相同。

证明：我们可以把$k\mathrm{m}\mathrm{o}\mathrm{d}(p-1)$看做$k-n(p-1)$，那么原方程化为：
$$\frac{a^k}{(a^{p-1}{)}^n}\equiv a^k(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
通过费马小定理，可知$a^{p-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$，所以原式显然成立。

正式的证明

其实知道引理后直接一个显然就好啦

由引理可知，即使在$(p,+\infty )$中存在解，我们依旧可以通过$\mathrm{m}\mathrm{o}\mathrm{d}(p-1)$操作在$ [0,p]$内找到解。

---

代码

-1表示无解

ll bsgs(ll a, ll b, ll p){
    a %= p, b %= p;
    if(!a && !b) return 1;
    if(!a || !b) return -1;
    mp.clear();		//初始化map
    ll m = ceil(sqrt(p*1.0)), tmp = 1;
    mp[b] = 0;
    for(int j = 1; j <= m; j++){			//枚举a^j*b
        tmp = tmp*a % p;
        if(!mp[tmp*b%p])
            mp[tmp*b%p] = j;
    }//循环完成后tmp=a^m
    ll t = 1, ans;
    for(int i = 1; i <= m; i++){		//枚举a^(im)
        t = t*tmp%p;
        if(mp[t]){			//判断是否存在
            ans = i*m-mp[t];
            return (ans%p+p)%p;
        }
    }
    return -1;
}

---

例题（水经验）

洛谷2485 [SDOI2011]计算器

前两个操作是快速幂和扩欧的模板，第三个操作是BSGS模板。

代码

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

int n, t;

ll qpow(ll a, ll n, ll mod){
    ll res = 1, base = a;
    while(n){
        if(n&1) res = (res*base)%mod;
        base = (base*base)%mod;
        n >>= 1;
    }
    return res;
}

ll exgcd(ll a, ll b, ll &x, ll &y){
    if(!b) {
        x = 1, y = 0;
        return a;
    }
    ll res = exgcd(b, a%b, x, y);
    ll t = x;
    x = y;
    y = t-a/b*y;
    return res;
}

map<ll, ll> mp;

ll bsgs(ll a, ll b, ll p){
    a %= p, b %= p;
    if(!a && !b) return 1;
    if(!a || !b) return -1;
    mp.clear();
    ll m = ceil(sqrt(p*1.0)), tmp = 1;
    mp[b] = 0;
    for(int j = 1; j <= m; j++){
        tmp = tmp*a % p;
        if(!mp[tmp*b%p])
            mp[tmp*b%p] = j;
    }
    ll t = 1, ans;
    for(int i = 1; i <= m; i++){
        t = t*tmp%p;
        if(mp[t]){
            ans = i*m-mp[t];
            return (ans%p+p)%p;
        }
    }
    return -1;
}

int main()
{
    cin >> n >> t;
    ll y, p, z;
    for(int i = 1; i <= n; i++){
        scanf("%lld%lld%lld", &y, &z, &p);
        if(t == 1){
            printf("%lld\n", qpow(y, z, p));
        }
        else if(t == 2){
            ll x, t;
            ll g = exgcd(y, p, x, t);
            if(z % g){
                puts("Orz, I cannot find x!");
                continue;
            }
            y /= g, z /= g, p /= g;
            x = (x*z%p+p)%p;
            printf("%lld\n", x);
        }
        else{
            ll ans = bsgs(y, z, p);
            if(ans == -1){
                puts("Orz, I cannot find x!");
            }
            else{
                printf("%lld\n", ans);
            }
        }
    }
    
    return 0;
}