BSGS与扩展BSGS

线性同余方程$bx\equiv 1(modp)$，可以用exgcd或者费马小定理解出。但对于高次同余方程$$a^x\equiv b(modp)$$，需要用BSGS或者扩展BSGS解出。

BSGS

问题描述：对于整数a,b,p，其中a,p互质，求一个最小的非负整数x，满足$a^x\equiv b(modp)$。

BSGS是Baby Step，Giant Step的缩写，又叫大小步算法。因为a,p互质，根据欧拉定理有$a^{\phi (p)}\equiv 1(modp)$成立，也就是说当x>=$\phi (p)$时，会出现循环节，在模p的意义下，$a^x$和$a^{x+\phi (p)}$意义是一样的。因此最小非负整数解如果存在，一定在$[0,\phi (p)-1]$之间。

如果p值不大，完全可以暴力枚举求解。p值若比较大，如接近$2^{32}$，要用到BSGS算法。

设$x=i\ast t-j$，其中$t=\lceil \sqrt{p}\rceil$，0$\le$$j$$\le$$t-1$，0$\le$$i$$\le$$t$。t取定值，通过$i$和$j$的搭配，$x$就可以取遍$[0,\phi (p)-1]$之间所有的数。带入：

$$\begin{array}{rl}{a}^{i\ast t-j}& \equiv b(modp)\\ (a^t{)}^i& \equiv b\ast a^j(modp)\end{array}$$
对于所有的0$\le$$j$$\le$$t-1$，把$b\ast a^{j}modp$的结果插入一个hash表。可以用map或者拉链法实现hash表。

从小到大枚举i的可能值，即0$\le$$i$$\le$$t$，计算$(a^t{)}^i\mathrm{m}\mathrm{o}\mathrm{d}p$，在hash表中查找是否存在对应的$j$，找到的第一个$i$和$j$的组合就是答案$i\ast t-j$。时间复杂度$O(\sqrt{p})$。用map实现hash表时间复杂度还要带一个查找时的log。

怎么保证求出的是最小非负整数解？把$b\ast a^{j}modp$的结果插入hash表时，$j$值是从小到大的，如果前后两个的$val$值相同，则保存的是较大的$j$。枚举$i$是也是从小到大。对于表达式$i\ast t-j$，$i$取最小值，$j$取最大值，结果就是最小的。

核心代码：

int BSGS(int a, int b, int p)
{
//	if (a == 0 && b == 1) return 0;  // 0^0定义不存在，这种情况不需要特判  
	if (a == 0 && b == 0) return 1;
	if (a == 0 && b != 0) return -1;
	q.clear();
	int t = ceil(sqrt(p));
	b %= p;
	for (int j = 0; j < t; ++j)
	{
		int val = (long long)b * szpow(a, j, p) % p;  // 两个int相乘可能会溢出  
		q[val] = j;
	}
	a = szpow(a, t, p);
	for (int i = 0; i <= t; ++i)
	{
		int val = szpow(a, i, p);
		if (q.find(val) != q.end() && i * t - q[val] >= 0) return i * t - q[val];
	}
	return -1;
}

注意事项：在用快速幂求解$a^{b}modp$时，也要开long long，否则两个int相乘可能会溢出；map实现hash表，key值可以取很大，也可以取负值，这点很方便。

int szpow(int a, int b, int p)
{
	long long r = 1, base = a;  // 不开long long会溢出 
	while (b)
	{
		if (b & 1) r = r * base % p;
		base = base * base % p;
		b >>= 1;
	}
	return r % p;
}

拉链法实现的hash表：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <map>
#include <algorithm>
using namespace std;
const int PP = 2887777;
int p, a, b;
int szfirst[PP+10], sznext[70000], cnt;
struct node
{
	int val, j;
}szdate[70000];
inline void szadd(int val, int j)
{
	int u = val % PP;
	++cnt;	
	sznext[cnt] = szfirst[u];
	szfirst[u] = cnt;
	szdate[cnt].val = val, szdate[cnt].j = j;
}
int szfind(int val)
{
	int u = val % PP;
	
	int d = szfirst[u];
	while (d)
	{
		if (szdate[d].val == val) return szdate[d].j;
		d = sznext[d];
	}
	return -1;
}
int szpow(int a, int b, int p)  
{
	long long r = 1, base = a;  // 不开long long会溢出  
	while (b)
	{
		if (b & 1) r = r * base % p;
		base = base * base % p;
		b >>= 1;
	}
	return r % p;
}
int BSGS(int a, int b, int p)
{
	if (a == 0 && b == 1) return 0;
	if (a == 0 && b == 0) return 1;
	if (a == 0 && b > 1) return -1;

	int t = ceil(sqrt(p));  // 2^16 
	b %= p;
	for (int j = 0; j < t; ++j)
	{
		int val = (long long)b * szpow(a, j, p) % p;  // 两个int相乘可能会溢出  
		szadd(val, j);
	}
	a = szpow(a, t, p);
	for (int i = 0; i <= t; ++i)
	{
		int val = szpow(a, i, p);
		int j = szfind(val);
		if (j != -1 && i * t - j >= 0) return i * t - j;
	}
	return -1;
}
int main()
{
	scanf("%d%d%d", &p, &a, &b);
	int ans;
	ans = BSGS(a, b, p);
	if (ans == -1) printf("no solution\n");
	else printf("%d\n", ans);
	return 0;
}

题目：
1.洛谷：P3846 【模板】BSGS
2.洛谷：P2485 计算器

扩展BSGS

当a,p不一定互质时，需要用到扩展BSGS。总体思路：将a,p约分，使其互质，然后用BSGS求解。

设$d=gcd(a,p)$，根据同余性质，下式成立：
$$\frac{a^x}{d}\equiv \frac{b}{d}(mod\frac{p}{d})$$
此时如果$d\nmid b$，要么此时$b=1$，$x=0$，否则根据裴蜀定理可知，上式无解。因此，在最开始的时候判断$b$是否为1，以后凡是遇到$d\nmid b$的情况，说明无解。
$$\begin{array}{rl}\frac{a}{d}{a}^{x-1}& \equiv \frac{b}{d}(mod\frac{p}{d})\\ t=\frac{a}{d}b& =\frac{b}{d}p=\frac{p}{d}\\ t{a}^{x-1}& \equiv b(modp)\end{array}$$
继续判断此时a,p是否互质，如果不互质，按照上面的方法继续约分。直到互质。得到同余方程：
$$\frac{a^k}{{\prod }_{i=1}^k{d}_i}{a}^{x-k}\equiv \frac{b}{{\prod }_{i=1}^k{d}_i}(mod\frac{p}{{\prod }_{i=1}^k{d}_i})$$
化减为：
$$t{a}^{x1}\equiv b1(modp1)$$
按照普通的BSGS解即可。最后求出的答案再加上k。此时求出的答案是$\ge k$的，还需要特判$[1,k]$的解的情况。

while (1)
	{
		scanf("%d%d%d", &a, &p, &b);
		if (!a && !p && !b) break;
		if (b == 1)  // 特判b为1的情况 
		{
			printf("0\n"); continue;  
		}
		int ans = 0, k = 0, d, t = 1;
		a %= p;
		b %= p;
		while (1)
		{
			d = gcd(a, p);
			if (d == 1) break; // 如果a、p互质，按照普通BSGS求解  
			++k;
			if (b % d) // 如果d不能整除b，根据裴蜀定理，无解  
			{
				printf("No Solution\n");
				ans = -1; 				
				break;
			}
			t = (long long) t * (a / d) % p; 
		//	t = a / d * t;  // 这样写有可能溢出，只能过40%的数据  
			b /= d; p /= d;
		}
		// 首先特判[1, k]是否满足 
		if (ans != -1) 
		for (int i = 1; i <= k; ++i)
		{
			if (szpow(a, i, p) == b % p) 
			{
				printf("%d\n", i);
				ans = -1;
				break;
			}
		} 
		if (ans == -1) continue;
		ans = BSGS(a, b, p, t);
		if (ans == -1) printf("No Solution\n");
		else printf("%d\n", ans + k);  // 最终结果还要加k  
	}	

注意事项：在求解BSGS的时候，不要忘记乘上t。

题目：
1.洛谷：P4195 【模板】扩展BSGS