SGU 261 Discrete Roots（原根+高次同余方程+线性同余方程）

Description
给出三个整数p,k,a，输出方程x^k=a(mod p)的所有解
Input
三个整数p,k,a (2<=p<=10^9,2<=k<=100000,0<=a< p)
Output
首先输出方程解的个数，然后升序输出所有解（解的范围为[0,p-1]）
Sample Input
11 3 8
Sample Output
1
2
Solution
首先暴力枚举求出p的最小原根root，那么root^t构成p的完全剩余系，即我们可以用root^t表示出a和x，设x=root^x*，a=root^b，解高次同余方程root^b=a(mod p)可以解出b，那么原方程化为root^(kx*)=root^b(mod p)，根据原根的性质知这个方程与方程kx*=b(mod p-1)等价，所以解这个线性同余方程可以得到所有的x*，那么root^x*即为所有的解，注意a=0时需要特判
Code

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<map>
#include<vector>
using namespace std;
#define maxn 2222
#define INF 0x3f3f3f3f
typedef long long ll;
vector<ll>vec;
ll mod_pow(ll a,ll b,ll p)//快速幂，求a^b(mod p)
{
    ll ans=1ll;
    a%=p;
    while(b)
    {
        if(b&1) ans=(ans*a)%p;
        a=(a*a)%p;
        b>>=1;
    } 
    return ans;
}
ll find_root(ll p)//求p的最小原根
{
    ll factor[maxn];
    int res=0;
    ll pp=p-1;
    for(int i=2;i*i<=pp;i++)
        if(pp%i==0)
        {
            factor[res++]=i;
            while(pp%i==0)
                pp/=i;
        }
    if(pp!=1)factor[res++]=pp;
    for(ll k=1;;k++)
    {
        int flag=1;
        for(int i=0;i<res&&flag;i++)
            if(mod_pow(k,(p-1)/factor[i],p)==1)
                flag=0;
        if(flag) return k;
    }
}
ll mul_linear(ll a,ll b,ll p)//求解高次同余方程a^x=b(mod p) 
{
    ll n=ceil(sqrt((double)(p+1)));
    map<ll,ll>Map;
    Map[1]=n; 
    ll m=mod_pow(a,p-2,p); 
    m=mod_pow(m,n,p); 
    ll temp=1;
    for(ll i=1;i<n;i++)
    {
        temp=(temp*a)%p;
        if(!Map[temp])
            Map[temp]=i;
    }
    for(ll i=0;i<n;i++)
    {
        if(Map[b])
        {
            ll r=Map[b];
            Map.clear();
            return i*n+(r==n?0:r);
        }
        b=(b*m)%p;
    }
    return -1;
}
ll extend_gcd(ll a,ll b,ll &x,ll &y)//扩展欧几里得
{ 
    ll d=a;
    if(b!=0)
    {
        d=extend_gcd(b,a%b,y,x);
        y-=(a/b)*x;
    }
    else 
    {
        x=1;
        y=0;
    }
    return d;
}
void linear(ll a,ll b,ll c)//求解一元线性同余方程ax=b(mod c)
{
    ll x,y;
    ll g=extend_gcd(a,c,x,y);
    if(b%g)
        return ;
    x=(x%c+c)%c;
    x=x*(b/g)%(c/g);
    for(ll i=0;i<g;i++)
        vec.push_back(((x+i*c/g)%c));
}
int main()
{
    ll p,k,a;
    while(~scanf("%lld%lld%lld",&p,&k,&a))
    {
        vec.clear();
        if(a==0)
        {
            printf("1\n0\n");
            continue;
        }
        ll root=find_root(p);
        ll b=mul_linear(root,a,p);
        linear(k,b,p-1);
        printf("%d\n",vec.size());
        for(int i=0;i<vec.size();i++)
            vec[i]=mod_pow(root,vec[i],p);
        sort(vec.begin(),vec.end());
        for(int i=0;i<vec.size();i++)
            printf("%lld%c",vec[i],i==vec.size()-1?'\n':' ');
    }
    return 0;
}