「BSGS」学习笔记
BSGS
拔山盖世、北上广深……
实际上叫大步小步，用于解决高次同余方程，形如：
\[a^x\equiv b\pmod p
\]
求 \(x\)。
设 \(x=i\times t-j\)，有：
\[a^{i\times t-j}\equiv b\pmod p
\]
\[a^{i\times t}\equiv b\times a^j\pmod p
\]
预处理每个 \(j\)，枚举 \(i\) 处理，\(t\) 取 \(\sqrt n\) 最优，复杂度为 \(O(\sqrt n)\)。
点击查看代码
ll BSGS(ll a,ll b,ll P)
{
if(1%p==b%p) return 0;//特判。
unordered_map<int,int>vis;
b%=P;
ll t=sqrt(P)+1,sum=b;
for(int i=0;i<=t-1;i++)
{
vis[sum]=i;
(sum*=a)%=P;
}
a=qpow(a,t,P);
if(!a) return b==0?1:-1;
sum=1;
for(int i=1;i<=t;i++)
{
(sum*=a)%=P;
if(vis.find(sum)!=vis.end()&&i*t-vis[sum]>=0)
return i*t-vis[sum];
}
return -1;
}
例题
P3846 [TJOI2007] 可爱的质数/【模板】BSGS
板子。
P4884 多少个 1？
\(\sum\limits_{i=0}^{n-1}10^i\equiv b\pmod p\) 等价于 \(10^n\equiv b\times 9+1\pmod P\)。
P4454 [CQOI2018] 破解D-H协议。
原根的性质是为了表示 \(a,p\) 互质的，之后直接跑 BSGS 即可。
P3306 [SDOI2013] 随机数生成器
当 \(a=0\) 的时候直接判掉就好了。
先忽略模数，若 \(a\ne 1\)，推式子有：
\[\begin{aligned}
x_n&=ax_{n-1}+b \\
&=a^2x_{n-2}+(a+1)b\\
&=a^{n-1}x_1+(a^{n-2}+a^{n-1}+…+a_0)b
\end{aligned}
\]
右面是一个等比数列，套求和公式 \(s_{n}=\dfrac{x_1(a^n-1)}{a-1}\)。
\[\begin{aligned}
x_n&=a^{n-1}x_1+\dfrac{b(a^{n-1}-1)}{a-1}\\
&=a^{n-1}(x_1+\dfrac{b}{a-1})-\dfrac{b}{a-1}\\
\end{aligned}
\]
故此有 \(a^{n-1}(x_1+\dfrac{b}{a-1})\equiv x_n+\dfrac{b}{a-1}\pmod p\)
需满足 \(a,p\) 互质，发现 \(p\) 为质数，\(\gcd(a,p)\) 只可能等于 \(1\) 或 \(p\)，若等于 \(1\) 不用管，等于 \(p\) 显然无解了，然后跑 \(BSGS\) 即可。
\(a=1\) 的特殊处理即可。
点击查看代码
#include<bits/stdc++.h>
#define ll long long
#define endl '\n'
#define sort stable_sort
using namespace std;
const int N=2e5+10;
template<typename Tp> inline void read(Tp&x)
{
x=0;register bool z=true;
register char c=getchar();
for(;c<'0'||c>'9';c=getchar()) if(c=='-') z=0;
for(;'0'<=c&&c<='9';c=getchar()) x=(x<<1)+(x<<3)+(c^48);
x=(z?x:~x+1);
}
template<typename Tp> inline void wt(Tp x)
{if(x>9)wt(x/10);putchar((x%10)+'0');}
template<typename Tp> inline void write(Tp x)
{if(x<0)putchar('-'),x=~x+1;wt(x);}
ll T,p,a,b,x,t;
ll gcd(ll a,ll b) {return b?gcd(b,a%b):a;}
ll qpow(ll a,ll b,ll P)
{
ll ans=1;
for(;b;b>>=1)
{
if(b&1) (ans*=a)%=P;
(a*=a)%=P;
}
return ans;
}
ll BSGS(ll a,ll b,ll P)
{
if(1%P==b%P) return 0;
unordered_map<ll,ll>vis;
b%=P;
ll t=sqrtl(P)+1,sum=b;
for(int i=0;i<=t-1;i++)
{
vis[sum]=i;
(sum*=a)%=P;
}
a=qpow(a,t,P);
if(!a) return b==0?1:-1;
sum=1;
for(int i=0;i<=t;i++)
{
if(vis.find(sum)!=vis.end()&&i*t-vis[sum]>=0)
return i*t-vis[sum];
(sum*=a)%=P;
}
return -1;
}
signed main()
{
read(T);
while(T--)
{
read(p),read(a),read(b),read(x),read(t);
if(x==t) {puts("1"); continue;}
if(a==0)
{
if(t==b) puts("2");
else puts("-1");
}
else if(a==1)
{
ll g=gcd(b,p);
if(g==p) puts("-1");
else
{
ll ans=((t-x+p)%p)*qpow(b,p-2,p)%p;
write(ans+1),puts("");
}
}
else
{
ll c=b*qpow(a-1,p-2,p)%p;
ll g=gcd(x+c,p);
if(g==p) puts("-1");
else
{
ll ans=BSGS(a,(t+c)%p*qpow(x+c,p-2,p)%p,p);
write(ans==-1?-1:ans+1),puts("");
}
}
}
}
BZOJ4128 Matrix
矩阵乘法 + BSGS。
点击查看代码
#include<bits/stdc++.h>
#define ll long long
#define endl '\n'
#define sort stable_sort
using namespace std;
const int N=75;
template<typename Tp> inline void read(Tp&x)
{
x=0;register bool z=true;
register char c=getchar();
for(;c<'0'||c>'9';c=getchar()) if(c=='-') z=0;
for(;'0'<=c&&c<='9';c=getchar()) x=(x<<1)+(x<<3)+(c^48);
x=(z?x:~x+1);
}
template<typename Tp> inline void wt(Tp x)
{if(x>9)wt(x/10);putchar((x%10)+'0');}
template<typename Tp> inline void write(Tp x)
{if(x<0)putchar('-'),x=~x+1;wt(x);}
ll n,P,a[N][N],b[N][N];
struct aa
{
ll x[N][N];
aa() {memset(x,0,sizeof(x));}
bool operator < (aa another) const
{
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
{
if(x[i][j]<another.x[i][j]) return 1;
if(x[i][j]>another.x[i][j]) return 0;
}
return 0;
}
};
aa mul(aa a,aa b,ll P)
{
aa c;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
for(int k=1;k<=n;k++)
(c.x[i][j]+=(a.x[i][k]*b.x[k][j])%P)%=P;
return c;
}
aa qpow(aa a,ll b,ll P)
{
aa ans;
for(int i=1;i<=n;i++) ans.x[i][i]=1;
for(;b;b>>=1)
{
if(b&1) ans=mul(ans,a,P);
a=mul(a,a,P);
}
return ans;
}
ll BSGS(aa a,aa b,ll P)
{
map<aa,ll>vis;
ll t=sqrt(P)+1;
for(int i=0;i<=t-1;i++)
{
vis[b]=i;
b=mul(b,a,P);
}
a=qpow(a,t,P);
aa sum;
for(int i=1;i<=n;i++) sum.x[i][i]=1;
for(int i=0;i<=t;i++)
{
if(vis.find(sum)!=vis.end()&&i*t-vis[sum]>=0)
return i*t-vis[sum];
sum=mul(sum,a,P);
}
return -1;
}
signed main()
{
read(n),read(P);
aa a,b;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
read(a.x[i][j]);
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
read(b.x[i][j]);
write(BSGS(a,b,P));
}
exBSGS
当 \(a,p\) 不互质的时候，可以不断除去 \(\gcd(a,p)\)，直至其互质，依据 \(a\equiv b\pmod p\) 等价于 \(\dfrac{a}{d}\equiv \dfrac{b}{d}\pmod {\dfrac{p}{d}}\)。
设 \(d_i\) 表示直至 \(a,p\) 互质前第 \(i\) 个 \(\gcd(a,p)\)，总共进行 \(t\) 次，最后有：
\[a^{x-t}\times \dfrac{a^t}{\prod\limits_{i=1}^td_i}\equiv \dfrac{b}{\prod\limits_{i=1}^td_i}\pmod {\dfrac{p}{\prod\limits_{i=1}^td_i}}
\]
之后直接跑 BSGS 即可，中间过程中一旦出现 \(b\) 无法被 \(d_i\) 整除则无解。
发现最后 \(p\) 不一定为质数，所以用 exgcd 求逆元。
点击查看代码
ll BSGS(ll a,ll b,ll P)
{
if(1%P==b%P) return 0;
unordered_map<ll,ll>vis;
b%=P;
ll t=sqrt(P)+1,sum=b;
for(int i=0;i<=t-1;i++)
{
vis[sum]=i;
(sum*=a)%=P;
}
a=qpow(a,t,P);
if(!a) return b==0?1:-1;
sum=1;
for(int i=0;i<=t;i++)
{
if(vis.find(sum)!=vis.end()&&i*t-vis[sum]>=0)
return i*t-vis[sum];
(sum*=a)%=P;
}
return -1;
}
ll exgcd(ll a,ll b,ll &x,ll &y)
{
if(b==0) {x=1,y=0; return a;}
ll d=exgcd(b,a%b,y,x);
y-=a/b*x;
return d;
}
ll inv(ll a,ll P)
{
ll x,y;
exgcd(a,P,x,y);
x=(x%P+P)%P;
return x;
}
ll exBSGS(ll a,ll b,ll P)
{
if(b==1||P==1) return 0;
ll d=gcd(a,P),k=0,mul=1;
while(d!=1)
{
if(b%d!=0) return -1;
k++;
b/=d,P/=d;
mul=mul*(a/d)%P;
if(mul==b) return k;
d=gcd(a,P);
}
ll ans=BSGS(a,b*inv(mul,P)%P,P);
if(ans==-1) return -1;
return ans+k;
}
例题
P4195 【模板】扩展 BSGS/exBSGS
板子。
浙公网安备 33010602011771号