[转载][学习笔记]大步小步法解决离散对数问题
注：本文除代码以外，其他部分均转载于 这位大佬
离散对数（Discrete Logarithm）问题是这样一个问题，它是要求解模方程
\[ a^x \equiv b \pmod m
\]
这个问题是否存在多项式算法目前还是未知的，这篇文章先从 \(m\) 是质数开始介绍大步小步法（Baby Step Giant Step）来解决它，之后再将其应用到 \(m\) 是任意数的情况。这个算法可以在 \(\mathcal O(\sqrt m)\) 的时间内计算出最小的 \(x\) ，或者说明不存在这样一个 \(x\) 。
题目链接：BZOJ-2480、SPOJ-MOD、BZOJ-3239
首先解决 \(m=p\) 是质数的情况，我们可以设 \(x = A \lceil \sqrt{p} \rceil + B\) ，其中 \(0 \leq B < \lceil \sqrt{p} \rceil\), \(0 \leq A < \lceil \sqrt{p} \rceil\) ，这样的话就变成求解 \(A\) 和 \(B\) 了，方程也变成
\[a^{A\lceil \sqrt{p} \rceil + B} \equiv b \pmod p
\]
可以在两边同时乘以 \(a^B\) 的逆元，由于 \(p\) 是质数，这个逆元一定存在，于是方程变成
\[ a^{A\lceil \sqrt{p} \rceil} \equiv b\cdot a^{-B} \pmod p
\]
由于 \(A, B\) 都是 \(\mathcal O(\sqrt p)\) 级别的数，可以先计算出右边这部分的值，存入 Hash 表，然后计算左边的值，在 Hash 表中查找，只要按照从小到大的顺序如果有解就能够找到最小的解，由于两边都只有 \(\mathcal O(\sqrt p)\) 个数，因此时间复杂度是 \(\mathcal O(\sqrt p)\) 的，这样 \(m\) 是质数的情况就解决了
一个优化：我们可以设 \(x = A \lceil \sqrt{p} \rceil - B\) ，其中 \(0 \leq B < \lceil \sqrt{p} \rceil\) , \(0 < A \leq \lceil \sqrt{p} \rceil + 1\) ，这样的话化简后的方程就是
\[a^{A\lceil \sqrt{p} \rceil} \equiv b\cdot a^B \pmod p
\]
就可以不用求出逆元，要注意只是不用求出逆元，而不是没有用到逆元的存在
现在来看 \(m\) 不是质数的情况，同样可以设 \(x = A \lceil \sqrt{m} \rceil + B\) ，根据上面的推导，会发现需要用到的性质就是 \(a^B\) 的逆元存在，所以当 \(m\) 和 \(a\) 互质的时候这个方法仍然有效！
如果 \((m, a) \neq 1\) 该怎么办呢？我们要想办法把方程转化为 \((m, a) = 1\) 的情况
把要求的模方程写成另外一种形式
\[a^x + km = b, k \in \mathbb Z
\]
设 \(g = (a, m)\) ，这样的话可以确定如果 \(g \nmid b\) 那么该方程一定无解，所以当 \(g \mid b\) 的时候，在方程左右两边同时除以 \(g\)
\[\frac{a}{g} a^{x-1} + k \frac{m}{g} = \frac{b}{g}, k \in \mathbb{Z}
\]
这样便消去了一个因子，得到方程
\[ \frac{a}{g} a^{x-1} \equiv \frac{b}{g} \pmod {\frac{m}{g}}
\]
令 \(m' = \frac{m}{g}, b' = \frac{b}{g} \left (\frac{a}{g}\right)^{-1}\) （这里不可以把 \(g\) 消掉），就可以得到新的方程
\[ a^{x'} \equiv b' \pmod {m'}
\]
得到解之后原方程的解 \(x = x' + 1\) ，不断重复这个过程最后一定会得到一个可以解的方程，套用刚刚的大步小步法解出后即可。要注意的是在这个过程中如果某一步发现 \(b'=1\) ，那么就可以直接退出，因为这时候已经得到了解
NOTE：上面这个过程是可能执行多次的比如说 \((a, m) = (6, 16) \rightarrow (6, 8) \rightarrow (6, 4) \rightarrow (6, 2)\)
代码
- 普通 bsgs
#include<bits/stdc++.h>
using namespace std;
#define rep(i,__l,__r) for(signed i=(__l),i##_end_=(__r);i<=i##_end_;++i)
#define fep(i,__l,__r) for(signed i=(__l),i##_end_=(__r);i>=i##_end_;--i)
#define writc(a,b) fwrit(a),putchar(b)
#define mp(a,b) make_pair(a,b)
#define ft first
#define sd second
#define LL long long
#define ull unsigned long long
#define uint unsigned int
#define pii pair< int,int >
#define Endl putchar('\n')
// #define int long long
// #define int unsigned
// #define int unsigned long long
#ifdef _GLIBCXX_CSTDIO
#define cg (c=getchar())
template<class T>inline void qread(T& x){
char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
if(f)x=-x;
}
template<class T>inline T qread(const T sample){
T x=0;char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
return f?-x:x;
}
#undef cg
template<class T>void fwrit(const T x){//just short,int and long long
if(x<0)return (void)(putchar('-'),fwrit(-x));
if(x>9)fwrit(x/10);
putchar(x%10^48);
}
#endif
// template<class T,class... Args>inline void qread(T& x,Args&... args){qread(x),qread(args...);}
template<class T>inline T Max(const T x,const T y){return x>y?x:y;}
template<class T>inline T Min(const T x,const T y){return x<y?x:y;}
template<class T>inline T fab(const T x){return x>0?x:-x;}
inline int gcd(const int a,const int b){return b?gcd(b,a%b):a;}
inline void getInv(int inv[],const int lim,const int MOD){
inv[0]=inv[1]=1;for(int i=2;i<=lim;++i)inv[i]=1ll*inv[MOD%i]*(MOD-MOD/i)%MOD;
}
inline LL mulMod(const LL a,const LL b,const LL mod){//long long multiplie_mod
return ((a*b-(LL)((long double)a/mod*b+1e-8)*mod)%mod+mod)%mod;
}
const int sqMAXP=1<<17;
inline int qkpow(int a,int n,const int mod){
int ret=1;
for(;n>0;n>>=1){
if(n&1)ret=1ll*ret*a%mod;
a=1ll*a*a%mod;
}return ret;
}
map<int,int>hsh;
inline int bsgs(const int a,const int b,const int p){//2 3 5
if(a%p==0 && b!=0)return -1;
if(b==1)return 0;
int sqp=ceil(sqrt(p*1.0)),base=1;
hsh.clear();
for(int i=1;i<=sqp;++i){
base=1ll*base*a%p;
hsh[base]=i;
}
int now=1,tmp;
rep(i,0,sqp-1){
tmp=1ll*b*qkpow(now,p-2,p)%p;
if(hsh.count(tmp))return i*sqp+hsh[tmp];
now=1ll*now*base%p;//放在循环的最后, 需要注意一下
}return -1;
}
int T,a,b,p;
signed main(){
ios::sync_with_stdio(false);
cin>>T;
while(T--){cin>>p>>a>>b;
int ret=bsgs(a,b,p);
if(ret==-1)puts("Couldn't Produce!");
else printf("%d\n",ret);
}
return 0;
}
- 扩展 bsgs
写法一
#include<bits/stdc++.h>
using namespace std;
#define rep(i,__l,__r) for(signed i=(__l),i##_end_=(__r);i<=i##_end_;++i)
#define fep(i,__l,__r) for(signed i=(__l),i##_end_=(__r);i>=i##_end_;--i)
#define writc(a,b) fwrit(a),putchar(b)
#define mp(a,b) make_pair(a,b)
#define ft first
#define sd second
#define LL long long
#define ull unsigned long long
#define uint unsigned int
#define pii pair< int,int >
#define Endl putchar('\n')
// #define int long long
// #define int unsigned
// #define int unsigned long long
#ifdef _GLIBCXX_CSTDIO
#define cg (c=getchar())
template<class T>inline void qread(T& x){
char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
if(f)x=-x;
}
template<class T>inline T qread(const T sample){
T x=0;char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
return f?-x:x;
}
#undef cg
template<class T>void fwrit(const T x){//just short,int and long long
if(x<0)return (void)(putchar('-'),fwrit(-x));
if(x>9)fwrit(x/10);
putchar(x%10^48);
}
#endif
// template<class T,class... Args>inline void qread(T& x,Args&... args){qread(x),qread(args...);}
template<class T>inline T Max(const T x,const T y){return x>y?x:y;}
template<class T>inline T Min(const T x,const T y){return x<y?x:y;}
template<class T>inline T fab(const T x){return x>0?x:-x;}
inline int gcd(const int a,const int b){return b?gcd(b,a%b):a;}
inline void getInv(int inv[],const int lim,const int MOD){
inv[0]=inv[1]=1;for(int i=2;i<=lim;++i)inv[i]=1ll*inv[MOD%i]*(MOD-MOD/i)%MOD;
}
inline LL mulMod(const LL a,const LL b,const LL mod){//long long multiplie_mod
return ((a*b-(LL)((long double)a/mod*b+1e-8)*mod)%mod+mod)%mod;
}
int a,p,b;
void exgcd(const int a,const int b,int& x,int& y){
if(!b)x=1,y=0;
else{
exgcd(b,a%b,y,x);
y-=(a/b)*x;
}
}
inline int inv(const int a,const int mod){
int x,y;
exgcd(a,mod,x,y);
if(x<0)x+=mod;
return x;
}
inline int bsgs(const int a,const int b,const int c){
if(b==1)return 0;
if(a%p==0 && b!=0)return -1;
int sqc=ceil(sqrt(c*1.0)),base=1;
map<int,int>hsh;hsh.clear();
rep(i,1,sqc){
base=1ll*base*a%c;
if(!hsh.count(base))hsh[base]=i;
//注意当余数相同的时候出现的冲突, 需要重新判断一下
}
int now=1,tmp;
rep(i,0,sqc){
tmp=1ll*b*inv(now,c)%c;
if(hsh.count(tmp))return i*sqc+hsh[tmp];
now=1ll*now*base%c;
}return -1;
}
inline int exbsgs(int a,int b,int c){
a%=c,b%=c;
int cnt=0,g,base=1;
while((g=gcd(a,c))!=1){
if(b%g)return -1;
b/=g,c/=g,++cnt;
b=1ll*b*inv(a/g,c)%c;
a%=c;
if(b==1)return cnt;
}
int ret=bsgs(a,b,c);
if(ret==-1)return -1;
return ret+cnt;
}
signed main(){
ios::sync_with_stdio(false);
return 0;
}
写法二
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long LL;
LL read(){
LL x,f=1;char c;
while((c=getchar())<'0'||'9'<c)
if(c=='-')f=-1;
for(x=(c^48);'0'<=(c=getchar())&&c<='9';x=(x<<3)+(x<<1)+(c^48));
return x*f;
}
const LL MAXN=1e5+5;
LL gcd(LL a,LL b){
return b?gcd(b,a%b):a;
}
struct node{
LL i,r;
node(){
}
node(LL I,LL R):i(I),r(R){
}
friend bool operator < (const node& a,const node& b){
if(a.r==b.r)return a.i>b.i;
return a.r<b.r;
}
}lis[MAXN];
bool cmpNode_for_lower_bound(const node& a,const node& b){
return a.r<b.r;
}
bool check(LL a,LL b,LL c){
a%=c,b%=c;
if(b==1)return printf("0\n");
LL cnt=0,base=1;
for(LL g=gcd(a,c);g!=1;g=gcd(a,c)){
if(b%g)return 0;
b/=g,c/=g;
base=base*(a/g)%c,++cnt;
a%=c;
if(b==base)return printf("%lld\n",cnt);
}
LL m=ceil(sqrt(c)+1),tmp=1;
for(LL i=1;i<=m;++i){
tmp=tmp*a%c;
lis[i]=node(i,tmp*b%c);
}
sort(lis+1,lis+1+m);
for(LL i=1;i<=m;++i){
base=base*tmp%c;
node *p=lower_bound(lis+1,lis+m+1,node(0,base-1)/*,cmpNode_for_lower_bound*/);
if(p!=lis+m+1&&(*p).r==base)return printf("%lld\n",i*m-(*p).i+cnt);
}
return 0;
}
int main(){
LL a,b,c;
while(~scanf("%lld%lld%lld",&a,&c,&b)&&a&&b&&c){
// if(b>=c)puts("Orz,I can’t find D!");
if(!check(a,b,c))puts("Orz,I can’t find D!");
}
return 0;
}
