HDU 6632 discrete logarithm problem （Pohlig-hellman算法）-优快云博客

本文链接：https://blog.youkuaiyun.com/oampamp1/article/details/104077651

原题题面

给定三个数 $p,a,b(p∈prime,65537≤p≤10^{18}, 2≤a,b≤p−1)$ , $p - 1$ 的质因子只有 $2$ 和/或 $3$ . 请算出最小的 $x$ 使得 $a^x\equiv b(mod\ p)$ ,若无解输出 $- 1$ .

样例输入

6
65537 2 3
65537 2 4
65537 3 4
65537 4 4
65537 5 4
666334875701477377 2 3

样例输出

-1
2
45056
1
36864
1957714645490451

题面分析

对于离散对数问题，我们最常用的是 $B S G S$ 算法，其复杂度为 $O(\sqrt p)$ 。
但题中的 $p$ 最大可以是 $1 e 18$ 左右，因此不能采用这种方法。
注意到 $p - 1$ 的质因子很少，可以采用 $p o h l i g - h e l l m a n$ 算法，具体可以参考我写的这篇文章：pohlig-hellman算法

AC代码(93ms)

#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=1;;g++)
    {
        int flag=1;
        for (int j=1;j<=numprime;j++)
        {
            if (quick_pow (g,(p-1)/prime[j],p)==1)
            {
                flag=0;
                break;
            }
        }
        if (flag)
            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
        {
            while (xx<0)
            {
                xx+=tt;
            }
            xx=quick_mul (xx,(x2/rr),tt);
        }
        if (xx==0)
            xx=tt;
        printf ("%lld\n",xx);
    }
}

后记

改进了之前用了笨笨的快速乘之后，从2.5s直接降到0.07s，在看到题目时间限制1.5s，整个人直接傻掉…
一道半年后才来切掉的 $p o h l i g - h e l l m a n$ 裸题 ~~（君子报仇，半年不晚）~~ ，学懂算法就可以了。
不过，标程给了个和原根相关的算法，也可以做。
DrGilbert 2020.1.23