BSGS算法(北上广深算法
用于求解 ax≡b(modp)a^x\equiv b~(mod~p)ax≡b (mod p) 的非负整数解 xxx，其中 (a,p)=1(a,p)=1(a,p)=1。
算法过程
令 x=A⌈p⌉−Bx=A\lceil\sqrt{p}\rceil-Bx=A⌈p⌉−B，其中 A,B≤⌈p⌉A,B\leq\lceil\sqrt{p}\rceilA,B≤⌈p⌉。
那么问题也就转化为找到一组 A,BA,BA,B，满足 aA⌈p⌉−B≡b(modp)a^{A\lceil\sqrt{p}\rceil-B}\equiv b~(mod~p)aA⌈p⌉−B≡b (mod p)。
由于 (a,p)=1(a,p)=1(a,p)=1。
所以 aA⌈p⌉≡baB(modp)a^{A\lceil\sqrt{p}\rceil}\equiv ba^B~(mod~p)aA⌈p⌉≡baB (mod p)
暴力枚举 BBB，计算 baBba^BbaB 的所有值，用哈希表记下该值所对应的最大 BBB。然后暴力枚举 AAA，看 aA⌈p⌉a^{A\lceil\sqrt{p}\rceil}aA⌈p⌉ 能否对应到一个 baBba^BbaB。那么一个 xxx 的值就为 A⌈p⌉−BA\lceil\sqrt{p}\rceil-BA⌈p⌉−B。
大概是一个分块的思想，每 p\sqrt{p}p 个数为一块进行考虑。
这样做可行是因为 A⌈p⌉−BA\lceil\sqrt{p}\rceil-BA⌈p⌉−B 覆盖了 aaa 所有的幂次。
第一个出现的 A⌈p⌉−BA\lceil\sqrt{p}\rceil-BA⌈p⌉−B 即是最小非负整数解。因为随着 AAA 增大，答案是递增的。
复杂度是 O(p)O(\sqrt{p})O(p)。
exBSGS算法
当 a,pa,pa,p 不互质时，怎么求解 xxx？
当 (a,p)=1(a,p)=1(a,p)=1 时，aaa 关于 ppp 的逆元存在，因此可以用 BSGSBSGSBSGS 求解。
所以，我们要想办法让 a,pa,pa,p 互质。
设 d1=gcd(a,p)d_1=\gcd(a,p)d1=gcd(a,p)，如果 d1∤bd_1\nmid bd1∤b，那么方程无解，否则我们将方程同时除掉一个 d1d_1d1，即：
ad1ax−1≡bd1(modpd1)\frac{a}{d_1}a^{x-1}\equiv \frac{b}{d_1}~(mod~\frac{p}{d_1})d1aax−1≡d1b (mod d1p)
如果 a,pd1a,\frac{p}{d_1}a,d1p 仍然不互质，则进行进行类似操作：
设 d2=gcd(a,pd1)d_2=\gcd(a,\frac{p}{d_1})d2=gcd(a,d1p)，如果 d2∤bd_2\nmid bd2∤b，那么方程无解，否则我们将方程同时除掉一个 d2d_2d2，即：
a2d1d2ax−2≡bd1d2(modpd1d2)\frac{a^2}{d_1d_2}a^{x-2}\equiv \frac{b}{d_1d_2}~(mod~\frac{p}{d_1d_2})d1d2a2ax−2≡d1d2b (mod d1d2p)
假设这样进行 kkk 次操作后，a,pa,pa,p 终于互质！
记 D=∏i=1kdiD=\prod\limits_{i=1}^{k}d_iD=i=1∏kdi，方程的形式为：
akDax−k≡bD(modpD)\frac{a^k}{D}a^{x-k}\equiv \frac{b}{D}~(mod~\frac{p}{D})Dakax−k≡Db (mod Dp)
这时候，(a,pD)=1(a,\frac{p}{D})=1(a,Dp)=1，于是 (akD,pD)=1(\frac{a^k}{D},\frac{p}{D})=1(Dak,Dp)=1，因此 akD\frac{a^k}{D}Dak 关于 pD\frac{p}{D}Dp 的逆元存在。
于是我们将它丢到方程右边，然后就可以 BSGSBSGSBSGS 求解了！！
这样 BSGSBSGSBSGS 找到的一个最小正整数解，然后加上 kkk，就是可能的答案了。
不过还可能存在 x≤kx\leq kx≤k，使得 ax≡b(modp)a^x\equiv b~(mod~p)ax≡b (mod p) 的情况，这种情况我们 O(k)O(k)O(k) 枚举 xxx 判断一下就可以了。
总的复杂度还是 O(p)O(\sqrt{p})O(p)
这是两道模板题
BSGS模板题
exBSGS模板题
代码如下
#include <bits/stdc++.h>
#include<ext/pb_ds/hash_policy.hpp>
#include<ext/pb_ds/assoc_container.hpp>
using namespace __gnu_pbds;
using namespace std;
typedef long long LL;
typedef unsigned long long uLL;
struct custom_hash {static uint64_t splitmix64(uint64_t x) {x += 0x9e3779b97f4a7c15;x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9;x = (x ^ (x >> 27)) * 0x94d049bb133111eb;return x ^ (x >> 31);}size_t operator()(uint64_t x) const {static const uint64_t FIXED_RANDOM = chrono::steady_clock::now().time_since_epoch().count();return splitmix64(x + FIXED_RANDOM);}
};
LL z = 1;
int read(){int x, f = 1;char ch;while(ch = getchar(), ch < '0' || ch > '9') if(ch == '-') f = -1;x = ch - '0';while(ch = getchar(), ch >= '0' && ch <= '9') x = x * 10 + ch - 48;return x * f;
}
int ksm(int a, int b, int p){int s = 1;while(b){if(b & 1) s = z * s * a % p;a = z * a * a % p;b >>= 1;}return s;
}
unordered_map<int, int, custom_hash> f;
void exgcd(int a, int b, int &x, int &y){if(!b) x = 1, y = 0;else{exgcd(b, a % b, y, x);y -= a / b * x;}
}
int inv(int a, int b){int x, y;exgcd(a, b, x, y);return (x + b) % b;
}
int BSGS(int a, int b, int p){int i, j, t, R = sqrt(p - 1) + 1;j = b;for(i = 0; i <= R; i++){f[j] = i;j = z * j * a % p;}//求出 b * a^i 的值所对应的最大 i j = 1, t = ksm(a, R, p);for(i = 0; i <= R; i++){if(f.count(j) && i * R > f[j]) return i * R - f[j];//第一次出现即是最小正整数解 j = z * j * t % p;}return -1;//如果无解返回 -1
}
int exBSGS(int a, int b, int p){f.clear();int i, ans = -1, g, k = 0, t = 1;while(1){g = __gcd(a, p);if(b % g) return -1;if(g == 1) break;p /= g; b /= g;//每次都除掉 gcd t = z * t * (a / g) % p;//记录 a^k / D k++;}b = z * b * inv(t, p) % p;ans = BSGS(a, b, p);if(ans != -1) ans += k;//求得一个解然后加 k t = 1;//可能存在小于等于 k 的解 for(i = 0; i <= k; i++){if(i % p == b){if(ans == -1) ans = i;else ans = min(ans, i);}}return ans;
}
int main(){int i, j, a, b, p;while(scanf("%d", &a) != EOF){p = read(); b = read();if(!a && !b && !p) break;j = exBSGS(a, b, p);if(j == -1) printf("No Solution\n");else printf("%d\n", j);}return 0;
}