前言:
今天更BSGS算法。俗称大步小步算法(Big-Step G…-Step),又称拔山盖世、北上广深、白色狗屎 。
问题:
求解指数同余方程:
a
x
≡
b
(
m
o
d
p
)
a^{x}\equiv b(mod\; p)
ax≡b(modp)的最小自然数解。
BSGS算法:模板题
这个算法只考虑
g
c
d
(
a
,
p
)
=
1
gcd(a,p)=1
gcd(a,p)=1的情况。
那么,怎么做呢。
其实做法很简单,就是枚举。
先来一种暴力枚举。
就是枚举
x
x
x从
[
0
,
p
−
1
]
[0,p-1]
[0,p−1],然后检验答案是否正确。
知道欧拉定理的应该知道只要枚举到
φ
(
p
)
−
1
\varphi(p)-1
φ(p)−1,但
p
p
p为质数时欧拉函数仍为
p
−
1
p-1
p−1。
看过阶与原根的应该知道只要枚举到
p
p
p的阶即可,但算
p
p
p的阶仍然要暴力去枚举。
所以这里只要一个一个枚举过去便好。
但是
o
(
p
)
o(p)
o(p)的暴力枚举配不上这么高大上的名字,所以还有优化的余地。
我们可以考虑折半枚举。
怎么折半,就是这个算法的精髓了。
我们令
x
=
k
n
+
i
x=kn+i
x=kn+i,其中
n
n
n为我们设置的一个定值,
k
k
k和i
i
i
i均为变量。
这样原式就变为
a
k
n
+
i
≡
b
(
m
o
d
p
)
a^{kn+i}\equiv b(mod\;p)
akn+i≡b(modp)。
由于
g
c
d
(
a
,
p
)
=
1
gcd(a,p)=1
gcd(a,p)=1,我们可以将
a
i
a^{i}
ai的逆元乘到右边,变成
a
k
n
≡
b
(
a
i
)
−
1
(
m
o
d
p
)
a^{kn}\equiv b(a^{i})^{-1}(mod\;p)
akn≡b(ai)−1(modp)。
这样我们可以枚举所有的
i
i
i,将其对应右边的值存入哈希表(也可用unordered_map)。然后枚举
k
k
k的过程中,去表中找满足的
i
i
i,找到后
k
n
+
i
kn+i
kn+i就是答案。当然因为要找最小的,枚举过程中动点手脚就好了。
那么这个
n
n
n的设置嘛,利用分块的思想,
n
n
n设为
p
\sqrt{p}
p就好了,这样
k
k
k只用枚举到
p
\sqrt{p}
p,
i
i
i也只用枚举到
p
\sqrt{p}
p。总的复杂度为
n
\sqrt{n}
n。
上代码:
#include <bits/stdc++.h>
#include <unordered_map>
using namespace std;
typedef long long ll;
const ll mod=1e9+7;
ll fpow(ll a,ll n,ll mod)
{
ll sum=1,base=a%mod;
while(n!=0)
{
if(n%2)sum=sum*base%mod;
base=base*base%mod;
n/=2;
}
return sum;
}
ll ex_gcd(ll a,ll b,ll& x,ll& y)
{
if(b==0)
{
x=1;y=0;
return a;
}
ll ans=ex_gcd(b,a%b,x,y);
ll tmp=x;
x=y;
y=tmp-a/b*y;
return ans;
}
ll inv(ll a,ll mod)//存在逆元条件:gcd(a,mod)=1
{
ll x,y;
ll g=ex_gcd(a,mod,x,y);
if(g!=1)return -1;
return (x%mod+mod)%mod;
}
ll BSGS(ll a,ll b,ll p)
{
b%=p;
if(b==1||p==1)return 0;
ll n=sqrt(p);
static unordered_map<ll,ll>Bmp;
Bmp.clear();
ll inva=inv(fpow(a,n-1,p),p)*b%p;
for(ll i=n-1;i>=0;i--)
{
Bmp[inva]=i;inva=inva*a%p;
}
ll ta=1,powa=fpow(a,n,p);
for(ll k=0;k<=p;k+=n)
{
if(Bmp.count(ta))return k+Bmp[ta];
ta=ta*powa%p;
}
return -1;
}
int main()
{
ll p,b,n;
scanf("%lld%lld%lld",&p,&b,&n);
ll ans=BSGS(b,n,p);
if(ans==-1)printf("no solution\n");
else printf("%lld\n",ans);
return 0;
}
扩展BSGS:模板题
既然BSGS只能解决
g
c
d
(
a
,
p
)
=
1
gcd(a,p)=1
gcd(a,p)=1的情况,那么当然还有
g
c
d
(
a
,
p
)
gcd(a,p)
gcd(a,p)可能不等于
1
1
1的情况。
那么这个时候就需要对算法进行扩展了。
因为
g
c
d
(
a
,
p
)
=
1
gcd(a,p)=1
gcd(a,p)=1,所以不能再像原来一样逆元了。
所以我们考虑把
g
c
d
gcd
gcd先提取出来。
我们令
d
=
g
c
d
(
a
,
p
)
d=gcd(a,p)
d=gcd(a,p)。
我们将原式
a
x
≡
b
(
m
o
d
p
)
a^{x}\equiv b(mod\;p)
ax≡b(modp)展开成
a
x
+
k
p
=
b
a^{x}+kp=b
ax+kp=b。
然后同时除以
d
d
d变成
a
x
d
+
k
p
d
=
b
d
\frac{a^{x}}{d}+k\frac{p}{d}=\frac{b}{d}
dax+kdp=db。
再同时取模
p
d
\frac{p}{d}
dp就得到了
a
x
d
≡
b
d
(
m
o
d
p
d
)
\frac{a^{x}}{d}\equiv \frac{b}{d}(mod\;\frac{p}{d})
dax≡db(moddp)
容易发现
b
b
b必须是
d
d
d的倍数,否则直接返回无解。
然后此时
g
c
d
(
a
d
,
p
d
)
=
1
gcd(\frac{a}{d},\frac{p}{d})=1
gcd(da,dp)=1,所以可以乘上它的逆元。
变成
a
x
−
1
≡
b
d
(
a
d
)
−
1
(
m
o
d
p
d
)
a^{x-1}\equiv \frac{b}{d}(\frac{a}{d})^{-1}(mod\;\frac{p}{d})
ax−1≡db(da)−1(moddp)。
就变成了一个新的问题。
那么观察此时
g
c
d
(
a
,
p
d
)
gcd(a,\frac{p}{d})
gcd(a,dp)是否一定为
1
1
1。
当然不一定,举例
a
=
2
,
p
=
8
a=2,p=8
a=2,p=8,那么
d
=
2
d=2
d=2,将
p
p
p除以
d
d
d后,仍能除以
d
d
d,没事,那再除一次嘛。
笔者懒,写了递归版,这个过程可以一步到位的,递归可能会慢一些,但是比起固有
p
\sqrt{p}
p的复杂度,多个
l
o
g
log
log也无伤大雅。
那么递归的思路就很清晰了。
如果
g
c
d
(
a
,
p
)
=
1
gcd(a,p)=1
gcd(a,p)=1,那么用BSGS求解。
如果
g
c
d
(
a
,
p
)
≠
1
gcd(a,p)\neq 1
gcd(a,p)=1,那么用上面的办法将模
p
p
p的问题变为模
p
g
c
d
(
a
,
p
)
\frac{p}{gcd(a,p)}
gcd(a,p)p的问题,然后再递归调用这个函数。
非递归思路:
如果
g
c
d
(
a
,
p
)
=
1
gcd(a,p)=1
gcd(a,p)=1,那么用BSGS求解。
如果
g
c
d
(
a
,
p
)
≠
1
gcd(a,p)\neq 1
gcd(a,p)=1,那么用上面的方法给
p
p
p除以
g
c
d
(
a
,
p
)
gcd(a,p)
gcd(a,p),只要除以多少个弄清楚,就可以一次弄好了。
最后是我的递归版代码:
#include <bits/stdc++.h>
#include <unordered_map>
using namespace std;
typedef long long ll;
const ll mod=1e9+7;
ll fpow(ll a,ll n,ll mod)
{
ll sum=1,base=a%mod;
while(n!=0)
{
if(n%2)sum=sum*base%mod;
base=base*base%mod;
n/=2;
}
return sum;
}
ll ex_gcd(ll a,ll b,ll& x,ll& y)
{
if(b==0)
{
x=1;y=0;
return a;
}
ll ans=ex_gcd(b,a%b,x,y);
ll tmp=x;
x=y;
y=tmp-a/b*y;
return ans;
}
ll inv(ll a,ll mod)//存在逆元条件:gcd(a,mod)=1
{
ll x,y;
ll g=ex_gcd(a,mod,x,y);
if(g!=1)return -1;
return (x%mod+mod)%mod;
}
ll gcd(ll a,ll b)
{
return b?gcd(b,a%b):a;
}
ll BSGS(ll a,ll b,ll p)
{
b%=p;
if(b==1||p==1)return 0;
ll n=sqrt(p);
static unordered_map<ll,ll>Bmp;
Bmp.clear();
ll inva=inv(fpow(a,n-1,p),p)*b%p;
for(ll i=n-1;i>=0;i--)
{
Bmp[inva]=i;inva=inva*a%p;
}
ll ta=1,powa=fpow(a,n,p);
for(ll k=0;k<=p;k+=n)
{
if(Bmp.count(ta))return k+Bmp[ta];
ta=ta*powa%p;
}
return -1;
}
ll exBSGS(ll a,ll b,ll p)
{
b%=p;
if(a==0&&b==0)return 1;
else if(a==0&&b!=0)return -1;
if(b==1||p==1)return 0;
ll d=gcd(a,p);
if(b%d!=0)return -1;
p=p/d;
b=b/d*inv(a/d,p)%p;
if(d!=1)
{
ll ans=exBSGS(a,b,p);
if(ans==-1)return -1;
return ans+1;
}
ll ans=BSGS(a,b,p);
if(ans==-1)return -1;
return ans+1;
}
int main()
{
ll a,p,b;
while(scanf("%lld%lld%lld",&a,&p,&b)!=EOF)
{
if(a==0&&p==0&&b==0)break;
ll ans=exBSGS(a,b,p);
if(ans==-1)printf("No Solution\n");
else printf("%lld\n",ans);
}
return 0;
}
最后再补充一个例题,Mod Tree。注意题目中的取模后相同并非模意义下相同。
高次同余方程:
不知道原根的出门左转。
最后再将一个BSGS的应用。
求解方程
x
n
≡
a
(
m
o
d
p
)
x^{n}\equiv a(mod\;p)
xn≡a(modp)。
我们只考虑
p
p
p为素数的求解办法。
做法:
我们先求出
p
p
p的原根
g
g
g。
根据原根的性质,当
p
p
p为素数时,
g
g
g的幂次可以得到一个完全剩余系。
所以我们可以令
x
=
g
k
x=g^{k}
x=gk,同样可以令
a
=
g
x
0
a=g^{x_{0}}
a=gx0。
此时,我们只要求解
k
k
k,就可以求解
x
x
x了。
方程变为
g
k
n
≡
g
x
0
(
m
o
d
p
)
g^{kn}\equiv g^{x_{0}}(mod\;p)
gkn≡gx0(modp)。
等价于方程
k
n
≡
x
0
(
m
o
d
p
−
1
)
kn\equiv x_{0}(mod\;p-1)
kn≡x0(modp−1)。(这个可以根据费马小定理,也可以根据
o
r
d
p
g
=
p
−
1
ord_{p}g=p-1
ordpg=p−1)
因为不能确保
g
c
d
(
n
,
p
−
1
)
=
1
gcd(n,p-1)=1
gcd(n,p−1)=1,所以不能乘逆元。
其实方程等价于求方程
n
k
+
(
p
−
1
)
y
=
x
0
nk+(p-1)y=x_{0}
nk+(p−1)y=x0,用扩展欧几里得即可。
到这一步答案就解出来了。
嗯?我要说的BSGS呢?
噢,你会发现这个
x
0
x_{0}
x0还没求。
g
x
0
≡
a
(
m
o
d
p
)
g^{x_{0}}\equiv a(mod\;p)
gx0≡a(modp),所以求
x
0
x_{0}
x0就是求指数同余方程了。
到这里,求毕。
这个没写代码,我懒了,想都别想。
那么关于
p
p
p不为素数的。
其实我不会,但是我还可以鬼扯,可能可行也可能不行,往下阅读时请保留自己的理智,免得被笔者误导。
考虑分解合数
p
=
p
1
k
1
p
2
k
2
…
p
n
k
n
p=p_{1}^{k_{1}}p_{2}^{k_{2}}…p_{n}^{k_{n}}
p=p1k1p2k2…pnkn。
噢,好像不可行。
但是当
p
=
p
1
p
2
…
p
n
p=p_{1}p_{2}…p_{n}
p=p1p2…pn,即每个素数的次数都为
1
1
1次时,还是有办法可行的。
我们只要求解出对于每个模数的方程
x
n
≡
a
(
m
o
d
p
i
)
x^{n}\equiv a(mod\;p_{i})
xn≡a(modpi)的解
x
≡
a
0
(
m
o
d
p
)
x\equiv a_{0}(mod\;p)
x≡a0(modp)。
然后用中国剩余定理合并即可。
这样想的话,如果有办法求解
x
n
≡
a
(
m
o
d
p
i
k
i
)
x^{n}\equiv a(mod\;p_{i}^{k_{i}})
xn≡a(modpiki)的话,就可以解决任意模数的问题了,暂时还没想好。
扯毕。
补充:
新学的。。源自于网安RSA加解密算法。
补充一个解方程
x
n
≡
a
(
m
o
d
p
)
x^{n}\equiv a(mod\;p)
xn≡a(modp)在满足
g
c
d
(
n
,
φ
(
p
)
)
=
1
gcd(n,\varphi (p))=1
gcd(n,φ(p))=1条件下的做法。
我们在方程两边都开一个幂次
b
b
b,得到
x
n
b
≡
a
b
(
m
o
d
p
)
x^{nb}\equiv a^{b}(mod\;p)
xnb≡ab(modp)。
发现当
n
b
≡
1
(
m
o
d
φ
(
p
)
)
nb\equiv 1(mod\;\varphi(p))
nb≡1(modφ(p))时,
a
b
a^{b}
ab就是方程的一个解。
所以,我们只要找到一个
b
b
b,使得
n
b
≡
1
(
m
o
d
φ
(
p
)
)
nb\equiv 1(mod\;\varphi(p))
nb≡1(modφ(p))成立,那么问题就变成了拓展欧几里得的经典问题了。
这样就可以在
o
(
l
o
g
p
)
o(logp)
o(logp)的时间复杂度下解决,比起
o
(
p
)
o(\sqrt{p})
o(p)的做法真是优秀太多了,当然前提要满足
g
c
d
(
n
,
φ
(
p
)
)
=
1
gcd(n,\varphi (p))=1
gcd(n,φ(p))=1,这个条件也是来自于拓欧方程有解的充要条件。