bzoj1420&bzoj1319 Discrete Root

题目链接：bzoj1420
题目大意：
已知k，a，p，求x^k=a (mod p)的所有根（根的范围[0,p-1]
Input
三个整数p,k,a。
Output
第一行一个整数，表示符合条件的x的个数。 第二行开始每行一个数，表示符合条件的x,按从小到大的顺序输出。

题解
数论
$x^k\equiv a(\mod p)$
设p的原根为g。(下面的I(x)就是x的指标
有$g^{I(x)}\equiv x(\mod p)$；$g^{I(a)}\equiv a(\mod p)$
则原方程变为$(g^{I(x)})^k\equiv g^{I(a)}(\mod p)$
即$k\times I(x)\equiv I(a) (\mod \phi(p))$
然后就可以用扩展欧几里德求I(x)的所有解。
代进原方程求出所有0到p-1范围里的根就好了
题目没有说，但是据discuss里的、bzoj1319和sgu里的原题来看，p应给是个质数。
a=0的情况下要特判

心好累啊WA了4次才发现是拓欧里的a/b忘了加括号让它向下取整QwQ
果然一点细节都不能有错啊

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long LL;
#define maxn 101000

struct node
{
    LL j,x;
    bool operator < (const node y)const {return x<y.x;}
}aj[maxn];
LL ans[maxn],pri[maxn],cnt;
void get_pri(LL pn)
{
    for (LL i=2;i*i<=pn;i++)
     if (pn%i==0)
     {
         pri[++cnt]=i;
         while (pn%i==0) pn/=i;
     }
    if (pn!=1) pri[++cnt]=pn;
}
void exgcd(LL a,LL b,LL &d,LL &x,LL &y)
{
    if (b==0)
    {
        d=a;x=1;y=0;
        return;
    }
    LL xx,yy;exgcd(b,a%b,d,xx,yy);
    x=yy;y=xx-yy*(a/b);
}
LL qpow(LL x,LL t,LL mod)
{
    LL ret=1;
    while (t)
    {
        if (t&1) ret=1LL*ret*x%mod;
        x=(1LL*x*x)%mod;t>>=1;
    }return ret;
}
LL twocut(LL r,LL x)
{
    LL l=0,ret=-1;
    while (l<=r)
    {
        LL mid=(l+r)>>1;
        if (aj[mid].x==x) {ret=aj[mid].j;break;}
        if (aj[mid].x<x) l=mid+1;
        else r=mid-1;
    }
    return ret;
}
LL BSGS(LL a,LL b,LL p)//a^x=b(mod p)
{
    LL i,m=ceil(sqrt(p));
    for (i=0;i<m;i++)
    {
        aj[i].j=i;
        aj[i].x=qpow(a,i,p)%p;
    }sort(aj,aj+m);
    LL am=qpow(qpow(a,p-2,p),m,p);
    LL ret=-1,num=0,k=b;aj[num]=aj[0];
    for (i=1;i<m;i++) if (aj[i].x!=aj[i-1].x) aj[++num]=aj[i];
    for (i=0;i<m;i++)
    {
        LL w=twocut(num,k);
        if (w!=-1) {ret=w+i*m;break;}
        k=((LL)(k*am))%p;
    }
    return ret;
}
LL get_r(LL p)
{
    if (p==2) return 1;
    cnt=0;get_pri(p-1);
    LL g=2;
    while (g<p)
    {
        bool bk=false;
        for (LL i=1;i<=cnt;i++)
         if (qpow(g,(p-1)/pri[i],p)==1) {bk=true;break;}
        if (!bk) break; g++;
    }
    if (g<p) return g;
    return 0;
}
int main()
{
    //freopen("a.in","r",stdin);
    //freopen("a.out","w",stdout);
    LL p,k,a,g,i;
    scanf("%lld%lld%lld",&p,&k,&a);
    if (a==0) {printf("1\n0\n");return 0;}
    g=get_r(p);
    LL Ix,Ia,yy,d;
    Ia=BSGS(g,a,p);
    //k*Ix+phi(p)*u=Ia
    exgcd(k,p-1,d,Ix,yy);
    if (Ia%d!=0) printf("0\n"); 
    else
    {
        LL sp=(p-1)/d,tot=0;
        Ix*=Ia/d%sp;Ix=(Ix%sp+sp)%sp;
        if (Ix==0) Ix=sp;
        for (i=Ix;i<p;i+=sp) ans[++tot]=qpow(g,i,p);
        sort(ans+1,ans+1+tot);
        printf("%lld\n",tot);
        for (i=1;i<=tot;i++) printf("%lld\n",ans[i]);
    }
    return 0;
}