exBSGS算法学习笔记

BSGS算法(北上广深算法

用于求解 $a^x\equiv b\left(modp\right)$ax≡b (mod p) 的非负整数解 $x$x，其中 $(a,p)=1$(a,p)=1。

算法过程

令 $x=A\lceil \sqrt{p}\rceil -B$x=A⌈p ​⌉−B，其中 $A,B\le \lceil \sqrt{p}\rceil$A,B≤⌈p ​⌉。
那么问题也就转化为找到一组 $A,B$A,B，满足 $a^{A\lceil \sqrt{p}\rceil -B}\equiv b\left(modp\right)$aA⌈p ​⌉−B≡b (mod p)。
由于 $(a,p)=1$(a,p)=1。
所以 $a^{A\lceil \sqrt{p}\rceil }\equiv b{a}^B\left(modp\right)$aA⌈p ​⌉≡baB (mod p)
暴力枚举 $B$B，计算 $b{a}^B$baB 的所有值，用哈希表记下该值所对应的最大 $B$B。然后暴力枚举 $A$A，看 $a^{A\lceil \sqrt{p}\rceil }$aA⌈p ​⌉ 能否对应到一个 $b{a}^B$baB。那么一个 $x$x 的值就为 $A\lceil \sqrt{p}\rceil -B$A⌈p ​⌉−B。
大概是一个分块的思想，每 $\sqrt{p}$p ​ 个数为一块进行考虑。
这样做可行是因为 $A\lceil \sqrt{p}\rceil -B$A⌈p ​⌉−B 覆盖了 $a$a 所有的幂次。
第一个出现的 $A\lceil \sqrt{p}\rceil -B$A⌈p ​⌉−B 即是最小非负整数解。因为随着 $A$A 增大，答案是递增的。
复杂度是 $O\left(\sqrt{p}\right)$O(p ​)。

exBSGS算法

当 $a,p$a,p 不互质时，怎么求解 $x$x？
当 $(a,p)=1$(a,p)=1 时， $a$a 关于 $p$p 的逆元存在，因此可以用 $BSGS$BSGS 求解。
所以，我们要想办法让 $a,p$a,p 互质。
设 $d_1=\gcd (a,p)$d1​=gcd(a,p)，如果 $d_1\nmid b$d1​∤b，那么方程无解，否则我们将方程同时除掉一个 $d_1$d1​，即：
$\frac{a}{d_1}{a}^{x-1}\equiv \frac{b}{d_1}\left(mod\frac{p}{d_1}\right)$d1​a​ax−1≡d1​b​ (mod d1​p​)
如果 $a,\frac{p}{d_1}$a,d1​p​ 仍然不互质，则进行进行类似操作：
设 $d_2=\gcd (a,\frac{p}{d_1})$d2​=gcd(a,d1​p​)，如果 $d_2\nmid b$d2​∤b，那么方程无解，否则我们将方程同时除掉一个 $d_2$d2​，即：
$\frac{a^2}{d_1{d}_2}{a}^{x-2}\equiv \frac{b}{d_1{d}_2}\left(mod\frac{p}{d_1{d}_2}\right)$d1​d2​a2​ax−2≡d1​d2​b​ (mod d1​d2​p​)
假设这样进行 $k$k 次操作后， $a,p$a,p 终于互质！
记 $D={\prod }_{i=1}^k{d}_i$D=i=1∏k​di​，方程的形式为：
$\frac{a^k}{D}{a}^{x-k}\equiv \frac{b}{D}\left(mod\frac{p}{D}\right)$Dak​ax−k≡Db​ (mod Dp​)
这时候， $(a,\frac{p}{D})=1$(a,Dp​)=1，于是 $(\frac{a^k}{D},\frac{p}{D})=1$(Dak​,Dp​)=1，因此 $\frac{a^k}{D}$Dak​ 关于 $\frac{p}{D}$Dp​ 的逆元存在。
于是我们将它丢到方程右边，然后就可以 $BSGS$BSGS 求解了！！
这样 $BSGS$BSGS 找到的一个最小正整数解，然后加上 $k$k，就是可能的答案了。
不过还可能存在 $x\le k$x≤k，使得 $a^x\equiv b\left(modp\right)$ax≡b (mod p) 的情况，这种情况我们 $O\left(k\right)$O(k) 枚举 $x$x 判断一下就可以了。
总的复杂度还是 $O\left(\sqrt{p}\right)$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;
}