BSGS及其扩展
之前写了一篇关于BSGS的学习笔记。因为太过老旧，就想修改一些错误，顺便添上扩展BSGS的部分。可惜博客园不能对已发布的随笔修改编辑器，索性重新发出来。旧文已删。
定义
Baby-Step-Giant-Step算法，简称BSGS算法，又称大步小步算法，用于求方程\(a^x\equiv b(\text{mod }c)\)的最小非负整数解，此过程又称离散对数。
普通BSGS算法只能解决\((a,c)=1\)的情况，对其进行扩展之后则没有这个限制。
原理
当\((a,c)=1\)时，若\(c|a\)且\(b\not=0\)，则显然无解。若\(c\not|\ \ a\)，有\(a^{c-1}\equiv 1(\text{mod }c)\)，即\(c-1\)是一个循环节。因此如果方程有解，必定在\([0,c-2]\)中。
朴素算法
由此得到一个朴素的算法。枚举\([0,c-2]\)中每一个数，如果相等就输出答案。时间复杂度\(O(c)\)。考虑对其进行优化，引入数论分块的思想。
数论分块
设\(x=i\times m+j(0\leq j<m)\)，原方程改写为\(a^{im+j}\equiv b(\text{mod }c)\)。令\(D(i)=a^{im}\)，原方程进一步改写为\(D(i)\cdot a^j\equiv (\text{mod }c)\)。
暴力枚举每个可能的\(i\)，用扩展欧几里得可以求出\(a^j\)。可是我们依然无法求得\(j\)。此时注意到\(j\)的所有可能取值只有\(m\)个，因此我们预先将所有的\(a^j\)存入一个哈希表，当我们确定\(a^j\)时，直接查表可以\(O(1)\)得到是否有对应的\(j\)，以及对应的\(j\)是多少。暴力枚举复杂度\(O(ilog_2c)\)，预处理复杂度\(O(m)\)，总复杂度\(O(ilog_2c+m)\)。注意到\(x\approx i\times m\)，因此当\(m\)取\(\sqrt c\)时（事实上由于\(\log_2c\)的存在不是很严谨）总复杂度最小，为\(O(\sqrt c\ \log_2c)\)。
至此，我们采用以空间换时间，最后通过分块的思想将该算法优化至根号级别。
例题
Luogu2485 [SDOI2011]计算器
题解
操作1：快速幂；操作2：扩展欧几里得；操作3：普通BSGS。
代码
#include<bits/stdc++.h>
#include<tr1/unordered_map>
#define LL long long
using namespace std;
using namespace std::tr1;
unordered_map<int,int>mp;
int t,k;
inline LL qpow(LL a,LL k,LL p){
LL ans=1;
for(;k;a=a*a%p,k>>=1){if(k&1){ans=ans*a%p;}}
return ans;
}
inline LL inv(LL a,LL p){return qpow(a,p-2,p);}
inline LL BSGS(LL a,LL b,LL p){
if(p==1&&b==0)return 0;
if(b==1)return 0;
if(a==0||b==0){
if(a==0&&b==0)return 1;
return -1;
}
LL i,y=(LL)sqrt(p),tmp,ans;
mp.clear();
for(i=0,tmp=1;i<y;i++,tmp=tmp*a%p){
mp[b*tmp%p]=i;
}
for(i=1,ans=tmp;y*(i-1)+1<=p-2;i++,ans=ans*tmp%p){
if(mp.find(ans)==mp.end())continue;
if(i*y-mp[ans]<0)continue;
return i*y-mp[ans];
}
return -1;
}
int main()
{
int i,j;
LL y,z,p,ans;
cin>>t>>k;
while(t--){
scanf("%lld%lld%lld",&y,&z,&p);
y%=p;
if(k==1){printf("%lld\n",qpow(y,z,p));}
else if(k==2){
z%=p;
if(y==0&&z){cout<<"Orz, I cannot find x!\n";continue;}
printf("%lld\n",z*inv(y,p)%p);
}
else{
z%=p;
ans=BSGS(y,z,p);
if(ans==-1){cout<<"Orz, I cannot find x!\n";continue;}
printf("%lld\n",ans);
}
}
return 0;
}
扩展
此时\(c\)不一定是质数，无法保证解的循环性。因此试图将该问题转化为普通BSGS问题。
首先给出一个判定解是否存在的方法：
当\((a,c)\not|\ \ b \quad \text{and}\quad b\not=1\)时，方程无自然数解。
因此，当\((a,c)|b\)时，不妨设\(G=(a,c)\)。
令\(x'=x-1,b'=\frac bG(\frac aG)^{-1},c'=\frac cG\)，则有
递归进行以上过程，直至\((a,c)=1\)，使用普通BSGS即可。
例题
Luogu4195 【模板】exBSGS/Spoj3105 Mod
代码
#include<bits/stdc++.h>
#include<tr1/unordered_map>
#define LL long long
using namespace std;
using namespace std::tr1;
unordered_map<int,int>mp;
#define _BUFFERSIZE 65536
static size_t _pos = _BUFFERSIZE, _len; static char _buf[_BUFFERSIZE];
inline void _getc(char &c){
if(_pos==_BUFFERSIZE){
_pos=0;
_len=fread(_buf,1,_BUFFERSIZE,stdin);
}
c=_pos<_len?_buf[_pos++]:0;
}
template<typename T>
void read(T& x){
x=0;
char c;
do _getc(c); while(!isdigit(c));
do{
x=x*10+(c-'0');
_getc(c);
}while(isdigit(c));
}
void exgcd(LL a,LL b,LL &x,LL &y){
if(!b){x=1;y=0;return;}
exgcd(b,a%b,x,y);
LL tmp=x;x=y;y=tmp-a/b*y;
}
LL inv(LL a,LL p){
LL x,y;
exgcd(a,p,x,y);
return (x%p+p)%p;
}
LL gcd(LL a,LL b){return b?gcd(b,a%b):a;}
LL BSGS(LL a,LL b,LL p){
if(p==1&&b==0)return 0;
if(b==1)return 0;
if(a==0||b==0){
if(a==0&&b==0)return 1;
return -1;
}
LL i,y=(LL)sqrt(p),tmp,ans;
mp.clear();
for(i=0,tmp=1;i<y;i++,tmp=tmp*a%p){
mp[b*tmp%p]=i;
}
for(i=1,ans=tmp;y*(i-1)+1<=p-2;i++,ans=ans*tmp%p){
if(mp.find(ans)==mp.end())continue;
if(i*y-mp[ans]<0)continue;
return i*y-mp[ans];
}
return -1;
}
LL exBSGS(LL a,LL b,LL p){
a%=p;b%=p;
if(b==1)return 0;
if(a==0){if(b==0){return 1;}return -1;}
if(b==0){
LL ans=0,d;
while((d=gcd(a,p))!=1){
ans++;p/=d;
if(p==1)return ans;
}
return -1;
}
LL ans,d,cnt=0;
while((d=gcd(a,p))!=1){
if(b%d)return -1;
cnt++;
p/=d;b=b/d*inv(a/d,p)%p;
}
ans=BSGS(a,b,p);
if(ans==-1)return -1;
return cnt+ans;
}
int main()
{
int i,j;
LL a,p,b,ans;
for(;;){
read(a);read(p);read(b);
if(!(a||b||p))break;
ans=exBSGS(a,b,p);
if(ans==-1){printf("No Solution\n");}
else{printf("%lld\n",ans);}
}
return 0;
}
