离散对数问题
给出
a
,
b
,
p
a,b,p
a,b,p，求解
a
x
=
b
(
m
o
d
p
)
a^x=b\ (mod\ \ p)
ax=b (mod p),其中
p
p
p是大素数,且p-1是光滑数
光滑数：素因子都是小素数，例如素因子都小于7的叫7-smooth数
Pohlig-Hellman
该算法常用于密码学。我在国内网上真心找不到该算法。有也是用来求椭圆问题的。所以来学习总结一下。
水平有限，英语还不好，国内还没看见讲解。只能在fanqiang学习。英语还差，心态简直炸裂。
给一个样例网站Example of the Pohlig-Hellman，虽然不知道国内能不能进。
因为BSGS的复杂度是 p \sqrt{p} p，所以当p是一个大素数时，就无法在限定时间内解决问题。也就用到了该算法。我们令 p − 1 = p 1 e 1 p 2 e 2 ⋯ p k e k p-1=p_1^{e_1}p_2^{e_2}\cdots p_k^{e_k} p−1=p1e1p2e2⋯pkek,算法复杂度大约为 O ( ∑ i = 1 k e i ( l o g n + p i ) ) O(\displaystyle\sum_{i=1}^{k}e_i(log\ {n}+\sqrt{p_i})) O(i=1∑kei(log n+pi)),最坏的情况下时间也就和 B S G S BSGS BSGS差不多。
首先找到
a
a
a相对于
p
p
p的阶
n
n
n（阶不一定是p-1，直接用p-1会出错）
令
n
=
p
1
e
1
p
2
e
2
⋯
p
k
e
k
n=p_1^{e_1}p_2^{e_2}\cdots p_k^{e_k}
n=p1e1p2e2⋯pkek
能得到
k
k
k组
x
=
x
i
(
m
o
d
p
i
e
i
)
x=x_i\ (mod\ \ p_i^{e_i})
x=xi (mod piei)，然后用中国剩余定理(CRT)合并求出
x
(
m
o
d
n
)
x\ (mod\ \ n)
x (mod n)
对于
x
i
x_i
xi:
我们先知道迭代函数：
对于
N
=
p
1
e
1
p
2
e
2
⋯
p
k
e
k
N=p_1^{e_1}p_2^{e_2}\cdots p_k^{e_k}
N=p1e1p2e2⋯pkek
t
0
,
t
1
⋯
t_0,t_1\cdots
t0,t1⋯是一个大于1的素数序列
{
t
0
=
⋯
=
t
e
1
−
1
=
p
1
t
e
1
=
⋯
=
t
e
1
+
e
2
−
1
=
p
2
⋮
t
u
−
e
k
+
1
=
⋯
=
t
u
−
1
=
p
k
(1)
\begin{cases}t_0=\cdots=t_{e_1-1}=p_1\\ t_{e_1}=\cdots=t_{e_1+e_2-1}=p_2\\ \vdots\\ t_{u-e_k+1}=\cdots=t_{u-1}=p_k \end{cases}\tag{1}
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧t0=⋯=te1−1=p1te1=⋯=te1+e2−1=p2⋮tu−ek+1=⋯=tu−1=pk(1)
从而
N
=
t
0
t
1
⋯
t
u
−
1
N=t_0t_1\cdots t_{u-1}
N=t0t1⋯tu−1，其中
u
=
∑
i
=
1
k
e
i
u=\displaystyle\sum_{i=1}^{k}e_i
u=i=1∑kei
对于任意整数
x
,
(
0
≤
x
≤
N
−
1
)
x,(0\leq x\leq N-1)
x,(0≤x≤N−1)，有式
(
2
)
(2)
(2)成立
x
=
c
0
+
c
1
t
0
+
c
2
t
0
t
1
+
⋯
+
c
u
−
1
t
0
t
1
⋯
t
u
−
2
(2)
x=c_0+c_1t_0+c_2t_0t_1+\cdots+c_{u-1}t_0t_1\cdots t_{u-2}\tag{2}
x=c0+c1t0+c2t0t1+⋯+cu−1t0t1⋯tu−2(2)
其中系数
c
i
c_i
ci满足
0
≤
c
i
<
t
i
−
1
0\leq c_i < t_{i-1}
0≤ci<ti−1
我们令
s
j
=
c
0
+
c
1
t
0
+
c
2
t
0
t
1
+
⋯
+
c
j
−
1
t
0
t
1
⋯
t
j
−
2
,
(
0
≤
j
<
u
−
1
)
s_j=c_0+c_1t_0+c_2t_0t_1+\cdots+c_{j-1}t_0t_1\cdots t_{j-2},(0\leq j < u-1)
sj=c0+c1t0+c2t0t1+⋯+cj−1t0t1⋯tj−2,(0≤j<u−1)
则
s
j
+
1
=
s
j
+
z
j
t
0
t
1
⋯
t
j
−
1
s_{j+1}=s_j+z_jt_0t_1\cdots t_{j-1}
sj+1=sj+zjt0t1⋯tj−1
由式
(
1
)
,
(
2
)
(1),(2)
(1),(2)得：
{
x
≡
c
0
(
m
o
d
t
0
)
x
−
c
0
≡
c
1
t
0
(
m
o
d
t
1
)
⋮
x
−
s
u
−
2
≡
c
u
−
1
t
0
t
1
⋯
t
u
−
2
(
m
o
d
t
u
−
1
)
(3)
\begin{cases}x\equiv c_0\ (mod\ \ t_0)\\ x-c_0\equiv c_1t_0\ (mod\ \ t_1)\\ \vdots\\ x-s_{u-2}\equiv c_{u-1}t_0t_1\cdots t_{u-2}\ (mod\ \ t_{u-1})\\ \end{cases}\tag{3}
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧x≡c0 (mod t0)x−c0≡c1t0 (mod t1)⋮x−su−2≡cu−1t0t1⋯tu−2 (mod tu−1)(3)
式
(
3
)
(3)
(3)即为求解的迭代函数
因此我们可以通过分解后的模数
p
i
e
i
p_i^{e_i}
piei设
x
i
=
c
0
+
c
1
p
i
+
c
2
p
i
2
+
⋯
+
c
e
i
−
1
p
i
e
i
−
1
(4)
x_i=c_0+c_1p_i+c_2p_i^2+\cdots+c_{e_i-1}p_i^{e_i-1}\tag{4}
xi=c0+c1pi+c2pi2+⋯+cei−1piei−1(4)
所以由
x
=
x
i
(
m
o
d
p
i
e
i
)
x=x_i\ (mod\ \ p_i^{e_i})
x=xi (mod piei)得
x
=
x
i
+
s
⋅
p
i
e
i
(5)
x=x_i+s·p_i^{e_i}\tag{5}
x=xi+s⋅piei(5)
由式
(
3
)
(
4
)
(3)(4)
(3)(4)知我们此次要求的
{
x
i
≡
c
0
(
m
o
d
p
i
)
x
i
−
c
0
≡
c
1
p
i
(
m
o
d
p
i
)
⋮
x
i
−
s
e
i
−
2
≡
c
e
i
−
1
p
i
e
i
−
1
(
m
o
d
p
i
)
(6)
\begin{cases}x_i\equiv c_0\ (mod\ \ p_i)\\ x_i-c_0\equiv c_1p_i\ (mod\ \ p_i)\\ \vdots\\ x_i-s_{e_i-2}\equiv c_{e_i-1}p_i^{e_i-1}\ (mod\ \ p_i)\\ \end{cases}\tag{6}
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧xi≡c0 (mod pi)xi−c0≡c1pi (mod pi)⋮xi−sei−2≡cei−1piei−1 (mod pi)(6)
由式
(
5
)
(
6
)
(5)(6)
(5)(6)求出
x
≡
x
i
(
m
o
d
p
i
e
i
)
x\equiv x_i\ (mod\ p_i^{e_i})
x≡xi (mod piei)
其中
c
0
c_0
c0的求法：
b
n
/
p
i
≡
(
a
x
)
n
/
p
i
(
m
o
d
p
)
≡
(
a
c
0
+
c
1
p
i
+
c
2
p
i
2
+
⋯
+
c
e
i
−
1
p
i
e
i
−
1
+
s
⋅
p
i
e
i
)
n
/
p
i
(
m
o
d
p
)
≡
(
a
c
0
)
n
/
p
i
(
m
o
d
p
)
(7)
\begin{aligned} b^{n/p_i}&\equiv (a^x)^{n/p_i}\ (mod\ \ p)\\ &\equiv (a^{c_0+c_1p_i+c_2p_i^2+\cdots+c_{e_i-1}p_i^{e_i-1}+s·p_i^{e_i}})^{n/p_i}\ (mod\ \ p)\\ &\equiv (a^{c_0})^{n/p_i}\ (mod\ \ p)\\ \end{aligned}\tag{7}
bn/pi≡(ax)n/pi (mod p)≡(ac0+c1pi+c2pi2+⋯+cei−1piei−1+s⋅piei)n/pi (mod p)≡(ac0)n/pi (mod p)(7)
学到这里的时候我觉得有点懵，实在不知道怎么把
c
1
p
i
+
c
2
p
i
2
+
⋯
+
c
e
i
−
1
p
i
e
i
−
1
+
s
⋅
p
i
e
i
c_1p_i+c_2p_i^2+\cdots+c_{e_i-1}p_i^{e_i-1}+s·p_i^{e_i}
c1pi+c2pi2+⋯+cei−1piei−1+s⋅piei舍掉的。
没找到解释，我是通过原根推的，我推的是（不知道正确性？）：
由原根性质，设
g
g
g是
p
p
p的一个原根
令
b
=
g
k
,
a
=
g
t
b=g^k,a=g^t
b=gk,a=gt，所以有
g
k
≡
(
g
t
)
x
(
m
o
d
p
)
g^k\equiv (g^t)^x\ (mod\ \ p)
gk≡(gt)x (mod p),其中
t
x
≡
k
(
m
o
d
p
−
1
)
tx\equiv k\ (mod\ \ p-1)
tx≡k (mod p−1)
(
g
k
)
p
−
1
p
i
≡
[
(
g
t
)
c
0
+
c
1
p
i
+
c
2
p
i
2
+
⋯
+
c
e
i
−
1
p
i
e
i
−
1
+
s
⋅
p
i
e
i
]
p
−
1
p
i
(
m
o
d
p
)
g
k
⋅
(
p
−
1
)
p
i
≡
g
t
⋅
(
p
−
1
)
⋅
c
0
p
i
+
t
⋅
(
p
−
1
)
⋅
c
1
+
t
⋅
(
p
−
1
)
c
2
p
i
+
⋯
+
t
⋅
(
p
−
1
)
⋅
s
⋅
p
i
e
i
−
1
(
m
o
d
p
)
\begin{aligned} (g^k)^{\frac{p-1}{p_i}}&\equiv [(g^t)^{c_0+c_1p_i+c_2p_i^2+\cdots+c_{e_i-1}p_i^{e_i-1}+s·p_i^{e_i}}]^{\frac{p-1}{p_i}}\ (mod\ \ p)\\ g^{\frac{k·(p-1)}{p_i}}&\equiv g^{\frac{t·(p-1)·c_0}{p_i}+t·(p-1)·c_1+t·(p-1)c_2p_i+\cdots+t·(p-1)·s·p_i^{e_i-1}}\ (mod\ \ p) \end{aligned}
(gk)pip−1gpik⋅(p−1)≡[(gt)c0+c1pi+c2pi2+⋯+cei−1piei−1+s⋅piei]pip−1 (mod p)≡gpit⋅(p−1)⋅c0+t⋅(p−1)⋅c1+t⋅(p−1)c2pi+⋯+t⋅(p−1)⋅s⋅piei−1 (mod p)
因此
k
⋅
(
p
−
1
)
p
i
≡
t
⋅
(
p
−
1
)
⋅
c
0
p
i
+
t
⋅
(
p
−
1
)
⋅
c
1
+
t
⋅
(
p
−
1
)
c
2
p
i
+
⋯
+
t
⋅
(
p
−
1
)
⋅
s
⋅
p
i
e
i
−
1
(
m
o
d
p
−
1
)
\frac{k·(p-1)}{p_i}\equiv \frac{t·(p-1)·c_0}{p_i}+t·(p-1)·c_1+t·(p-1)c_2p_i+\cdots+t·(p-1)·s·p_i^{e_i-1}\ (mod\ \ p-1)
pik⋅(p−1)≡pit⋅(p−1)⋅c0+t⋅(p−1)⋅c1+t⋅(p−1)c2pi+⋯+t⋅(p−1)⋅s⋅piei−1 (mod p−1)
根据模运算显然
k
⋅
(
p
−
1
)
p
i
≡
t
⋅
(
p
−
1
)
⋅
c
0
p
i
(
m
o
d
p
−
1
)
\frac{k·(p-1)}{p_i}\equiv \frac{t·(p-1)·c_0}{p_i}\ (mod\ \ p-1)
pik⋅(p−1)≡pit⋅(p−1)⋅c0 (mod p−1)
显然式
(
7
)
(7)
(7)的阶是
n
n
n，原根的阶是
p
−
1
p-1
p−1，性质相同，故式
(
7
)
(7)
(7)成立。
因此我们可以通过
b
n
/
p
i
≡
(
a
c
0
)
n
/
p
i
(
m
o
d
p
)
b^{n/p_i}\equiv (a^{c_0})^{n/p_i}\ (mod\ \ p)
bn/pi≡(ac0)n/pi (mod p)求解出
c
0
c_0
c0
因为
c
0
∈
[
0
,
p
i
−
1
]
c_0\in [0,p_i-1]
c0∈[0,pi−1]，所以
p
i
p_i
pi不是很大的情况下完全可以解。
依旧是BSGS的思想，我们可以将
a
n
/
p
i
,
a
2
n
/
p
i
,
⋯
,
a
(
p
i
−
1
)
n
/
p
i
a^{n/p_i},a^{2n/p_i},\cdots,a^{(p_i-1)n/p_i}
an/pi,a2n/pi,⋯,a(pi−1)n/pi存储。
然后查找即可。
我们令
b
1
=
b
⋅
a
−
c
0
b_1=b·a^{-c_0}
b1=b⋅a−c0,所以可以通过
b
1
n
/
p
i
2
≡
(
a
c
1
)
n
/
p
i
(
m
o
d
p
)
b_1^{n/p_i^2}\equiv (a^{c_1})^{n/p_i}\ (mod\ \ p)
b1n/pi2≡(ac1)n/pi (mod p)求解出
c
1
c_1
c1
我们令
b
2
=
b
1
⋅
a
−
c
1
p
i
b_2=b_1·a^{-c_1p_i}
b2=b1⋅a−c1pi,所以可以通过
b
2
n
/
p
i
3
≡
(
a
c
2
)
n
/
p
i
(
m
o
d
p
)
b_2^{n/p_i^3}\equiv (a^{c_2})^{n/p_i}\ (mod\ \ p)
b2n/pi3≡(ac2)n/pi (mod p)求解出
c
2
c_2
c2
其他同理。
这样我们就已知了所有同余式
x
≡
x
i
(
m
o
d
p
i
e
i
)
x\equiv x_i (mod\ \ p_i^{e_i})
x≡xi(mod piei)
当然，如果不存在这样的同余式，也就代表着无解。
然后一发CRT即可。
代码
#include <bits/stdc++.h>
#include <iostream>
#include <algorithm>
#include <cstdio>
#include <queue>
#include <cmath>
#include <string>
#include <cstring>
#include <map>
#include <set>
#include <vector>
#include <cmath>
#include <ctime>
using namespace std;
#define me(x,y) memset(x,y,sizeof x)
#define MIN(x,y) (x) < (y) ? (x) : (y)
#define MAX(x,y) (x) > (y) ? (x) : (y)
#define SGN(x) ((x)>0?1:((x)<0?-1:0))
#define ABS(x) ((x)>0?(x):-(x))
typedef long long ll;
typedef unsigned long long ull;
const int maxn = 1e5+10;
const ll INF = 0x3f3f3f3f;
const int MOD = 1e9+7;
const int eps = 1e-8;
ll qmul(ll a,ll b,ll p){
return (a*b-(ll)((long double)a*b/p)*p+p)%p;
}
ll qpow(ll a,ll b,ll p){
ll ans=1;
while(b){
if(b&1) ans = qmul(ans,a,p);
a = qmul(a,a,p);
b >>= 1;
}
return ans;
}
const int S=5;
//以a为基,n-1=x*2^t a^(n-1)=1(mod n) 验证n是不是合数
//一定是合数返回true,不一定返回false
bool check(long long a,long long n,long long x,long long t)
{
long long ret=qpow(a,x,n);
long long last=ret;
for(int i=1;i<=t;i++)
{
ret=qmul(ret,ret,n);
if(ret==1&&last!=1&&last!=n-1) return true;//合数
last=ret;
}
if(ret!=1) return true;
return false;
}
// Miller_Rabin()算法素数判定
//是素数返回true.(可能是伪素数，但概率极小)
//合数返回false;
bool Miller_Rabin(long long n)
{
if(n<2)return false;
if(n==2)return true;
if((n&1)==0) return false;//偶数
long long x=n-1;
long long t=0;
while((x&1)==0){x>>=1;t++;}
for(int i=0;i<S;i++)
{
long long a=rand()%(n-1)+1;
if(check(a,n,x,t))
return false;//合数
}
return true;
}
long long factor[100];//质因数分解结果（刚返回时是无序的）
int tol;//质因数的个数。数组小标从0开始
long long gcd(long long a,long long b)
{
if(a==0)return 1;
if(a<0) return gcd(-a,b);
while(b)
{
long long t=a%b;
a=b;
b=t;
}
return a;
}
long long Pollard_rho(long long x,long long c)
{
long long i=1,k=2;
long long x0=rand()%x;
long long y=x0;
while(1)
{
i++;
x0=(qmul(x0,x0,x)+c)%x;
long long d=gcd(y-x0,x);
if(d!=1&&d!=x) return d;
if(y==x0) return x;
if(i==k){y=x0;k+=k;}
}
}
//对n进行素因子分解
map<ll,ll> mp;
void findfac(long long n)
{
if(Miller_Rabin(n))//素数
{
factor[tol++]=n;
mp[n]++;
return;
}
long long p=n;
while(p>=n) p=Pollard_rho(p,rand()%(n-1)+1);
findfac(p);
findfac(n/p);
}
ll get_n(ll x,ll a,ll p){
tol = 0;
mp.clear();
findfac(x);
for(int i= 0; i < tol; ++i){
while(x%factor[i] == 0 && qpow(a,x/factor[i],p) == 1)
x /= factor[i];
}
mp.clear();
findfac(x);
return x;
}
void ex_gcd(ll a,ll b,ll &d,ll &x,ll &y){
if(!b){
d = a;x = 1;y = 0;return;
}
ex_gcd(b,a%b,d,y,x);
y -= x*(a/b);
}
ll crt(int n,ll *a,ll *m){ // x%m[i]=a[i]
ll M = 1,ret = 0;
ll x,y,d;
for(int i = 1; i <= n; ++i) M *= m[i];
for(int i = 1; i <= n; ++i){
ll Mi = M / m[i];
ex_gcd(Mi,m[i],d,x,y);
x = (x%m[i]+m[i])%m[i];
ret = (ret+qmul(qmul(a[i],Mi,M),x,M))%M;
}
return (ret+M)%M;
}
map<ll,ll> hs;
ll get_c(ll b,ll p,ll pi){
for(int i = 0; i < pi; ++i){
if(b == hs[i]) return i;
}
return -1;
}
ll xx[maxn],pp[maxn];
ll Pohlig_Hellman(ll a,ll b,ll p){
ll n = get_n(p-1,a,p); //a对p的阶
ll x,y,inv = qpow(a,p-2,p),col = 0;
for(auto it : mp){ //枚举素因子
hs.clear();
int pi = it.first,count = it.second;
ll tmp = 1,cnt = qpow(a,n/pi,p);
ll bb = b,an = 0;
for(int i = 0; i < pi; ++i){ //哈希存储a^(n/pi),a^(2n/pi),a^(3n/pi)...
hs[i] = tmp;
// cout<<"hs["<<i<<"]:"<<hs[i]<<endl;
tmp = qmul(cnt,tmp,p);
}
ll now = n;
for(int i = 1; i <= count; ++i){
now = n/qpow(pi,i,p); // n / (pi^i)
ll res = qpow(bb,now,p),c;
c = get_c(res,p,pi);
// cout<<"c: "<<c<<" ";
if(c == -1) return -1;
ll rel = c*qpow(pi,i-1,p);
an += rel;
ll inva = qpow(qpow(a,rel,p),p-2,p);
// cout<<"inva:"<<inva<<" ";
bb = qmul(bb,inva,p);
// cout<<"bb: "<<bb<<endl;
}
// cout<<endl;
xx[++col] = an;
pp[col] = qpow(pi,count,p);
// cout<<"xx: " <<xx[col]<<" pp:"<<pp[col]<<endl;
}
if(col == 0) return -1;
return crt(col,xx,pp);
}
int main(){
ios::sync_with_stdio(false);
int t;
cin>>t;
while(t--){
ll p,a,b;
cin>>p>>a>>b;
ll ans = Pohlig_Hellman(a,b,p);
if(ans == -1) cout<<-1<<endl;
else cout<<ans<<endl;
}
return 0;
}
2835
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
