BSGS学习笔记

VI.BSGS（大步小步算法）

欢迎来到 北上广深 拔山盖世 比赛搞事 不算个事 算法学习现场。

BSGS，全名 Baby Step Giant Step 算法，是用于求解 $a^x\equiv b\pmod p$，其中 $\gcd(a,p)=0$ 的算法。

我们记 $K=\sqrt{p}$ 。则 $x$ 即可表示成 $cK-d$，其中 $0\leq d<K$ 的形式。

于是 $a^{cK-d}\equiv b\pmod p$，则 $(a^K)^c\equiv ba^d\pmod p$。于是我们只需预处理出 $b$ 乘上 $a$ 的所有 $d$ 次幂，扔进哈希表中，然后依次枚举 $a$ 的 $K$ 次幂的 $c$ 次幂进行求解即可。

若使用 map，复杂度为 $O(\sqrt n\log n)$；若使用手写哈希表或 unordered_map，复杂度为 $O(\sqrt n)$。

VI.I.[TJOI2007] 可爱的质数/【模板】BSGS

代码：

#include<bits/stdc++.h>
using namespace std;
int a,b,p,K;
unordered_map<int,int>mp;
int main(){
	scanf("%d%d%d",&p,&a,&b),K=sqrt(p);
	int P=1;
	for(int j=0;j<K;j++)mp[1ll*P*b%p]=j,P=1ll*P*a%p;
	for(int i=1,j=P;i<=K;i++,j=1ll*j*P%p)if(mp.find(j)!=mp.end()){printf("%lld\n",1ll*i*K-mp[j]);return 0;}
	puts("no solution");
	return 0;
}

VI.II.[SDOI2013] 随机数生成器

我们来考虑一下：

日期	原式子	化简后式子
$1$	$X$	$X$
$2$	$aX+b$	$aX+b$
$3$	$a(aX+b)+b$	$a^2X+(a+1)b$
$4$	$a\Big(a(ax+b)+b\Big)+b$	$a^3X+(a^2+a+1)b$
$i$	$\dots$	$a^{i-1}X+b\sum\limits_{j=0}^{i-2}a^j$

于是我们要寻找 $a^{i-1}X+b\sum\limits_{j=0}^{i-2}a^j\equiv T\pmod p$ 的解。

考虑上等比数列求和公式，得到 $a^{i-1}X+b\times\dfrac{a^{i-1}-1}{a-1}\equiv T\pmod p$

稍微整理一下 $a^{i-1}\Big(\dfrac{b}{a-1}+X\Big)\equiv T+\dfrac{b}{a-1}\pmod p$

扔过去，得到 $a^{i-1}\equiv\dfrac{T+\dfrac{b}{a-1}}{\dfrac{b}{a-1}+X}\pmod p$

OK！是BSGS的形式了！直接上就行了……个锤子。

还要特判一堆东西。

1. $X=T$。

这个东西虽然不很特殊，但是因为我们接下来众多特判都要判这个，所以就在开头统一判掉。答案：$1$。

2. $a=0$。

这时，除了第一天是 $X$，其余天都是 $b$。因为我们已经特判了 $X=T$ 的情形，所以只要特判 $b=T$ 的情形，是则答案为 $2$，否则答案为 $-1$。

3. $a=1$。

这时，无法使用等比数列求和，式子为 $X+(i-1)b\equiv T\pmod p$。整理得 $i=\dfrac{T+b-X}{b}$。

但是，还要继续特判！如果 $b=0$，显然只有 $X=T$ 时答案为 $1$（这种情况已在开头特判过），否则答案为 $-1$。

否则，若 $T+b-X=0$，上式出来的结果是 $0$，要手动将其判作 $p$，因为 $p$ 是质数，要 $p$ 个 $b$ 才能绕一圈。

4. $\dfrac{b}{a-1}+X=0$。

此种情形下，无法除过去。除非 $T+\dfrac{b}{a-1}$ 也为 $0$，否则答案即为 $-1$。但是这就意味着有 $X=T$，已经特判过了，所以不用再考虑这种情形，碰到直接判 $-1$ 即可。

5. $i-1=0$。

常规BSGS一般不会考虑这种情形。但是，这就意味着 $i=1$，我们已经在最开头特判过了，故不用考虑。

代码：

#include<bits/stdc++.h>
using namespace std;
int T_T,p,a,b,X,T;
int ksm(int x,int y=p-2){
	int z=1;
	for(;y;y>>=1,x=1ll*x*x%p)if(y&1)z=1ll*z*x%p;
	return z;
}
unordered_map<int,int>mp;
int main(){
	scanf("%d",&T_T);
	while(T_T--){
		scanf("%d%d%d%d%d",&p,&a,&b,&X,&T);
		if(X==T){puts("1");continue;}
		if(!a){
			if(T==b)puts("2");else puts("-1");
			continue;
		}
		if(a==1){
			T=(0ll+T-X+b+p)%p;
			if(!b)puts("-1");else{
				T=1ll*T*ksm(b)%p;
				if(T)printf("%d\n",T);
				else printf("%d\n",p);
			}
			continue;
		}
		b=1ll*b*ksm(a-1)%p;
		X=(X+b)%p,T=(T+b)%p;
		if(!X){puts("-1");continue;}
		T=1ll*T*ksm(X)%p;
//		printf("%d %d\n",a,T);
		mp.clear();
		int K=sqrt(p);
		int P=1;
		for(int j=0;j<K;j++)mp[1ll*P*T%p]=j,P=1ll*P*a%p;
		bool fd=false;
		for(int i=1,j=P;i<=K+1;i++,j=1ll*j*P%p)if(mp.find(j)!=mp.end()){printf("%d\n",i*K-mp[j]+1),fd=true;break;}
		if(!fd)puts("-1");
	}
	return 0;
}

VI.III.CF1106F Lunar New Year and a Recursive Sequence

首先，我们不妨设 $f_k=x$。则，因为递推式里面全都是一坨坨东西的幂再乘一起，初始的元素中又仅有 $f_k$ 可能不为 $1$，所以明显每个元素都可以被表示成 $x^a$ 的形式，其中 $a$ 是我们接下来要求的东西。

我们发现，因为底数是相乘的，所以放在指数上就是相加，可以使用矩阵快速幂处理，只不过因为在指数上，所以要模 $\varphi(p)$。

现在我们知道 $a$ 了，则要求解的是 $x^a\equiv b\pmod p$ 的形式。如何求解？

因为 $998244353$ 有着地球人都知道的原根 $g=3$，所以我们不妨令 $x=g^c$。于是原式转换为 $(g^c)^a\equiv b\pmod p$，继续转换得 $(g^a)^c\equiv b\pmod p$。因为 $g^a$ 已知，所以就是BSGS模板，随便套套就套出 $c$ 来了。

虽然原题上说的是求”任意的 $x$“，但是现行洛谷翻译上给的却是求”最小的 $x$“，我也就照着翻译写了。因此我们这里介绍一下求最小 $x$ 的做法。

首先，$x=g^c$ 一定是一个合法解；又因为 $g$ 最短的循环节是 $\varphi(p)$，所以就有 $\forall t\in\Z,x^a\equiv g^{ac+t\varphi(p)}\equiv b\pmod p$。明显所有这样的 $t$ 构成全部合法解的集合（但是会有重复的）。

这就得到了 $\forall a|t\varphi(p),x\equiv g^{c+\tfrac{t\varphi(p)}{a}}\equiv b\pmod p$。所有这样的 $t$ 得到唯一的解集合。

$a|t\varphi(p)$ 唯一等价于 $\dfrac{a}{\gcd(a,\varphi(p))}|t$。于是我们设 $t=\dfrac{a}{\gcd(a,\varphi(p))}\times i$，就得到了解集为 $\forall i\in\Big[0,\gcd\big(a,\varphi(p)\big)\Big)\cup\Z,g^{c+\tfrac{\varphi(p)}{\gcd(a,\varphi(p))}\times i}$。

现在，我们考虑找到上述集合中最小的数。显然，当 $\gcd\Big(a,\varphi(p)\Big)$ 不算很大（比如 $10^6$ 以内），可以枚举所有的 $i$ 取 $\min$；而当 $\gcd\Big(a,\varphi(p)\Big)$ 较大时，集合中的数又不会很稀疏，可以从 $1$ 开始一个一个枚举check它是否在集合内（当然，使用BSGS）。

时间复杂度比较玄学，取决于你取的分界点。反正我取了 $2^{20}$ 作为分界点。

代码：

#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
const int phi=998244352;
int n,m,A;
struct Matrix{
	int a[110][110];
	Matrix(){memset(a,0,sizeof(a));}
	int*operator[](const int&x){return a[x];}
	friend Matrix operator*(Matrix&u,Matrix&v){
		Matrix w;
		for(int i=1;i<=A;i++)for(int j=1;j<=A;j++)for(int k=1;k<=A;k++)(w[i][j]+=1ll*u[i][k]*v[k][j]%phi)%=phi;
		return w;
	}
	void print()const{for(int i=1;i<=A;i++){for(int j=1;j<=A;j++)printf("%d ",a[i][j]);puts("");}}
}M;
int ksm(int x,int y){
	int z=1;
	for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;
	return z;
}
int KSM(){
	Matrix I;
	for(int i=1;i<=A;i++)I[i][i]=1;
	for(int y=n-A;y;y>>=1,M=M*M)if(y&1)I=I*M;
	return I[A][A];
}
unordered_map<int,int>mp;
int BSGS(int a,int b){
	if(b==0)return 1;
	mp.clear();
	int K=sqrt(mod);
	int P=1;
	for(int i=0;i<K;i++,P=1ll*P*a%mod)mp[1ll*P*b%mod]=i;
	for(int i=1,j=P;i<=K+2;i++,j=1ll*j*P%mod)if(mp.find(j)!=mp.end())return i*K-mp[j];
	return -1;
}
int main(){
	scanf("%d",&A);
	for(int i=1;i<=A;i++)scanf("%d",&M[A-i+1][A]);
	for(int i=2;i<=A;i++)M[i][i-1]=1;
//	M.print();
	scanf("%d%d",&n,&m);
	int a=KSM();
	int c=BSGS(ksm(3,a),m);
	if(c==-1){puts("-1");return 0;}
	int K=__gcd(a,phi),R=phi/K;
	if(K<=1048576){
		int res=mod;
		for(int i=0,j=ksm(3,c),k=ksm(3,R);i<K;i++,j=1ll*j*k%mod)res=min(res,j);
		printf("%d\n",res);
	}else for(int i=1;;i++)if(BSGS(3,i)%R==c%R){printf("%d\n",i);return 0;}
	return 0;
}