【算法笔记】Baby Step Giant Step（BSGS）及其扩展
- 本文总计约 8000 字，阅读大约需要 30 分钟。
- 警告！警告！警告！本文有大量的 \(\LaTeX\) 公式渲染，可能会导致加载异常缓慢！
大手牵小手，大步小步走。
前言
BSGS，虽然在洛谷上是蓝题的水平，但实际上的难度远远不及蓝题。
然而，BSGS 在解决离散对数问题上非常有用，在历年省选和 NOI 中，考察 BSGS 的问题都非常的多。所以它的重要性又不止是蓝题。
更进一步地，BSGS 背后的思想— —分块思想，也是非常重要的思想了，所以，学习 BSGS，是一件很有用的事。
题目引入
给定 \(a,b,m\)，是否存在一个正整数 \(x\)，使得 \(a^x\equiv b\pmod m\)？若存在，求出 \(b\) 的最小值。
例如：\(a=2\)，\(b=3\)，\(p=5\)，那么显然，当 \(x=3\) 时，有 \(2^3\equiv 8 \equiv 3\pmod{5}\)。
再例如：\(a=2\)，\(b=3\)，\(p=10\) 时，不存在这样的 \(x\)。
\(\text{Subtask 1}\)：\(a,b\le 10^9\)，\(p\le 10^3\)；
\(\text{Subtask 2}\)：\(a,b\le 10^9\)，\(p\le 10^9\)，保证 \(a\perp p\)；
\(\text{Subtask 2}\)：\(a,b\le 10^9\)，\(p\le 10^9\)。
上面的这个问题，就是非常经典的离散对数问题（Discrete Logarithm），这个问题在密码学中起了很大的贡献。
暴力算法
依照惯例，我们应该问一下：暴力怎么做？
当然，我们可以暴力枚举 \(x\)，计算 \(a^x\bmod p\) 的值。根据 Euler 定理，这个最小的 \(x\) 应该是 \(\mathcal{O}(\varphi(p))\) 级别的。更具体地，应该有 \(x\le 2\varphi(p)<2p\)。
这也就是说，我们可以通过枚举 \(1\sim 2p\) 中的每一个 \(x\)，计算 \(a^x\equiv b\pmod p\) 是否成立，这个过程的时间复杂度是 \(\mathcal{O}(p)\) 的。
代码很简洁，就这么写：
#include <iostream>
using namespace std;
typedef long long ll;
ll a, p, b;
ll ans = 1;
int main(void) {
scanf("%lld%lld%lld", &p, &a, &b); //判断是否存在 x, 使得 a^x mod p = b
b %= p;
for(int i = 0; i <= (p << 1); ++i) { //暴力枚举所有的 x
if(ans == b) {
printf("%d", i); //判断是否满足
return 0;
}
ans = (ans * a) % p;
}
printf("No Solution"); //如果枚举完了都没有发现符合条件的 x，那就是无解
return 0;
}
//by CaO
BSGS
尽管 \(\mathcal{O}(p)\) 的枚举可以解决这个问题，但它的性能显然不够好！
如果把 \(p\) 开到 \(10^9\) 级别，那么上面的这个代码就会 TLE 了。所以，我们就要想一个更加优雅的暴力来枚举所有的 \(x\)。
『什么啊，所以还是要暴力枚举吗？』
确实如此。但是，如果我们改变了枚举的方法，那么它的速度就会有质的飞跃。
顺带提一下，在 OI 中，被称为『优雅的暴力』的算法是什么？
『唔……分块吗？』
对，所以，作为另外一种『优雅的暴力』，BSGS 也是基于分块的思想建立起来的。
BSGS 的思路
我们不妨把待求的 \(x\) 分解一下，给定一个 \(A\)，那么 \(x\) 可以写成 \(x=mA-n\) 的形式。
原来的式子，就可以变形为：
两边同时乘以 \(a^n\)，就有：
当然，这里的 \(n\) 应该不大于 \(A\)，\(m\le\left\lceil\dfrac{p}{A}\right\rceil\)。
我们可以先暴力枚举每一个 \(n\)，并把 \(ba^n\bmod p\) 的值用 Hash 表存储起来；然后，再枚举每一个 \(m\)，判断 \(a^{mA}\bmod p\) 在 Hash 表中有没有对应的元素。
那么，枚举 \(n\) 的复杂度为 \(\Theta(A)\)，枚举 \(m\) 的复杂度为 \(\mathcal{O}\left(\dfrac{p}{A}\right)\)。
故总复杂度为 \(\mathcal{O}\left(A+\dfrac{n}{A}\right)\)。由均值不等式，我们取 \(A=\sqrt{n}\)，这个时候时间复杂度取到最优 \(\mathcal{O}(\sqrt{n})\) 了！
看，我们在枚举右面的 \(ba^n\) 的时候，是细细地按照 \(a\) 的指数，一小步一小步枚举的；而枚举左面的 \(a^{mA}\) 的时候，则是以 \(a^A\) 的指数，大步大步走着枚举的。
所以，这个算法就被形象地称为大步小步算法（Baby Step, Giant Step，简称 BSGS）。
当然，也有人戏称为拔山盖世算法，或者北上广深算法 QwQ，但那都是后话了。
代码
这个代码并不好写，有的时候，还要加一些奇奇怪怪的特判，比如[SDOI2011]计算器，这道题如果不加特判的话，只能拿到 \(70\sim 90\text{ pts}\)。这就非常考察选手对算法的理解程度了。
#include <algorithm>
#include <iostream>
#include <cmath>
#include <map>
using namespace std;
typedef long long ll;
map<ll, ll> mp;
ll a, b, p, ans;
ll qpow(ll a, ll b, ll p) { //快速幂
ll ans = 1;
while(b) {
if(b & 1) ans = (ans * a) % p;
a = (a * a) % p; b >>= 1;
}
return ans;
}
ll bsgs(ll a, ll b, ll p) { //BSGS 主体
if(a % p == b % p) return 1; //这个时候 a 和 b 模 p 同余，那么 x = 1 应该是最小的整数解
if(a % p == 0 && b) return -1; //如果 a = 0 但是 b 不为 0，那么很显然无解
ll unit = (ll)ceil(sqrt(p)), tmp = qpow(a, unit, p);
for(int i = 0; i <= unit; ++i)
mp[b] = i, b = (b * a) % p; //Hash 记录每个小步的情况
b = 1;
for(int i = 1; i <= unit; ++i) { //枚举每个大步
b = (b * tmp) % p;
if(mp[b]) return i * unit - mp[b];
}
return -1; //如果没有找到满足题意的，返回无解
}
int main(void) {
scanf("%lld%lld%lld", &p, &a, &b);
ll ans = BSGS(a, b, p);
if(ans == -1) printf("No Solution");
else printf("%lld", ans);
return 0;
}
//by CaO
当然，如果用手写的 Hash，当然是 \(\mathcal{O}(\sqrt{n})\) 的。但，好多时候，手写 Hash 非常麻烦，而且容易出错，所以用 C++ 中的 std::map，因为 std::map 本质上是一棵红黑树，时间复杂度会比 \(\mathcal{O}(1)\) 的 Hash 要多个 \(\log\)，所以用 std::map 实现的 BSGS，时间复杂度实际上是 \(\mathcal{O}(\sqrt{n}\log n)\)。
BSGS 的拓展
然而，上面的 BSGS，是基于 \(a\perp p\) 时的情况给出的，因为 \(a\) 与 \(p\) 不互质的时候，Euler 定理不一定成立。所以这个时候，传统的 BSGS 就会失灵。
我们就需要想一下，如何把 BSGS 扩展到可以解决 \(a\) 与 \(p\) 不互质的情况下呢？
我们把题目所给的这个式子搬出来：
我们不妨设 \(\gcd(a,p)=g\)，那么一定有：
当然，\(b\) 不一定被 \(g\) 整除，这个时候同余方程就是无解的。
但是即使 \(b\) 能被 \(g\) 整除，我们也不能保证此时一定有 \(a\perp \dfrac{p}{g}\)，怎么办？
那我们再对这个方程除以 \(\gcd\left(a,\dfrac{p}{g}\right)\)，不断地重复这个操作，直到 \(a\perp \dfrac{p}{\prod g}\) 即可，即：
这个时候，两边同时乘以 \(\left(\dfrac{a}{\prod_k g_i}\right)^{-1}\)，得到：
这样一个方程，我们就可以用传统的 BSGS 解决了！
像这样，先将底数和模数通过不断除以最大公约数的方式使它们互质，然后用 BSGS 解决离散对数的算法，叫做扩展 BSGS，又称 exBSGS。
可以证明，这个算法的时间复杂度为 \(\mathcal{O}(\sqrt{p}+\log^2 p)=\mathcal{O}(\sqrt{p})\)。
代码
代码如下：
#include <algorithm>
#include <iostream>
#include <cmath>
#include <map>
using namespace std;
typedef long long ll;
map<ll, ll> mp;
ll a, b, p, ans;
ll exgcd(ll a, ll b, ll& x, ll& y) { //求逆元和 gcd 都可以使用扩欧，非常方便 QwQ
if(!b) {
x = 1, y = 0;
return a;
}
ll ans = exgcd(b, a % b, y, x);
y -= a / b * x;
return ans;
}
ll inv(ll a, ll p) { //求 a 模 p 的逆元
ll x, y;
exgcd(a, p, x, y);
return (x % p + p) % p;
}
ll qpow(ll a, ll b, ll p) { //快速幂
ll ans = 1;
while(b) {
if(b & 1) ans = (ans * a) % p;
a = (a * a) % p; b >>= 1;
}
return ans;
}
ll BSGS(ll a, ll b, ll p) { //BSGS 主体，不解释了
ll unit = (ll)ceil(sqrt(p)), tmp = qpow(a, unit, p);
for(int i = 0; i <= unit; ++i)
mp[b] = i, b = (b * a) % p;
b = 1;
for(int i = 1; i <= unit; ++i) {
b = (b * tmp) % p;
if(mp[b]) return i * unit - mp[b];
}
return -1;
}
ll exBSGS(ll a, ll b, ll p) {
ll x, y, g = exgcd(a, p, x, y), k = 0, tmp = 1;
while(g != 1) { //当 gcd(a, p) 不为 1 时，就不断除以 gcd(a, p) 直到 a 与 p 互质
if(b % g) return -1; //b 不能被 gcd(a, p) 整除，当然不互质
++k, b /= g, p /= g, tmp = tmp * (a / g) % p;
if(tmp == b) return k;
g = exgcd(a, p, x, y);
}
//用传统 BSGS 来解决问题
ll ans = BSGS(a, b * inv(tmp, p) % p, p);
if(ans == -1) return -1;
return ans + k;
}
int main(void) {
scanf("%lld%lld%lld", &p, &a, &b);
ll ans = exBSGS(a, b, p);
if(ans == -1) printf("no solution");
else printf("%lld", ans);
return 0;
}
//by CaO
例题
本题目列表会持续更新。
