1. 引入
Baby Step Giant Step算法（简称BSGS），用于求解形如\(a^x\equiv b\pmod p\)（\(a,b,p\in \mathbb{N}\)）的同余方程，即著名的离散对数问题。
本文分为 \((a,p)=1\) 和 \((a,p)\neq 1\) 两种情况讨论。
2. 方程 \(a^x\equiv b \pmod p\) 的解性
因为若 \(a^{x}\equiv a^{x+n}\pmod p\)，则 \(a^{x+i}\equiv a^{x+n+i}\)。由抽屉原理以及同余的性质得，\(a^x\) 对 \(p\) 的模具有周期性，最大周期不超过 \(p\)。由这个周期性可得，方程有解，等价于在 \([0,p)\) 中有解。由此我们只需要试图找到 \([0,p)\) 中的最小自然数解，就可以得到方程是否有解，以及根据指数模的周期性得到方程的通解。
下面我们只讨论求解方程的最小自然数解。
3. 求解 \(a^x\equiv b\pmod p\) \((a,p)=1\)
Baby Step Giant Step它的名字就告诉我们，这个算法要将问题规模由大化小。我们令 \(t=\lceil \sqrt{p}\rceil\)，那么我们所有 \([0,p)\) 的自然数都可以包含在集合 \(\{x|x=it-j,\ i\in [1,t],\ j \in (0,t]\}\) 中。我们只需要验证这个集合中是否存在解，即可验证原方程是否存在解。
我们考虑将方程表示为
\[
a^{it-j}\equiv b\pmod p
\]
根据 \((a,p)=1\)，我们有
\[a^{it}\equiv ba^j\pmod p\]
我们发现等号两边的数我们都可以在 \(O(\sqrt p)\) 的时间内枚举出来。
我们只需要 \(O(\sqrt p)\) 枚举出所有 \(ba^j\bmod p(\ j \in (0,t])\)，然后把它们插入到哈希表当中，哈希表的对应位置储存这个值对应的最大的 \(j\)（因为要求最小自然数解，所以在 \(i\) 相同的时候要使 \(j\) 最大）。
然后 \(O(\sqrt p)\) 从小到大枚举我们需要的 \(a^{it}\bmod p(\ i\in [1,t])\)，并在哈希表中查询是否有相等的值，若存在则取最小的 \(i\) 及对应最大的 \(j\)，并将 \(it-j\) 作为最小自然数解；否则若不存在则说明无自然数解。
inline int solve_BSGS(const int &a, const int &b, const int &p)
{
int t = ceil(sqrt(p));
std::map<int, int> hash;
//map实现hash表
hash.clear();
int tmp = b;
for (int i = 1; i <= t; ++i)
{
tmp = 1LL * tmp * a % p;
hash[tmp] = i;
}
//插入b*a^j
int pw = qpow(a, t, p);
tmp = pw;
for (int i = 1; i <= t; ++i)
{
if (hash.find(tmp) != hash.end())
return i * t hash[tmp];
tmp = 1LL * tmp * pw % p;
}
//查询a^(it)
return -1; //返回无解
}4. 求解 \(a^x\equiv b\pmod p\) \((a,p)\neq 1\)
对于这个方程，我们不能像上面那样求解的原因就是 \(a\) 在模 \(p\) 意义下不存在逆元，不能将 \(a^{it-j}\) 表示为 \(a^{it}\times a^{-j}\)。那么我们从不定方程的角度分析这个同余方程。
这个方程等价于
\[a^x+py=b\]
令 \(a_1=(a,p)\)，令原方程化为
\[a^{x-1}\frac{a}{a_1}+\frac{p}{a_1}y=\frac{b}{a_1}\]
若此时 \((a,\frac{p}{a_1})\neq 1\)，那么我们接着化成
\[a^{x-2}\frac{a^2}{a_1a_2}+\frac{p}{a_1a_2}y=\frac{b}{a_1a_2}\]
同理不断进行这样的操作，最后我们达到 \((a,\frac{p}{a_1a_2\dots a_n})=1\)的目标，并将方程化为
\[a^{x-n}\frac{a^n}{a_1a_2\dots a_n}+\frac{p}{a_1a_2\dots a_n}y=\frac{b}{a_1a_2\dots a_n}\]
然后记 \(a'=\frac{a}{a_1a_2\dots a_n}\)，\(p'=\frac{p}{a_1a_2\dots a_n}\)，\(b'=\frac{b}{a_1a_2\dots a_n}\)，那么原不定方程可以化为同余方程
\[a^{x-n}a'\equiv b'\pmod{p'}\]
显然 \((a',p')=1\)，因此我们可以写成
\[a^{x-n}\equiv b'(a')^{-1}\pmod{p'}\]
然后就可以用互质的方法解决了。可以发现，每次在 \(p\) 中除去一个最大公约数，每次都会有至少同一个质因子的次数减少 \(1\)，那么在int范围内，\(n\) 最多只会取到 \(30\)。以上两个情况的时间复杂度均为 \(O(\sqrt{p}\log_2{\sqrt p})\)，因为用map实现哈希，可以做到更优秀。
有一些值得注意的地方：
1. 在我们令 \(b\) 除以一个最大公约数 \(d\) 时，若 \(d\nmid b\)，结合求解二元不定方程的知识，我们判定方程无自然数解。
2. 我们在不互质情况的化简操作中，已经假定了 \(x\ge n\)，所以方程才能写成那样的互质形式。对于 \(x\le n\) 的情况，我们应当提前枚举判断。
3. 注意在求解方程之前将 \(a,b\) 对 \(p\) 取模，注意 \(a\) 取模后为 \(0\) 的情况。
4. 对于对时间限制要求较为紧的题目，应当使用更为优秀的哈希表实现方式。
模板题：SP3105 MOD Power Modulo Inverted
洛谷链接：https://www.luogu.org/problemnew/show/SP3105
代码：
//map水过
#include <map>
#include <cmath>
#include <cstdio>
#include <cctype>
#include <cstring>
#include <cstdlib>
#include <iostream>
#include <algorithm>
int a, p, b;
inline int qpow(int b, int p, const int &mod)
{
int res = 1;
for (; p; p >>= 1, b = 1LL * b * b % mod)
if (p & 1)
res = 1LL * res * b % mod;
return res;
}
inline int ex_gcd(const int &a, const int &b, int &x, int &y)
{
if (!b)
return x = 1, y = 0, a;
int res = ex_gcd(b, a % b, y, x);
return y -= a / b * x, res;
}
inline int solve_equ(const int &a, const int &b, const int &c)
{
int x, y;
int d = ex_gcd(a, b, x, y);
if (c % d != 0) return -1;
int mod = b / d;
return (1LL * c / d * x % mod + mod) % mod;
}
inline int solve_BSGS(const int &a, const int &b, const int &p)
{
int t = ceil(sqrt(p));
std::map<int, int> hash;
hash.clear();
int tmp = b;
for (int i = 1; i <= t; ++i)
{
tmp = 1LL * tmp * a % p;
hash[tmp] = i;
}
int pw = qpow(a, t, p);
tmp = pw;
for (int i = 1; i <= t; ++i)
{
if (hash.find(tmp) != hash.end())
return i * t - hash[tmp];
tmp = 1LL * tmp * pw % p;
}
return -1;
}
inline bool check()
{
int k = 1 % p;
for (int i = 0; i <= 40; ++i)
{
if (k == b)
return printf("%d\n", i), true;
k = 1LL * k * a % p;
}
if (!a)
return puts("No Solution"), true;
return false;
}
int main()
{
while (scanf("%d%d%d", &a, &p, &b), a || p || b)
{
a %= p, b %= p;
if (check())
continue;
int d;
int ap = 1, n = 0;
bool flg = false;
while ((d = std::__gcd(a, p)) != 1)
{
++n;
ap = 1LL * ap * (a / d) % p;
p /= d;
if (b % d)
{
flg = true;
break;
}
b /= d;
}
if (flg)
puts("No Solution");
else
{
int res = solve_BSGS(a, 1LL * b * solve_equ(ap, p, 1) % p, p);
if (res == -1)
puts("No Solution");
else
printf("%d\n", res + n);
}
}
return 0;
}