(ex)BSGS/(扩展)大步小步算法 学习笔记
(ex)BSGS/(扩展)大步小步算法 学习笔记
在即将暂时退役之际杀掉了P4195的毒瘤模板题,于是来写篇学习笔记。
谨此为我初中三年摆烂的OI生涯画上一个句号。(距离中考还有20天!)
BSGS
求\(a^x\equiv b\pmod p\)的非负整数解,其中\(a, p\)互质。
算法思路
我们不妨令\(t=\lceil{\sqrt{p}\rceil}\),\(j\lt t\),\(i\leq t\)
原式转化为\(a^{it-j} \equiv b \pmod p\)
即 \(\left(a^t\right)^i\equiv b\cdot a^j \pmod p\)
于是我们可以这么在\(\Theta\left(\sqrt{n}\right)\)(不考虑hash表)内求出解:
-
枚举\(j \in [0, t)\),求出所有的\(b\cdot a^j \mod p\),用hash表记录下来
-
枚举\(i \in [0, t]\),求出所有的\(\left(a^t\right)^i \mod p\),在hash表中查找是否有对应的\(j\)值
需要注意的是,当\(a^t \mod p=0\)时,若\(b=0\),则解为\(x=1\);
否则无解。
正确性讨论
由Euler定理,我们有\(a^{\varphi\left(p\right)}\equiv 1\pmod p\),其中\(\varphi\left(x\right)\)为Euler函数。
也就是说,\(\mod p\)意义下,\(a^x\mod p\)在\(x=n\cdot\varphi\left(p\right)\)(\(n\)遍历非负整数)后循环。
我们知道对于任意整数\(x \gt 1\),\(x>\varphi\left(x\right)\),而我们取的\(it-j\)可以遍历\([0, x]\),因此能够取到\(a^x \mod p\)的一切情况,故一定不会漏解。
代码实现
#include <bits/stdc++.h>
using namespace std;
#define int long long
int qpow(int a, int n, int p) {
a%=p;
int res=1;
while (n) {
if (n&1) res=res*a%p;
a=a*a%p;
n>>=1;
}
return res;
}
signed bsgs(int a, int b, int p) {
b%=p;
int t=sqrt(p)+1;
unordered_map<int, int> hash; hash.clear();
for (int j=0; j<t; ++j) {
int power=b*qpow(a,j,p)%p;
hash[power]=j;
}
a=qpow(a,t,p);
if (a==0) return b==0?1:-1;
for (int i=0; i<=t; ++i) {
int val=qpow(a,i,p);
int j=hash.find(val)==hash.end()?-1:hash[val];
if (j>=0 && i*t-j>=0) return i*t-j;
}
return -1;
}
signed main() {
int a, b, p;
cin>>p>>a>>b;
int res=bsgs(a,b,p);
if (res==-1) puts("no solution");
else cout<<res<<endl;
return 0;
}
这份代码是可以通过P3846的。我们可以考虑对它进行优化。
我们发现枚举\(a^j\)和\(\left(a^t\right)^i\)的时候,\(i, j\)是递增的,于是我们可以直接用一个变量来维护\(a^j\)和\(\left(a^t\right)^i\)。
另外,用unordered_map
来实现hash表似乎会快一些。
inline int bsgs(int a, int b, int p) {
int t=(int)(sqrt(p))+1;
unordered_map<int,int> h; h.clear();
int powa=1;
for (reg int j=0; j<t; ++j) {
int val=powa*b%p;
h[val]=j;
powa=powa*a%p;
}
a=qpow(a,t,p);
powa=1;
for (reg int i=0; i<=t; ++i) {
int val=powa%p;
int j=h.find(val)==h.end()?-1:h[val];
if (j>=0 && i*t-j>=0) return i*t-j;
powa=powa*a%p;
}
return -1;
}
为exBSGS进行修改
我们现在来考虑\(ka^x\equiv b\pmod p\)(\(a, p\)互质,\(k\)为正整数)的式子,我们可以同样地将它们变形为
于是稍微修改一下上述代码就可以了。
inline int bsgs(int a, int b, int k, int p) {
int t=(int)(sqrt(p))+1;
unordered_map<int,int> h; h.clear();
int powa=1;
for (reg int j=0; j<t; ++j) {
int val=powa*b%p;
h[val]=j;
powa=powa*a%p;
}
a=qpow(a,t,p);
powa=1;
for (reg int i=0; i<=t; ++i) {
int val=k*powa%p;
int j=h.find(val)==h.end()?-1:h[val];
if (j>=0 && i*t-j>=0) return i*t-j;
powa=powa*a%p;
}
return -1;
}
exBSGS
求\(a^x\equiv b\pmod p\)的非负整数解,其中\(a, p\)未必互质。
算法思路
考虑利用同余式的约化性质,转换成朴素的BSGS
来求解。
我们有如下同余式的约化性质:
若\(a\equiv b\pmod p\),\(d\mid a,d \mid b\),则
\(\frac{a}{d}\equiv\frac{b}{d}\pmod {\frac{p}{\gcd(d,p)}}\)
我们回到\(a^x\equiv b\pmod p\),令\(d_1=\gcd(a,p)\),
当\(d_1 \mid b\),我们可以变形为
\(\frac{a}{d_1}\cdot a^{x-1}\equiv \frac{b}{d_1} \pmod {\frac{p}{d_1}}\)
(若\(d_1 \nmid b\),立刻得出无解)
若\(a,\frac{p}{d_1}\)仍不互质,我们可以继续令\(d_2=\gcd(a,\frac{p}{d_1})\),\(\cdots\),直到\(a,\frac{p}{d_1d_2\cdots d_k}\)互质为止。
我们设一共进行了\(k\)次约化,\(D=d_1d_2\cdots d_k\)(此时\(a\),\(\frac{p}{D}\)互质),原式可变形为
\(\frac{a^k}{D}\cdot a^{x-k} \equiv \frac{b}{D} \pmod {\frac{p}{D}}\)
我们令\(k=\frac{a^k}{D}, X=x-k, B=\frac{b}{D}, P=\frac{p}{D}\),即
\(ka^X\equiv B \pmod P\)
于是可以利用修改后的BSGS
算法来求解。
注意求解之后得到\(X=x-k\),故\(x=X+k\)。
Trick
\(\frac{a^k}{D}=\frac{a}{d_1}\frac{a}{d_2}\cdots\frac{a}{d_k}\),于是可以在每一个循环内单独计算。
用cout<<'\n'
似乎会比用cout<<endl
快一些。
可以预先将b%=p, a%=p
,若取模过后\(b==1\)或者\(p==1\),那么显然\(x=0\)为解。
我们记\(D_k=\frac{a}{d1}\frac{a}{d2}\cdots \frac{a}{d_k}\),当\(\frac{a^k}{D^k}==\frac{b}{D_k}\)时,有
即
于是\(x=k\)为解。其中\(k\)是正在进行的第\(k\)次约化。
代码实现
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define reg register
inline int qpow(int a, int n, int p) {
a%=p; int res=1;
while (n) {
if (n&1) res=res*a%p;
a=a*a%p; n>>=1;
}
return res;
}
inline int bsgs(int a, int b, int k, int p) {
int t=(int)(sqrt(p))+1;
unordered_map<int,int> h; h.clear();
int powa=1;
for (reg int j=0; j<t; ++j) {
int val=powa*b%p;
h[val]=j;
powa=powa*a%p;
}
a=qpow(a,t,p);
powa=1;
for (reg int i=0; i<=t; ++i) {
int val=k*powa%p;
int j=h.find(val)==h.end()?-1:h[val];
if (j>=0 && i*t-j>=0) return i*t-j;
powa=powa*a%p;
}
return -1;
}
inline int exbsgs(int a, int b, int p) {
a%=p, b%=p;
if (b==1 || p==1) return 0;
int A=a, NA=1, B=b, P=p, k=0, D=1;
while (__gcd(a,P)>1) {
int d=__gcd(a,P);
if (B%d) return -1;
k++; A/=d, B/=d, P/=d, NA=NA*(a/d)%P, D=D*d%P;
if (NA==B) return k; // NA就是上文提到的(a^k)/(D^k)
}
int res=bsgs(a,B,NA,P);
if (res==-1) return res;
if ((qpow(a,res+k,p)-b)%p) return -1;
return res+k;
}
signed main() {
ios::sync_with_stdio(false);
cin.tie(NULL); cout.tie(NULL);
int a, b, p;
while (cin>>a>>p>>b && a) {
int res=exbsgs(a,b,p);
if (res==-1) cout<<"No Solution\n";
else cout<<res<<'\n';
}
return 0;
}
Record,\(2.46\rm{s}\),可以通过本题(包括Hack数据)。
后记
这道题算是比较毒瘤的,我是一共调了三天才过的(我太弱了)
还有\(20\)天就要中考了,而我还在这摸鱼(悲)
就谨此纪念一下我的初中OI生涯罢。
顺便在此祝福向宴酱中考顺利!