大步小步算法 BSGS(baby steps giant steps)
BSGS是用来解决离散对数问题的,即
a
x
≡
b
(
m
o
d
p
)
a^x \equiv b \pmod p
ax≡b(modp)。其中,a,b,p已知,且a和p互质,求x。根据欧拉定理,我们知道
a
ϕ
(
p
)
≡
1
(
m
o
d
p
)
a^{\phi(p)} \equiv 1 \pmod p
aϕ(p)≡1(modp),而
ϕ
(
p
)
<
p
\phi(p)<p
ϕ(p)<p,所以,采用枚举法,我们可以在
O
(
p
)
O(p)
O(p)时间复杂度求出x。
BSGS可以在
O
(
p
)
O(\sqrt p)
O(p)的时间复杂度内求出x,它是这么做的:
首先,令
m
=
⌈
p
⌉
,
r
=
x
m
o
d
m
m= \lceil \sqrt p \rceil ,\; r=x \bmod m
m=⌈p⌉,r=xmodm则
x
=
k
∗
m
+
r
(1)
x=k*m+r \tag 1
x=k∗m+r(1)
其中
0
≤
k
<
m
,
0
≤
r
<
m
0\leq k < m ,0 \leq r < m
0≤k<m,0≤r<m
于是有
a
k
∗
m
+
r
≡
b
(
m
o
d
p
)
a^{k*m+r} \equiv b \pmod p
ak∗m+r≡b(modp)
因为a与p互质,所以两边同时乘以
a
−
r
a^{-r}
a−r,则
a
k
∗
m
≡
b
∗
a
−
r
(
m
o
d
p
)
a^{k*m} \equiv b*a^{-r} \pmod p
ak∗m≡b∗a−r(modp)
对于
(
0
≤
i
≤
r
)
(0 \leq i \leq r)
(0≤i≤r),求出所有的
b
∗
a
−
i
m
o
d
p
b*a^{-i} \bmod p
b∗a−imodp,将值和对应的i都保存在map中.
然后再求左边
a
j
∗
m
m
o
d
p
(
0
≤
j
<
k
)
a^{j*m} \bmod p (0 \leq j < k)
aj∗mmodp(0≤j<k) ,并在map中查找是否出现过相同的值。如果有,则已经找出答案了,则答案为
j
∗
m
+
i
j*m+i
j∗m+i。如果没有,表示无解。
其实我们可以稍微改进一下,以避免求逆元。
将公式(1)变换一下:
x
=
k
∗
m
−
r
(2)
x=k*m-r \tag 2
x=k∗m−r(2)
此时,
1
≤
k
≤
m
,
0
<
r
≤
m
1 \leq k \leq m,0 < r \leq m
1≤k≤m,0<r≤m
则
a
k
∗
m
−
r
≡
b
(
m
o
d
p
)
(3)
a^{k*m-r} \equiv b \pmod p \tag 3
ak∗m−r≡b(modp)(3)
两边同时乘以
a
r
a^r
ar,可得:
a
k
∗
m
≡
b
∗
a
r
(
m
o
d
p
)
(4)
a^{k*m} \equiv b*a^r \pmod p \tag 4
ak∗m≡b∗ar(modp)(4)
先求出右边所有的
b
∗
a
i
(
m
o
d
p
)
(
1
≤
i
≤
r
)
b*a^i (\bmod p) (1 \leq i \leq r)
b∗ai(modp)(1≤i≤r) ,保存在map中,然后再求左边的
a
j
∗
m
m
o
d
p
a^{j*m} \bmod p
aj∗mmodp,并在map中查找是否出现过。如果出现过,左边枚举的是j,右边枚举的是i,则答案为
x
=
j
∗
m
−
i
x=j*m-i
x=j∗m−i.这样就可以不用求逆元了,但仍然用到了逆元。因为我们的推导必须是等价推导,只有当a与p互质时,即
a
r
a^r
ar的逆元存在时,式(3)和式(4)才是等价的。因为如果a不与p互质,式(4)是不能推出式(3)的。
例题1:(poj2417)
题意:给出一个质数P,一个整数B,2<=B<P,和一个整数N,1<=N<=P,求底为B,模为P下的离散对数,如果有多个解,则输出最小的非负整数解。即解方程
B
x
≡
N
m
o
d
P
B^x \equiv N \bmod P
Bx≡NmodP,求最小的x,x>=0
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<map>
#include<cmath>
using namespace std;
#define LL long long int
map<LL,int> mp;
int b,n,p,x,ans;
int bsgs()
{
mp.clear();
int m=ceil(sqrt(p));
LL tmp=1;
for(int i=1;i<=m;i++)
{
tmp=tmp*b%p;
mp[tmp*n%p]=i; //这里会覆盖掉之前的解,即保存的是最大的i。
}
int res=1;
for(int i=1;i<=m;i++)
{ res=res*tmp%p;
if(mp[res])
{
ans=m*i-mp[res]; //因为mp[res]保存的是最大的,所以ans是最小的解。
return ans;
}
}
return -1;
}
int main()
{
while(scanf("%d%d%d",&p,&b,&n)!=-1)
{
n%=p,b%=p;
if(n==1)
ans=0;
else if(b==0)
{if(n!=0)
ans=-1;
else ans=1;
}
else
ans=bsgs();
if(ans>=0)printf("%d\n",ans);
else printf("no solution\n");
}
return 0;
}
扩展BSGS
对于 a x ≡ b ( m o d p ) (1) a^x \equiv b \pmod p \tag 1 ax≡b(modp)(1)
如果a与p不互质,则不能直接套用BSGS。而扩展BSGS正是解决这种情况的。
可以将(1)变形:
a
x
+
k
p
=
b
(2)
a^x+kp=b \tag 2
ax+kp=b(2)
设a与p的最大公约数为g,即
g
=
(
a
,
p
)
g=\left(a,p\right)
g=(a,p).若
g
∤
b
g \nmid b
g∤b,则方程无解,不过这里有一个例外,即当b=1时,x=0。
将式(2)左右两边除以g,得到:
a
x
−
1
a
g
+
k
p
g
=
b
g
(3)
a^{x-1}\frac{a}{g}+k\frac{p}{g}=\frac{b}{g} \tag 3
ax−1ga+kgp=gb(3)
设
a
′
=
a
g
,
p
′
=
p
g
a'=\frac{a}{g},p'=\frac{p}{g}
a′=ga,p′=gp,
b
′
=
b
g
b'=\frac{b}{g}
b′=gb,
a
x
−
1
a
′
+
k
p
′
=
b
′
a^{x-1}a'+kp'=b'
ax−1a′+kp′=b′
若a与
p
g
\frac{p}{g}
gp仍然不互质,则继续找出a与p’的最大公约数g’,(3)两边继续除以g’,直到最后a与p’互质为止。
这里我们假设a与p’互质,此时,a’与p’也是互质的,于是可以转换为:
a
x
−
c
n
t
≡
b
′
(
a
′
)
−
1
(
m
o
d
p
′
)
a^{x-cnt} \equiv b'(a')^{-1}\pmod{p'}
ax−cnt≡b′(a′)−1(modp′)
其中cnt表示两边除以最大公约数g的次数。
此处右边有逆元。为了避免求逆元,可以将a’保留在左边。在用bsgs枚举左边时,初始值就设为a’。
如果在除以g的过程中,发现
b
′
(
a
′
)
−
1
=
1
b'(a')^{-1}=1
b′(a′)−1=1,则立马可以得到答案,即x-cnt=0,即x=cnt。
此时,直接套用基础的BSGS即可了。
不要忘了最后的答案要加上cnt。
#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cstring>
#include<map>
#include<cmath>
#include<vector>
using namespace std;
#define LL long long int
#define HASHL 10007
int x,z,k,m,ans;
int hashcode(int a)
{
a=(a*998244353LL+100007)%HASHL;
return a%HASHL;
}
vector<pair<int,int> > myv[HASHL+5];
int getphi(int x)
{
int res=x;
for(int i=2;i*i<=x;i++)
{
if(x%i==0)
{
while(x%i==0)
x/=i;
}
res=res*(i-1)/i;
}
if(x>1)
res=res*(x-1)/x;
return res;
}
bool flg[HASHL+5];
int gcd(int a,int b)
{
if(!b)return a;
else return gcd(b,a%b);
}
int exbsgs()
{
memset(flg,0,sizeof flg);
if(z==0)return -1;
if(z==1)
{
if(k==0)return 0;
else return -1;
}
if(k==1)return 0;
x%=z,k%=z;
int g=1,cnt=0,tmpk=1;
while((g=gcd(x,z))>1)
{
if(tmpk%z==k%z)return cnt;
if(k%g)return -1;
k=k/g;
z=z/g;
tmpk=1ll*tmpk*(x/g)%z;
cnt++;
}
int m=ceil(sqrt(double(z)));
int tmp=1,tmpid=0;
for(int i=1;i<=m;i++)
{
tmp=1ll*tmp*x%z;
tmpid=hashcode(1ll*tmp*k%z);
if(flg[tmpid]==0)
{
myv[tmpid].clear();
flg[tmpid]=1;
}
myv[tmpid].push_back(make_pair(1ll*tmp*k%z,i));
}
int res=tmpk%z,tmpres=1;
for(int i=1;i<=m;i++)
{
res=1ll*res*tmp%z;
tmpres=hashcode(res);
if(flg[tmpres]==1)
{
for(int j=myv[tmpres].size()-1;j>=0;j--)
{
if(myv[tmpres][j].first==res)
{
return (i*m-myv[tmpres][j].second)+cnt;
}
}
}
}
return -1;
}
int main()
{
int ans;
while(scanf("%d %d %d",&x,&z,&k)!=-1)
{
if(!(x||z||k))return 0;
ans=exbsgs();
if(ans==-1)printf("No Solution\n");
else printf("%d\n",ans);
}
return 0;
}