拓展BSGS

原创已于 2024-08-17 17:38:55 修改 · 440 阅读

2 ·

CC 4.0 BY-SA版权

文章标签：

#算法 #c++

于 2022-03-03 21:34:11 首次发布

数学专栏收录该内容

4 篇文章

订阅专栏

本文详细介绍了拓展BSGS算法，用于解决不互质情况下ax ≡ b (mod p)的问题。通过转化和分块策略，将问题转换为可应用BSGS的形式，降低枚举复杂度至平方根级别。文章包含算法流程、改进思路和实现中的注意事项，包括特判条件和结果调整。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

拓展BSGS

BSGS解决像这样的问题：
对于a,p,b,求解形如a^x ≡b(mod p)的问题，其中a,p要求互质
对于a,p不互质的情况，我们引入拓展BSGS算法

一.前置性质

a^x 对p的模随着x周期性变化，并且周期不大于p。根据鸽笼原理，由于 a^x %p∈[0,p-1]值域大小为p，故x取1~p+1时必有两个a^x 的余数是一样的，于是这两个x的差值绝对值就是周期。
若a^x ≡b(mod p)有解等价于 a^x ≡b(mod p)在x∈[0,p-1]范围内有解，由上面的性质显然成立。这告诉我们，枚举范围为[0,p-1]即可。

二.BSGS算法流程

1.思想：分块，块大小为B=p^1/2（上取整）,一部分预处理出一个hash表，一部分在hash里面查找，从而将枚举复杂度降为根号级别.
2.思路：

首先，令x=mB-n,其中n∈[0,B)，于是原问题转化成a^mB ≡b*aⁿ 的形式，我们接下来针对这一问题进行求解
接下来在[0,B)的范围内枚举n，通过hash表构造出b*aⁿ %p的值到n的映射，即Baby Step的预处理。
随后执行Giant Step,在[1,B]范围内枚举m,求出a^mB %p的值，在hash表里查找这个值是否有对应的n，如有，则x=mB-n就是一组合法的解。如果始终没有找到x，那么说明无解

3.说明
至于B的取值要是上取整，这是为了保证B² 能覆盖到[0,p-1]内的每一个数字。在hash表O(1)调取插入时，复杂度为根号级别

三.拓展BSGS

1.改进思路：我们对于a,p不互质的情况显然不能应用朴素BSGS算法，因为此时aⁿ 的逆元不一定存在，最开始的问题转化就没法保证正确性。那么我们接下来进行一些转化，使得该问题再次成为熟悉的朴素BSGS算法所能解决的问题
2.具体流程

首先，将原问题转化为：m* a^x + n* p=b。取g=gcd(a,p),令a’=a/g ,p’=p/g,b’=b/g。将a’加入一个栈中，并将p赋值为p’，b赋值为b’。
重复上述过程，直至a,p互质。如果过程中b出现b%g!=0的情况，则表明原方程无解。此时，设A表示栈中元素的乘积对p取模的结果，cnt表示栈中元素个数，则原问题转化为：mA a^x-cnt + n* p = b,此时a,p互质。显然这个问题等价于a^x-cnt≡b A^-1 (mod p),可以使用朴素BSGS算法求解。不过，如果上述过程中出现b==A的情况，直接返回当时的cnt值就好了。

注意：

实际写程序的时候不需要开栈统计a’,只需要两个变量分别维护A和cnt即可。
由于p在每一步操作中都变化，而最后A要求对最终的p取模的余数，所以你可能对这个方式的正确性产生怀疑。不过由于p始终是那个最终的p的倍数，所以都存起来统一取模和拿每一步的p当作模数得到的最终结果是一致的。
显然复杂度相较于朴素BSGS没有本质上的差别。

洛谷模板（费了好大劲才过最后一个hack测试点）
说一些坑点：
1.先把a,b对p取模，然后执行：
if(b== 1 || p==1) return 0;
hack数据：21 48 49，输出应为0
2.只有这一处需要添加特判，在做完预处理之后不需要再进行上述特判。
hack数据：18 14 16，输出应为2
3.执行BSGS后，得到的结果X是真正结果ans加上“栈中元素数量cnt”，所以不能直接输出。特别注意，得到的BSGS结果不能直接加上cnt,而需要先调正再加，因为需要保证X>cnt
hack数据：38 20 48，输出结果应为5

#include<bits/stdc++.h>
using namespace std;
#define il inline
#define re register
typedef long long LL;
il int read(){
	int s=0,w=1;char c=getchar();
	while(c<'0'||c>'9'){ if(c=='-') w=-1;c=getchar();}
	while(c>='0'&&c<='9'){ s=(s<<1)+(s<<3)+c-'0';c=getchar();}
	return s*w;
} 
LL a,b,p;
map<LL,LL>mp;
il LL gcd(LL x,LL y){
	if(x<y) swap(x,y);
	return y?gcd(y,x%y):x;
}
il int exgcd(int &x,int &y,int a,int b){
	if(b==0){ x=1,y=0;return a; }
	int g=exgcd(x,y,b,a%b);
	int z=x;x=y,y=1LL*z-1LL*(a/b)*y;
	return g;
}
il int Ni(int val,int mod){
	int x,y,g;
	g=exgcd(x,y,val,mod);
	return (x%mod+mod)%mod;
}
il LL qpow(LL a,LL x,LL mod){
	LL res=1LL;
	while(x){
		if(x&1) res=(res*a)%mod;
		a=(a*a)%mod,x>>=1;
	} return res%mod;
}
il int exBSGS(){
	//prework
	LL A=1LL,g=0;int cnt=0;
	a%=p,b%=p;
	if(b==1 || p==1) return 0;//必要的特判，此时答案为0，参看下面的hack数据 
	for(re int g=gcd(a,p);g>1;g=gcd(a,p)){
		if(b%g) return -1;
		//printf("A= %lld,b= %lld,p= %lld\n",A,b,p);
		p/=g,b/=g,A=a/g*A%p,++cnt;
		if(A==b) return cnt; 
	}
	b=b*Ni(A,p)%p;
	//printf("a= %lld,A= %lld,b= %lld,p= %lld,cnt= %d\n",a,A,b,p,cnt);
	//BSGS
	mp.clear(),mp[b]=0;
	LL SQ=ceil(sqrt(p)),G=qpow(a,SQ,p);
	for(re int i=1;i<=SQ;++i){
		b=b*a%p,mp[b]=i;
	}
	LL mul=1LL;
	for(re int i=1;i<=SQ;++i){
		mul=mul*G%p;
		if(mp.count(mul)) return (i*SQ-mp[mul]+p)%p+cnt;
	} return -1;
}
int main()
{
	//freopen("P4195_11.in","r",stdin);
	//freopen("LZ.out","w",stdout);
	while(1){
		a=read(),p=read(),b=read();
		if(a==0 && b==0 && p==0) break;
		int ans=exBSGS();
		if(ans==-1) printf("No Solution\n");
		else printf("%d\n",ans);
	}
}
/*
源自洛谷上#11 
21 48 49
18 14 16
38 20 48

0
2
5
*/