[转载][学习笔记]大步小步法解决离散对数问题

注：本文除代码以外，其他部分均转载于这位大佬

离散对数（Discrete Logarithm）问题是这样一个问题，它是要求解模方程

a^{x} \equiv b (\mod m)

$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 ⌈ \sqrt{p} ⌉ + B} \equiv b (\mod p)

$a^{A\lceil \sqrt{p} \rceil + B} \equiv b \pmod p$

可以在两边同时乘以 $a^B$ 的逆元，由于 $p$ 是质数，这个逆元一定存在，于是方程变成

a^{A ⌈ \sqrt{p} ⌉} \equiv b \cdot a^{- B} (\mod 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 ⌈ \sqrt{p} ⌉} \equiv b \cdot a^{B} (\mod p)

$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} + k m = b, k \in Z

$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 Z

$\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} (\mod \frac{m}{g})

$\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^{'} (\mod m^{'})

$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;
}

posted @ 2020-05-28 16:52 Arextre 阅读(751) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

阅读排行：
· 分享一个免费、快速、无限量使用的满血 DeepSeek R1 模型，支持深度思考和联网搜索！
· 基于 Docker 搭建 FRP 内网穿透开源项目（很简单哒）
· ollama系列01：轻松3步本地部署deepseek，普通电脑可用
· 25岁的心里话
· 按钮权限的设计及实现

公告

这位身穿漆黑长袍且头戴三角帽，胸口上别着象征星辰的胸针，灰色的发丝随风摇曳，在太阳照耀下散发出耀眼的光芒，琉璃色的双瞳看似朝向前方，实际上却是眺望着远方的某处。她所经之处树木挥动枝条发出低语，挥洒树叶为她献上祝福。这位具有如此美貌，任谁都只能以惹人怜爱形容的美少女究竟是谁呢？
我不知道诶，我只知道博客的密码是她的生辰~~~
另外，背景的小伊蕾娜是 Yuuki 老师画的，祂还很用心地发了我一份没有水印的版本，以作为我的博客背景~~~祂的 pixiv uid 是 36481656 .

昵称： Arextre
园龄： 5年2个月
粉丝： 38
关注： 24

+加关注

2025年3月

日

一

二

三

四

五

六

Arextre

[转载][学习笔记]大步小步法解决离散对数问题

公告

搜索

常用链接

最新随笔

积分与排名

随笔分类 (982)

随笔档案 (378)

只能膜拜の大佬们

阅读排行榜

评论排行榜

推荐排行榜

最新评论