非常不错的一篇博客 http://blog.miskcoo.com/2015/05/discrete-logarithm-problem
baby step giant step
用来解
一类的方程 O(p–√) O ( p )
其实想法也比较简单,对于a,p互质的情况:
定一个值
m=p–√
m
=
p
,x可以写为%m下的表达形式:
Am+r
A
m
+
r
。
把方程写为
aAm≡b⋅(ar)−1(mod p)
a
A
m
≡
b
⋅
(
a
r
)
−
1
(
m
o
d
p
)
小步: 因为A与r的范围都是
p–√
p
的,先预处理出所有
b⋅(ar)−1
b
⋅
(
a
r
)
−
1
的值,hash存下。
大步: 从0到m枚举A,
若
aAm=b⋅(ar)−1,r∈[0,m)
a
A
m
=
b
⋅
(
a
r
)
−
1
,
r
∈
[
0
,
m
)
,说明存在
x=Am+r
x
=
A
m
+
r
。
若需要最小整数解则r应该取最小值。
一个小优化
可以将
Am+r
A
m
+
r
的改为
Am−r
A
m
−
r
的形式,这样不需要处理逆元,但是还是需要用到逆元的存在!
若
aAm=b⋅(ar),r∈[0,m)
a
A
m
=
b
⋅
(
a
r
)
,
r
∈
[
0
,
m
)
,说明存在
x=Am−r
x
=
A
m
−
r
。
注意此时A的枚举范围从0..m改为1..m+1。
此时需要特判b=1的情况,否则因为A是从1开始枚举的,而
r<m
r
<
m
,这样是发现不了x=0的解的。
扩展BSGS
其实也不是什么扩展,只是处理了a,p不互质的情况
记
g=(a,p)
g
=
(
a
,
p
)
将
ax≡b(mod p)
a
x
≡
b
(
m
o
d
p
)
写为等式,即
a⋅ax−1−kp=b
a
⋅
a
x
−
1
−
k
p
=
b
由exgcd知识,b不是g的倍数时无解。
等式两边同除以g,可以得到
若
(a,p/g)=1
(
a
,
p
/
g
)
=
1
,这种形式也可以用bsgs解,方法是显然的,甚至不用处理逆元,只需要将大步与小步的初始值改变即可。
最后将解加一就是原方程的解了。
注意并不是将(a,p)变为(a/g,p/g),所以可能上述过程需要做多次,将a/g的积处理一下即可。
注意为了防止最终要解的方程出现负指数,上述过程中注意判断x=0的解。
具体实现见标吧,上面推荐的那篇博客。
Note
mod 1的情况需要特判
b=1的情况需要特判
exbsgs过程中注意判断x=0的解
下面是给自己参阅(co)的标
值得一提的是手打hash要比map快上1/2左右
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <map>
using namespace std;
typedef long long ll;
ll y,z,p;
ll ksm(ll x,ll y,ll m) {
ll ret=1;
for (; y; y>>=1) {
if (y&1) ret=ret*x%m;
x=x*x%m;
}
return ret;
}
ll gcd(ll x,ll y) {return y==0?x:gcd(y,x%y);}
map<ll,ll> ha;
ll exbsgs(ll a,ll b,ll p) {
a%=p,b%=p;
if (b==1||p==1) return 0; //1
ll cnt=0;
ll t=1;
for (ll g=gcd(a,p); g!=1; g=gcd(a,p)) {
if (b%g) return -1;
p/=g,b/=g,t=t*(a/g)%p;
cnt++;
if (b==t) return cnt;//2
}
ha.clear();
ll m=sqrt(p)+1;
for (ll i=0,s=b; i<m; i++,s=s*a%p) {
ha[s]=i;
}
ll step=ksm(a,m,p);
for (ll i=1; i<=m+1; i++) {
t=t*step%p;
if (ha.count(t)) return cnt + i*m - ha[t];
}
return -1;
}
int main() {
cin>>y>>p>>z;
for (; ; cin>>y>>p>>z) {
if (y==0 && z==0 && p==0) return 0;
ll ret=exbsgs(y,z,p);
if (ret==-1) cout<<"No Solution"<<endl;
else cout<<ret<<endl;
}
}