discrete logarithm problem
Problem Description
You are given three integers p,a,b, where p is a prime number and p−1 only has prime factors 2 and/or 3. Please find the minimum positive integer x such that ax≡b(modp).
Input
The first line contains an integer T indicating there are T tests. Each test consists of a single line containing three integers: p,a,b.
-
T≤200
-
65537 ≤ p ≤ 1 0 18 65537≤p≤10^{18} 65537≤p≤1018
-
the prime factors of p−1 can only be 2 or 3
-
2≤a,b≤p−1
Output
For each test, output a line containing an integer x, representing the minimum positive value such that ax≡b(modp). If there didn’t exist any such number x, please output −1.
Sample Input
6
65537 2 3
65537 2 4
65537 3 4
65537 4 4
65537 5 4
666334875701477377 2 3
Sample Output
-1
2
45056
1
36864
1957714645490451
题意
已知
p
,
a
,
b
p,a,b
p,a,b,求解
a
x
≡
b
(
m
o
d
p
)
a^x\equiv b\ (mod\ \ p)
ax≡b (mod p)
题解
因为
p
≤
1
0
18
p\leq 10^{18}
p≤1018,而离散对数求解算法BSGS的复杂度是
p
\sqrt{p}
p。
赛后题解说是CTF的简单算法, Pohlig Hellman。但是没学过密码学,等着看看再来补一发 Pohlig Hellman的。
补完了补完了,代码在下面,Pohlig Hellman讲解在这里,水平有限,不一定讲得清楚。
时间直接从1000+ms降到62ms,爽啊
回去重学了原根,然后反复研究才弄明白大佬不明不白的三言两语的题解。
该题做法是:
<1> 先求 g g g, g g g为 p p p的原根。原根是判断所有 g p − 1 p i % p ! = 1 g^{\frac{p-1}{p_i}}\%p!=1 gpip−1%p!=1,题中说了 p − 1 p-1 p−1素因子只有2,3。所以判断条件就成了 g p − 1 2 % p ! = 1 g^{\frac{p-1}{2}}\%p!=1 g2p−1%p!=1, g p − 1 3 % p ! = 1 g^{\frac{p-1}{3}}\%p!=1 g3p−1%p!=1
<2> 求解 p − 1 = 2 k 3 t p-1=2^k3^t p−1=2k3t中的 k , t k,t k,t(后面用)
<3> 求解
a
=
g
z
,
b
=
g
c
a=g^z,b=g^c
a=gz,b=gc中的
z
,
c
z,c
z,c。根据原根性质可以知道对于质数
p
p
p,其
φ
(
p
)
=
p
−
1
\varphi(p)=p-1
φ(p)=p−1,所以我们可以将
{
1
,
2
,
⋯
 
,
p
−
1
}
\{1,2,\cdots,p-1\}
{1,2,⋯,p−1}和
{
g
1
,
g
2
,
⋯
 
,
g
p
−
1
}
\{g^1,g^2,\cdots,g^{p-1}\}
{g1,g2,⋯,gp−1}建立一一对应关系。也就知道存在唯一的
z
,
c
,
(
z
,
c
≤
p
−
1
)
z,c,(z,c\leq p-1)
z,c,(z,c≤p−1)使得等式成立。
但是因为
p
≤
1
0
18
p\leq 10^{18}
p≤1018,无法暴力寻找
z
,
c
z,c
z,c,也无法用离散对数求解。需要其他方法(下面)
<4> 因此式子由 a x = b ( m o d p ) a^x=b\ (mod\ \ p) ax=b (mod p) 转换为 g z x = g c ( m o d p ) g^{zx}=g^c\ (mod\ \ p) gzx=gc (mod p)。因为<3>我们又可以将其转化为 z x = c ( m o d p − 1 ) zx=c\ (mod\ \ p-1) zx=c (mod p−1),这是一个单变元模线性方程,可以用exgcd求解。因为模数为素数所以方程为 z x + ( p − 1 ) y = g c d ( z , p − 1 ) zx+(p-1)y=gcd(z,p-1) zx+(p−1)y=gcd(z,p−1),无解输出-1,有解即可求出最小正整数解x
大佬原话说(没错,大佬口头说的题解就这么长):
对于一个数x,先判断其有没有模p平方根,没有的话就乘以原根g,再判断有没有四次方根,没有的话乘以
g
2
g^2
g2,再判断有没有八次方根,没有的话就乘以
g
4
g^4
g4。不断循环,就能找到x乘以g的几次方等于1,也就由
g
p
−
1
=
1
g^{p-1}=1
gp−1=1知道了x等于g的几次方。
我的理解是:(表达能力实在有限。。。)
本题因为p-1只有素因子2,3,
g
p
−
1
%
p
=
g
2
k
3
t
%
p
=
1
g^{p-1}\%p=g^{2^k3^t}\%p=1
gp−1%p=g2k3t%p=1,所以只需要判断
g
1
,
g
2
,
g
2
2
,
⋯
 
,
g
2
k
3
,
g
2
k
3
2
,
⋯
 
,
g
2
k
3
t
g^1,g^2,g^{2^2},\cdots,g^{2^k3},g^{2^k3^2},\cdots,g^{2^k3^t}
g1,g2,g22,⋯,g2k3,g2k32,⋯,g2k3t。
对于
x
=
g
k
(
m
o
d
p
)
x=g^k\ (mod\ \ p)
x=gk (mod p),如果它没有平方根,就代表k为奇数,此时给两边乘上一个g,就成了g的偶数次,也就存在了平方根,存在平方根后,判断存不存在四次方根,如果不存在就代表
k
=
2
,
6
,
10
,
14
⋯
k=2,6,10,14\cdots
k=2,6,10,14⋯,所以乘以
g
2
g^2
g2也就存在了四次方根。后面依次继续八次方根。。。其他次方根同理,但是判断次数及判断条件会发生变化,例如:三次方根不存在则k可能会是1,2,而乘一次是k+1,所以三次方根需要多判断一次,举一个本题的例子说:对于
g
2
3
3
g^{2^33}
g233,我们按顺序判断完其存在八次方根后,判断其存不存在24次方根,不存在则
k
=
8
,
16
k=8,16
k=8,16(我只枚举了最小的情况),此时我们给其乘以
g
8
g^8
g8,不能一次保证存在24次方根,所以需要再判断一次。
该题因此可以通过从
2
k
3
t
~
1
2^k3^t~1
2k3t~1判断,可以确认x是g的几次方。
也就求出了<3>中需要的
z
,
c
z,c
z,c。需要注意的是全程可能两数相乘爆long long,需要快速乘
代码
#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;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
const int maxn = 1e6;
const ll INF = 0x3f3f3f3f;
const int MOD = 1e9+7;
const int eps = 1e-8;
ll kk,tt;
ll qmul(ll a,ll b,ll p){
a%=p,b%=p;
ll ans=0;
while(b){
if(b&1) ans=(ans+a)%p;
a=(a+a)%p;
b>>=1;
}
return ans%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;
}
bool judge(ll x,ll p){
if(qpow(x,(p-1)/2,p)==1) return 0;
if(qpow(x,(p-1)/3,p)==1) return 0;
return 1;
}
void p_power(ll p){ //p-1=2^kk*3^tt
kk=0,tt=0;
ll pp = p-1;
while(pp%2 == 0){
kk++;
pp /= 2;
}
while(pp%3 == 0){
tt++;
pp /= 3;
}
}
ll get_power(ll x,ll g,ll p){
ll pp = p-1;
ll ans=0;
ll rec=g,po=1; //rec = g ^ po
for(int i = 0; i < kk; ++i){
pp /= 2;
if(qpow(x,pp,p) != 1){
ans += po;
x = qmul(x,rec,p);
}
rec = qmul(rec,rec,p);
po = po*2;
}
for(int i = 0; i < tt; ++i){
pp /= 3;
if(qpow(x,pp,p) != 1){
ans += po;
x = qmul(x,rec,p);
if(qpow(x,pp,p) != 1){
ans += po;
x = qmul(x,rec,p);
}
}
rec = qmul(rec,qmul(rec,rec,p),p);
po = po*3;
}
return p-1-ans;
}
ll exgcd(ll a,ll b,ll &x,ll &y){
if(b == 0){
x=1,y=0;
return a;
}
else{
ll d= exgcd(b,a%b,y,x);
y-=x*(a/b);
return d;
}
}
int main(){
#ifndef ONLINE_JUDGE
freopen("1in.in","r",stdin);
freopen("1out.out","w",stdout);
#endif
int t;cin>>t;
while(t--){
ll g,p,a,b;
scanf("%lld%lld%lld",&p,&a,&b);
for(int i= 2;;i++){ //求p的原根g
if(judge(i,p)){
g=i;break;
}
}
// cout<<"g:"<<g<<endl;
p_power(p); //p-1 = 2^k * 3^t
a=get_power(a,g,p);
b=get_power(b,g,p);
// cout<<"a:"<<a<<" "<<"b:"<<b<<endl;
if(a == 0 && b == 0) printf("1\n");
else if(a == 0) printf("-1\n");
else{
ll x,y;
ll d = exgcd(a,p-1,x,y); // ax=b (mod p)-> ax+(p-1)y = b
if(b%d) printf("-1\n");
else{
ll s = (p-1)/d;
if(x < 0) x += s;
x = qmul(x,b/d,s);
printf("%lld\n",x);
}
}
}
}
Pohlig Hellman:
#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;
}