sgu261：Discrete Roots(原根+离散对数+扩展欧几里得)

题目大意：
$~~~~~~$给你两个质数$p,k$，和整数$a(0\leq a <p)$，求解所有满足$x^k\equiv a(mod~p)$的$x$。

分析：
$~~~~~~$假设我们已经会了原根和离散对数$...$
$~~~~~~$设$p$的原根为$r$，$x\equiv r^i(mod~p),a\equiv r^j(mod~p)$，其中$j$可以通过大步小步求出来，问题就转化成$r^{ik}\equiv r^j(mod~p)\Rightarrow ik\equiv j(mod~(p-1))$，这个很显然可以用扩展欧几里得做，然后就没有然后了。

$($以下为学术内容，说明原根和离散对数的定义和求法，大神请忽略$)$
$~~~~~~$$1.$阶：
$~~~~~~~~~$对于两个互质的正整数$a,m$，定义$a$对模$m$的阶为最小的正整数$d$满足：
$~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~$$a^d \equiv 1 (mod~m)$
$~~~~~~~~~$记作$\delta_m(a)$。
$~~~~~~$$2.$原根：
$~~~~~~~~~$满足$\delta_m(a)=\varphi(m)$的$a$。
$~~~~~~$$3.$原根的求法：
$~~~~~~~~~$由于最小的原根通常很小，我们枚举原根$r$，满足下列两个即可：
$~~~~~~~~~$$①r^{\varphi(m)}\equiv 1(mod~m)$；
$~~~~~~~~~$$②$对于$\varphi(m)$的任意一个质因子$p$，满足$r^{\varphi(m)/p}\not \equiv 1(mod~m)$。
$~~~~~~$$4.$离散对数：
$~~~~~~~~~$对于整数$a,b,m$，其中$a,b$均与$m$互质，定义$log_{m,a}b$为最小的整数$d$满足：
$~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~$$a^d\equiv  b (mod~m)$
$~~~~~~$$5.$离散对数的求法$($大步小步$)$：
$~~~~~~~~~$设定一个参数$l$。
$~~~~~~~~~$将$a^0,a^1,...,a^{l-1}$存入哈希表。
$~~~~~~~~~$枚举$k:0,1,..., \left \lfloor \frac{\varphi(m)-1}{l} \right \rfloor$，判断是否有$a^{kl+j}\equiv b(mod~m),0\leq j <l$。
$~~~~~~~~~$即判断是否存在$0\leq j<l$满足$a^j\equiv ba^{-kl}(mod~m)$
$~~~~~~~~~$将$l$设为$\sqrt{\varphi(m)}$，时间复杂度就是$O(\sqrt{\varphi(m)})$。

AC code：

#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <algorithm>
#include <string>
#include <sstream>
#include <iostream>
#include <map>
#include <set>
#include <list>
#include <stack>
#include <queue>
#include <vector>
#define debug(...) fprintf(stderr, __VA_ARGS__)
#define pb push_back
#define mp make_pair
typedef long long LL;
typedef double DB;
typedef long double LD;
using namespace std;

const LL SP = 1e5;

LL p, k, a;

bool hashp[SP];
LL prime[SP], tot;

vector<LL> fact;

LL r, x, y;

map<LL,LL> Map;

vector<LL> ans;

void euler_prime(LL lim)
{
    for(int i = 2; i <= lim; ++i)
    {
        if(!hashp[i]) prime[++tot] = i;
        for(int j = 1; j <= tot && prime[j]*i <= lim; ++j)
        {
            hashp[i*prime[j]] = true;
            if(!(i%prime[j])) break;
        }
    }
}

void fact_decomposition(LL p)
{
    for(int i = 1; i <= tot; ++i)
    {
        if(!(p%prime[i]))
            fact.pb(prime[i]);
        while(!(p%prime[i]))
            p /= prime[i];
    }
    if(p > 1) fact.pb(p);
}

LL pow_mod(LL a, LL b, LL p)
{
    LL ret = 1, base = a;
    while(b)
    {
        if(b&1) ret = ret*base%p;
        base = base*base%p;
        b >>= 1;
    }
    return ret;
}

LL find_root(int p)
{
    LL ret = 0;
    for(int i = 1; i <= 200; ++i)
        if(pow_mod(i, p-1, p) == 1)
        {
            bool flag = true;
            for(int j = 0, sz = fact.size(); j < sz; ++j)
                if(pow_mod(i, (p-1)/fact[j], p) == 1)
                {
                    flag = false;
                    break;
                }
            if(flag) 
            {
                ret = i;
                break;
            }
        }
    return ret;
}

LL inv(LL a, LL p)
{
    return pow_mod(a, p-2, p);
}

LL baby_step_giant_step(LL a, LL b, LL p)
{
    LL l = sqrt(p-1), ret;
    for(int i = 0, j = 1; i < l; ++i)
    {
        Map[j] = i;
        j = (LL)j*a%p;
    }
    LL ia = pow_mod(inv(a, p), l, p), ta = 1;
    for(int k = 0, lim = (p-2)/l; k <= lim; ++k)
    {
        if(Map.find(b*ta%p) != Map.end())
        {
            ret = k*l+Map[(LL)b*ta%p];
            break;
        }
        ta = (LL)ta*ia%p;
    }
    return ret;
}

namespace math
{
    LL gcd(LL a, LL b)
    {
        while(b) b^=a^=b^=a%=b;
        return a;
    }
    void ext_gcd(LL a, LL b, LL &x, LL &y)
    {
        if(!b)
        {
            x = 1, y = 0;
            return ;
        }
        ext_gcd(b, a%b, x, y);
        LL tmp = x;
        x = y;
        y = tmp-a/b*y;
        return ;
    }
}

void solve(vector<LL> &ret, LL a, LL b, LL p)
{
    LL x, y, add;
    LL gcd = math::gcd(a, p);
    if(b%gcd) 
    {
        puts("0");
        exit(0);    
    }
    b /= gcd;math::ext_gcd(a, p, x, y);
    if(x < 0) x += ((-x-1)/p+1)*p;
    x *= b;add = p/gcd;x %= add;
    while(x < p) ret.pb(pow_mod(r, x, p+1)), x += add;
}

int main()
{
    #ifndef ONLINE_JUDGE
    freopen("input.txt", "r", stdin);
    freopen("output.txt", "w", stdout);
    #endif

    cin >> p >> k >> a;
    if(!a) 
    {
        puts("1");puts("0");
        return 0;
    }
    euler_prime(sqrt(p));
    fact_decomposition(p-1);
    r = find_root(p);
    y = baby_step_giant_step(r, a, p);
    solve(ans, k, y, p-1);
    cout << ans.size() << endl;
    sort(ans.begin(), ans.end());
    for(int i = 0, sz = ans.size(); i < sz-1; ++i)
        cout << ans[i] << ' ';
    cout << ans[ans.size()-1] << endl;

    #ifndef ONLINE_JUDGE
    fclose(stdin);
    fclose(stdout);
    #endif
    return 0;
}