写在前面
一切都要从去年8月的多校开始说起。
笔者在多校某场比赛中，遇到一个离散对数的题：HDU 6632
给定
a
,
b
,
p
(
p
∈
p
r
i
m
e
,
65537
≤
p
≤
1
e
18
,
2
≤
a
,
b
≤
p
−
1
a,b,p(p∈prime,65537\leq p \leq 1e18,2\leq a,b\leq p-1
a,b,p(p∈prime,65537≤p≤1e18,2≤a,b≤p−1且
p
−
1
p-1
p−1的质因子只能由
2
2
2(和/或)
3
3
3组成)，求
a
x
≡
b
(
m
o
d
p
)
a^x\equiv b(mod\ p)
ax≡b(mod p)的最小正整数解，若无解则输出
−
1
-1
−1。
因为质数
p
p
p很大，且通常使用的
B
S
G
S
BSGS
BSGS复杂度为
O
(
n
)
O(\sqrt n)
O(n)，显然无法通过正常方法解题。于是很理所应当地翻车了。
随后出题人给的题解中提到了pohlig-hellman算法，想去学，但国内对于该算法的相关资料实在是少之又少。于是在理应不存在的404网站们找到了相关视频学习后才醒悟（一醒就是半年，咕咕咕）。于是决定整理一下这个问题。
注：笔者是一名十八线蒟蒻ACMer，因此下列文章可能有许多细节上的问题，请您指正，我将尽快回复和改正，非常感谢！
前置技能
快速幂
基础数论知识
扩展gcd
欧拉函数
原根
阶
中国剩余定理（CRT）
质因子分解
pohlig-hellman算法
需要注意的是，pohlig-hellman算法的复杂度在一般情况下比
B
S
G
S
BSGS
BSGS高！
因此，使用pohlig-hellman的场合只能是较为特殊的情况，即:
1.p是质数，且p-1包含的质因子较少&较小。
和
B
S
G
S
BSGS
BSGS算法一样，pohlig-hellman算法也是用于解决离散对数问题（也有很多文献提到是解决椭圆曲线之类的）。
即给定a,b,p,求
a
x
≡
b
(
m
o
d
p
)
a^x \equiv b(mod\ p)
ax≡b(mod p)。
在开始讲算法前，首先要提一下欧拉定理。
若 ( a , p ) = 1 , 那么 a φ ( p ) ≡ 1 ( m o d p ) ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ( ∗ ) (a,p)=1,那么a^{φ(p)} \equiv1(mod\ p)······(*) (a,p)=1,那么aφ(p)≡1(mod p)⋅⋅⋅⋅⋅⋅(∗)
证明略。
如果
p
p
p是质数，那
φ
(
p
)
=
p
−
1
φ(p)=p-1
φ(p)=p−1，即费马小定理。
对于
a
x
≡
b
(
m
o
d
p
)
a^x\equiv b(mod\ p)
ax≡b(mod p)，记其原根为
g
g
g。
则
a
=
g
a
i
,
b
=
g
b
i
a=g^{a_i},b=g^{b_i}
a=gai,b=gbi(原根的次幂可以在[1,p-1]中一一对应，故
a
i
,
b
i
a_i,b_i
ai,bi必定存在)
即
g
a
i
x
≡
g
b
i
(
m
o
d
p
)
g^{a_ix}\equiv g^{b_i}(mod\ p)
gaix≡gbi(mod p)
所以
a
i
x
≡
b
i
(
m
o
d
p
−
1
)
a_ix\equiv b_i(mod\ p-1)
aix≡bi(mod p−1)
(因为
g
p
−
1
≡
1
(
m
o
d
p
)
g^{p-1}\equiv 1(mod\ p)
gp−1≡1(mod p),
p
−
1
p-1
p−1是最小循环节，即阶)
故利用扩展gcd求出
x
x
x即可。
于是我们的问题就变成了如何求
a
i
,
b
i
a_i,b_i
ai,bi。
接下来，开始将pohlig-hellman算法的步骤（此处
a
a
a是原根，
p
p
p是质数）。
以下步骤建议读者自行手写一遍！
以下步骤建议读者自行手写一遍！
以下步骤建议读者自行手写一遍！
重要的话说三遍！
1.将
φ
φ
φ
(
p
)
(p)
(p)=p-1质因子分解
即
p
−
1
=
p
1
k
1
p
2
k
2
.
.
.
p
m
k
m
(
m
为质因子个数
)
即p-1=p^{k_1}_1p^{k_2}_2...p^{k_m}_m(m为质因子个数)
即p−1=p1k1p2k2...pmkm(m为质因子个数)
2.对于每个质因子
p
i
p_i
pi，列出方程
x
=
p
i
0
a
0
+
p
i
1
a
1
+
p
i
2
a
2
.
.
.
+
p
i
k
i
−
1
a
k
i
−
1
(
m
o
d
p
i
k
i
)
x=p^0_ia_0+p^1_ia_1+p^2_ia_2...+p^{k_{i-1}}_ia_{k_{i-1}}(mod\ p_i^{k_i})
x=pi0a0+pi1a1+pi2a2...+piki−1aki−1(mod piki)
(即把x写成
p
i
p_i
pi进制，每个系数都小于
p
i
p_i
pi)。
3.
令
r
=
1
,
求
(
a
x
)
p
−
1
p
i
r
≡
b
p
−
1
p
i
r
(
m
o
d
p
)
令r=1,求(a^x)^{\frac{p-1}{p_i^r}}\equiv b^{\frac{p-1}{p_i^r}}(mod\ p)
令r=1,求(ax)pirp−1≡bpirp−1(mod p)
（注意这个
p
−
1
p
i
r
{\frac{p-1}{p_i^r}}
pirp−1，这是我们求出
a
i
a_i
ai各项系数的关键！）
4.展开 x x x,得到 ( a p i 0 a 0 + p i 1 a 1 + p i 2 a 2 . . . + p i k i − 1 a k i − 1 ) p − 1 p i r ≡ b p − 1 p i r ( m o d p ) (a^{p^0_ia_0+p^1_ia_1+p^2_ia_2...+p^{k_{i-1}}_ia_{k_{i-1}}})^{\frac{p-1}{p_i^r}}\equiv b^{\frac{p-1}{p_i^r}}(mod\ p) (api0a0+pi1a1+pi2a2...+piki−1aki−1)pirp−1≡bpirp−1(mod p)
5.展开，得(以下以r=1为例)
a
a
0
∗
p
−
1
p
i
∗
a
(
p
−
1
)
a
1
∗
(
a
p
−
1
)
p
a
2
∗
(
a
p
−
1
)
p
2
a
3
.
.
.
(
a
p
−
1
)
p
(
k
i
−
2
)
a
k
i
−
1
≡
b
p
−
1
p
i
r
(
m
o
d
p
)
a^{a_0*\frac{p-1}{p_i}}*a^{(p-1)a_1}*(a^{p-1})^{pa_2}*(a^{p-1})^{p^2a_3}...(a^{p-1})^{p^{(k_{i}-2)}a_{k_{i-1}}}\equiv b^{\frac{p-1}{p_i^r}}(mod\ p)
aa0∗pip−1∗a(p−1)a1∗(ap−1)pa2∗(ap−1)p2a3...(ap−1)p(ki−2)aki−1≡bpirp−1(mod p)
6.注意到从第二项开始，后面每一项都是1(因为由欧拉定理得 a p − 1 ≡ 1 ( m o d p ) a^{p-1}\equiv1(mod\ p) ap−1≡1(mod p)，1的任何次方都是1)，故我们可以得到式子
a
a
0
∗
p
−
1
p
i
≡
b
p
−
1
p
i
r
(
m
o
d
p
)
a^{a_0*\frac{p-1}{p_i}} \equiv b^{\frac{p-1}{p_i^r}}(mod\ p)
aa0∗pip−1≡bpirp−1(mod p)
因为a_i的所有系数∈
[
0
,
p
i
−
1
]
[0,p_i-1]
[0,pi−1](前文提到了
p
i
p_i
pi进制)，故可以在
O
(
p
i
)
O(p_i)
O(pi)的时间内暴力求出
a
0
a_0
a0。
7.得到
a
0
a_0
a0后，令
r
=
r
+
1
r=r+1
r=r+1,回到步骤3。
直到所有a_i被求出后，可得到：
x
=
p
i
0
a
0
+
p
i
1
a
1
+
p
i
2
a
2
.
.
.
+
p
i
k
i
−
1
a
k
i
−
1
(
m
o
d
p
i
k
i
)
。
x=p^0_ia_0+p^1_ia_1+p^2_ia_2...+p^{k_{i-1}}_ia_{k_{i-1}}(mod\ p_i^{k_i})。
x=pi0a0+pi1a1+pi2a2...+piki−1aki−1(mod piki)。
(
注意模数不是
p
！
)
(注意模数不是p！)
(注意模数不是p！)
8.在算出m个关于x的式子后,来一发CRT（中国剩余定理）即可(不需要扩展CRT,因为模数互质)。
需要注意！CRT方程组必定有解，因此，如果原根变换后的最后的扩展gcd无解，也说明原方程无解！！！
举例
下面，我们通过一个例子来帮助理解以上的步骤。
eg1:
7
x
≡
12
(
m
o
d
41
)
7^x\equiv12(mod\ 41)
7x≡12(mod 41)
解：首先
41
41
41的最小原根是
6
6
6。
故令
7
=
6
x
1
,
12
=
6
x
2
7=6^{x_1},12=6^{x_2}
7=6x1,12=6x2
先计算
6
x
≡
7
(
m
o
d
41
)
6^{x}\equiv 7(mod\ 41)
6x≡7(mod 41)
分解
φ
(
41
)
=
40
=
2
3
5
。
φ(41)=40=2^35。
φ(41)=40=235。
对于质因子2，我们首先列出方程
x
=
2
0
a
0
+
2
1
a
1
+
2
2
a
2
(
a
i
∈
[
0
,
1
]
)
x=2^0a_0+2^1a_1+2^2a_2(a_i∈[0,1])
x=20a0+21a1+22a2(ai∈[0,1])
在等式两边同乘以
p
−
1
2
=
20
\frac{p-1}{2}=20
2p−1=20,得
6
x
20
≡
7
20
(
m
o
d
41
)
。
6^{x^{20}}\equiv 7^{20}(mod\ 41)。
6x20≡720(mod 41)。
将 x x x展开，得
6 ( 2 0 a 0 + 2 1 a 1 + 2 2 a 2 ) 20 ≡ 7 20 ( m o d 41 ) 6^{(2^0a_0+2^1a_1+2^2a_2)^{20}}\equiv 7^{20}(mod\ 41) 6(20a0+21a1+22a2)20≡720(mod 41)
化简，得
6 20 a 0 ∗ 6 40 a 1 ∗ 6 8 0 2 a ∗ ≡ 7 20 ( m o d 41 ) 6^{20a_0}*6^{40a_1}*6^{80^a_2}*\equiv 7^{20}(mod\ 41) 620a0∗640a1∗6802a∗≡720(mod 41)
由费马小定理/欧拉定理得后面两项都是1，故：
6 20 a 0 ≡ 7 20 ( m o d 41 ) 6^{20a_0} \equiv 7^{20}(mod\ 41) 620a0≡720(mod 41)
快速幂优化，得
4 0 a 0 ≡ 40 ( m o d 41 ) 40^{a_{0}}\equiv 40(mod\ 41) 40a0≡40(mod 41)
在 [ 0 , 1 ] [0,1] [0,1]内枚举 a 0 a_0 a0,得到 a 0 = 1 。 a_0=1。 a0=1。
有了 a 1 a_1 a1后，我们将原式同乘以 p − 1 2 2 = 10 \frac{p-1}{2^2}=10 22p−1=10,得到 6 x 10 ≡ 7 10 ( m o d 41 ) 。 6^{x^{10}}\equiv 7^{10}(mod\ 41)。 6x10≡710(mod 41)。
重复上述步骤，得
6 ( 1 ∗ 2 1 + a 1 + 2 2 a 2 ) 10 ≡ 7 10 ( m o d 41 ) 6^{(1*2^1+a_1+2^2a_2)^{10}}\equiv 7^{10}(mod\ 41) 6(1∗21+a1+22a2)10≡710(mod 41)
化简，得
6 10 ∗ 6 20 a 1 ∗ 6 40 a 2 ≡ 7 10 ( m o d 41 ) 6^{10}*6^{20a_1}*6^{40a_2}\equiv 7^{10}(mod\ 41) 610∗620a1∗640a2≡710(mod 41)
同样道理，后面一项是1，得
6 10 ∗ 6 2 0 a 1 ≡ 7 10 ( m o d 41 ) 6^{10}*6^{20^{a_1}}\equiv 7^{10}(mod\ 41) 610∗620a1≡710(mod 41)
快速幂优化，得
32 ∗ 4 0 a 1 ≡ 9 ( m o d 41 ) 32*40^{a_1}\equiv9(mod\ 41) 32∗40a1≡9(mod 41)
在 [ 0 , 1 ] [0,1] [0,1]内枚举 a 1 a_1 a1，得 a 1 = 1 。 a_1=1。 a1=1。
同理我们得到 a 2 = 1 a_2=1 a2=1
于是 x ≡ 2 0 ∗ 1 + 2 1 ∗ 1 + 2 2 ∗ 1 ≡ 7 ( m o d 8 ) x\equiv2^0*1+2^1*1+2^2*1\equiv7(mod\ 8) x≡20∗1+21∗1+22∗1≡7(mod 8)
即 x ≡ 7 ( m o d 8 ) 即x\equiv7(mod\ 8) 即x≡7(mod 8)
对于质因子
5
5
5，同理，设
x
=
5
0
a
0
(
a
i
∈
[
0
,
4
]
)
x=5^0a_0(a_i∈[0,4])
x=50a0(ai∈[0,4])，得到
a
0
=
4
a_0=4
a0=4，
即
x
=
4
(
m
o
d
5
)
。
x=4(mod\ 5)。
x=4(mod 5)。
综上所述，我们得到：
{
x
≡
7
(
m
o
d
8
)
x
≡
4
(
m
o
d
5
)
\begin{cases} x\equiv 7(mod\ 8)\\ x\equiv 4(mod\ 5)\\ \end{cases}
{x≡7(mod 8)x≡4(mod 5)
由中国剩余定理可以得到
x
≡
39
(
m
o
d
40
)
。
x\equiv39(mod\ 40)。
x≡39(mod 40)。
即
6
39
≡
7
(
m
o
d
41
)
6^{39}\equiv 7(mod\ 41)
639≡7(mod 41)
同理可以求出
6
27
≡
12
(
m
o
d
41
)
6^{27}\equiv12(mod\ 41)
627≡12(mod 41)
因此
6
3
9
x
≡
6
27
(
m
o
d
41
)
6^{39^x}\equiv 6^{27}(mod\ 41)
639x≡627(mod 41)
39
x
≡
27
(
m
o
d
40
)
39x\equiv27(mod\ 40)
39x≡27(mod 40)
解得
x
≡
13
(
m
o
d
40
)
x\equiv13(mod\ 40)
x≡13(mod 40)
故
7
x
≡
12
(
m
o
d
41
)
7^x\equiv 12(mod\ 41)
7x≡12(mod 41)的解是
x
≡
13
(
m
o
d
40
)
x\equiv13(mod\ 40)
x≡13(mod 40)
代码
需要用快速乘优化快速幂，对于模数是1e18左右时还要考虑到先除后乘等防溢出的操作。
UPD：代码内的质数因子是提前筛出来的，可以用质数筛法存入prime数组内或提前筛出写入数组。
#include<bits/stdc++.h>
using namespace std;
const long long numprime=2;//因子个数
long long prime[numprime+1]={0,2,3};
//质因子，需要提前筛出，此处只写了两个，请根据需要进行添加或跑一遍任意筛法存储。
long long factor[numprime+1]={0};//因子次数
long long xishu[200]={0};//系数
long long pow_prime[200];//存储pi^k
long long numc=0LL;//crt方程总数
long long crta[numprime+1],crtb[numprime+1];//a是余数，b是模数
long long quick_mul (long long a,long long b,long long c)//快速乘
{
return (a*b-(long long)((long double)a*b/c)*c+c)%c;
}
long long quick_pow (long long a,long long b,long long c)//快速幂
{
long long ans=1,base=a;
while (b!=0)
{
if (b&1)
ans=quick_mul (ans,base,c);
base=quick_mul (base,base,c);
b>>=1;
}
return ans%c;
}
void fenjie(long long p)//质因子分解
{
long long p0=p-1;
for (int i=1;i<=numprime;i++)
{
factor[i]=0;
while (p0%prime[i]==0)
{
factor[i]++;
p0/=prime[i];
}
}
}
bool check(long long a,long long p)//判断是否是原根
{
for(int i=1;i<=numprime;i++)
{
if(quick_pow (a,(p-1)/prime[i],p)==1)
return 0;
}
return 1;
}
long long Root (long long p)//遍历寻找原根
{
for (long long g=2;g<p;g++)
{
if(check(g,p))
return g;
}
}
long long extended_gcd (long long a,long long b,long long &x,long long &y)
//扩展gcd
{
long long r,t;
if (b==0)
{
x=1;
y=0;
return a;
}
r=extended_gcd (b,a%b,x,y);
t=x;
x=y;
y=t-(a/b)*y;
return r;
}
long long Inv (long long a,long long mod)//逆元
{
long long p,x,y;
p=extended_gcd (a,mod,x,y);//p是gcd=1
long long t=mod;
if (x>=0)
x=x%t;
else
{
while (x<0)
x+=t;
x%t;
}
if (x==0)
x+=t;
return x;
}
long long CRT (long long numc,long long a[],long long b[],long long p)
//中国剩余定理
{
long long ans=0;
long long mod=p-1;
for (long long i=1;i<=numc;i++)
{
long long x1=quick_mul (a[i]%mod,mod/b[i]%mod,mod);
x1=quick_mul (x1%mod,Inv (mod/b[i],b[i])%mod,mod);
ans=(ans+x1+mod)%mod;
}
return ans;
}
long long pohlig_hellman (long long g,long long b,long long p)
{
numc=0;
//g^x=b(mod p)
memset (crta,0,sizeof (crta));
memset (crtb,0,sizeof (crtb));
for (int i=1;i<=numprime;i++)//枚举每个质因子
{
memset (xishu,0,sizeof (xishu));
pow_prime[0]=1;
for (int j=1;j<=factor[i]+1;j++)
{
pow_prime[j]=pow_prime[j-1]*prime[i];//预处理
}
long long sum=1;
long long yushu=0,moshu=pow_prime[factor[i]];
for (int j=1;j<=factor[i];j++)//求出每个系数
{
long long a0=quick_pow (g,(p-1)/prime[i],p);
long long b0=quick_pow (b,(p-1)/pow_prime[j],p);
for (long long x=0;x<=prime[i]-1;x++)//遍历找出系数x
{
if (quick_mul (sum,quick_pow (a0,x,p),p)==b0)
{
xishu[j]=x;
long long xs=0;//已求出的系数按权求和
for (long long k=1;k<=j;k++)
{
xs=(xs+(pow_prime[k-1]*xishu[k]%(p-1))+p-1)%(p-1)+p-1;
}
sum=quick_pow (g,quick_mul (xs,(p-1)/pow_prime[j+1],p-1)+p-1,p);
break;
}
}
yushu+=pow_prime[j-1]*xishu[j];
}
if (factor[i])
{
numc++;
crta[numc]=yushu;
crtb[numc]=moshu;
}
}
return CRT (numc,crta,crtb,p);
}
int main ()
{
int t;
scanf ("%d",&t);
while (t--)
{
long long a,b,p,g,xx,yy,rr,x1,x2;
scanf ("%lld%lld%lld",&p,&a,&b);
fenjie (p);
if(check(a,p))// 如果a是原根，那只要做一遍pohlig-hellman
{
printf("%lld\n",pohlig_hellman (a,b,p));
continue;
}
g=Root (p);
// printf ("%lld's Primitive Root :%lld\n",p,g);
x1=pohlig_hellman (g,a,p);
x2=pohlig_hellman (g,b,p);
// printf("solve %lldx=%lld(mod %lld)\n",x1,x2,p-1);//x1x=x2(mod p-1)
rr=extended_gcd (x1,p-1,xx,yy);
// printf("gcd=%lld\n",rr);
if (x2%rr!=0)
{
printf ("-1\n");
continue;
}
long long tt=(p-1)/rr;
if (xx>=0)
{
xx=quick_mul (xx,(x2/rr),tt);
}
else
{
xx=xx%tt+tt;
xx=quick_mul (xx,(x2/rr),tt);
}
printf ("%lld\n",xx);
}
}
后记
写这篇文章的动机很简单：我没看懂国内关于pohlig-hellman算法的文章（肯定是我语文差），因此在某些网站上学懂之后想着自己写一篇文章，一来是方便日后复习，二来帮助下和我语文一样差的程序猿…
顺便一提，模数是有原根的合数时，可能也适用，具体还需要进一步确认，先咕咕咕。
UPD:感谢群聚njust-poorwill的指点。
UPD2:添加了部分描述，避免产生理解上的问题
DrGilbert 2020.1.21
1731
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
