【算法笔记】Baby Step Giant Step（BSGS）及其扩展

合集 - OI算法笔记(13)

1.【算法笔记】K.M.P.字符串匹配2022-01-182.【算法笔记】Tarjan 算法 · 上2022-01-213.【算法笔记】Manacher 算法2022-01-314.【算法笔记】二维凸包2022-01-305.【算法笔记】扩展 Euclid 算法2022-01-276.【算法笔记】Tarjan 算法 · 下2022-01-247.【算法笔记】Lagrange 插值法2022-04-148.【算法笔记】中国剩余定理（CRT）与扩展中国剩余定理（exCRT）2022-04-12

9.【算法笔记】Baby Step Giant Step（BSGS）及其扩展2022-04-09

10.【算法笔记】Gauss 消元 · 行列式2022-04-0311.【算法笔记】Euler-Fermat 定理2022-03-1912.【算法笔记】快速 Fourier 变换（FFT）2022-03-1313.【算法笔记】乘法逆元2022-02-01

收起

• 本文总计约 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^{mA-n}\equiv b\pmod p.$$

两边同时乘以 $a^n$，就有：

$$a^{mA}\equiv ba^n\pmod p.$$

当然，这里的 $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$ 不互质的情况下呢？

我们把题目所给的这个式子搬出来：

$$a^x\equiv b\pmod p.$$

我们不妨设 $\gcd(a,p)=g$，那么一定有：

$$\dfrac{a}{g}\cdot a^{x-1}\equiv\dfrac{b}{g}\quad\left(\bmod {\dfrac{p}{g}}\right).$$

当然，$b$ 不一定被 $g$ 整除，这个时候同余方程就是无解的。

但是即使 $b$ 能被 $g$ 整除，我们也不能保证此时一定有 $a\perp \dfrac{p}{g}$，怎么办？

那我们再对这个方程除以 $\gcd\left(a,\dfrac{p}{g}\right)$，不断地重复这个操作，直到 $a\perp \dfrac{p}{\prod g}$ 即可，即：

$$\dfrac{a}{\prod_k g_i}\cdot a^{x-k}\equiv \dfrac{b}{\prod_k g_i}\quad \left(\bmod \dfrac{p}{\prod_k g_i}\right).$$

这个时候，两边同时乘以 $\left(\dfrac{a}{\prod_k g_i}\right)^{-1}$，得到：

$$a^{x-k}\equiv \left(\dfrac{a}{\prod_k g_i}\right)^{-1}\cdot\dfrac{b}{\prod_k g_i}\quad \left(\bmod \dfrac{p}{\prod_k g_i}\right).$$

这样一个方程，我们就可以用传统的 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

例题

本题目列表会持续更新。

• [TJOI2007]可爱的质数
• [SDOI2011]计算器
• [SDOI2013]随机数生成器
• [CQOI2018]破解 D-H 协议
• [P4195]【模板】扩展 BSGS/exBSGS
• [P4884]多少个 1？