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的形式了!直接上就行了……个锤子。
还要特判一堆东西。
- \(X=T\)。
这个东西虽然不很特殊,但是因为我们接下来众多特判都要判这个,所以就在开头统一判掉。答案:\(1\)。
- \(a=0\)。
这时,除了第一天是 \(X\),其余天都是 \(b\)。因为我们已经特判了 \(X=T\) 的情形,所以只要特判 \(b=T\) 的情形,是则答案为 \(2\),否则答案为 \(-1\)。
- \(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\) 才能绕一圈。
- \(\dfrac{b}{a-1}+X=0\)。
此种情形下,无法除过去。除非 \(T+\dfrac{b}{a-1}\) 也为 \(0\),否则答案即为 \(-1\)。但是这就意味着有 \(X=T\),已经特判过了,所以不用再考虑这种情形,碰到直接判 \(-1\) 即可。
- \(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;
}