bzoj1319&1420 Sgu261 Discrete Roots （原根+BSGS）

bzoj1319&1420 Sgu261 Discrete Roots

原题地址：
http://www.lydsy.com/JudgeOnline/problem.php?id=1319
http://www.lydsy.com/JudgeOnline/problem.php?id=1420

题意：
给出三个整数p,k,a,其中p为质数，求出所有满足x^k=a (mod p),0<=x<=p-1的x。

数据范围
0 < = a < p < = 10^9, 2 < = k < = 100000

题解：
xk≡a (mod p)
由于p为质数，可以求出p的原根。
1.a=0：
那么x只能为p的倍数，又因为0<=x<=p-1，则只有一个解 0。
2.a!=0：
因为p为质数，a与p互素，则a可以用ginda表示
同时，x也与p互素。
由 (gindx)a≡gindb (mod p)
得 a×indx≡indb (mod φ(p))
由于gindx唯一确定一个x，
所以上述a×indx≡indb (mod φ(p))解的数量就是原方程解的数量
快速幂一下得到每个x。

代码：

#include<cstdio>#include<iostream>#include<cmath>#include<cstring>#include<algorithm>#define LL long longusing namespace std;const int N=100005;const LL S=100003;LL k,a,p,phi,g,inda,indx,top;LL stack[N];LL modpow(LL A,LL B,LL mod){    LL base=A; LL ans=1LL;    for(;B;B>>=1)    {        if(B&1) ans=(ans*base)%mod;        base=(base*base)%mod;    }    return ans;}struct H{    LL head[S],dest[S][2];    int nxt[N],num;    void init() {memset(head,0,sizeof(head)); num=0;}    void insert(LL A,LL B)    {        LL key=B%S;        for(int i=head[key];i;i=nxt[i])        if(dest[i][1]==B) return;        num++;        dest[num][0]=A;        dest[num][1]=B;        nxt[num]=head[key];        head[key]=num;    }    LL find(LL val)    {        LL key=val%S;        for(int i=head[key];i;i=nxt[i])        if(dest[i][1]==val) return dest[i][0];        return -1;    }}Hash;LL BSGS(LL A,LL B,LL P){    Hash.init(); LL cur=1; LL C=sqrt(P-1)+1;    for(LL i=0;i<C;i++)    {        if(cur==B) return i;        Hash.insert(i,cur);        cur=cur*A%P;    }    LL base=modpow(cur,P-2,P);    cur=B*base%P;    for(LL i=C;i<P;i=i+C)    {        LL ret=Hash.find(cur);        if(ret!=-1) return ret+i;        cur=cur*base%P;    }    return -1;}LL getG(LL m){    LL x=phi; top=0;    for(LL i=2;i*i<=x;i++)    {        if(x%i==0)        {            stack[++top]=i;            while(x%i==0) x=x/i;        }    }    if(x>1) stack[++top]=x;    for(LL i=1;i<=m;i++)    {        int j;         for(j=1;j<=top;j++)        if(modpow(i,phi/stack[j],m)==1) break;        if(j>top) return i;    }    return -1;}LL getgcd(LL A,LL B) {return B==0? A: getgcd(B,A%B);}void exgcd(LL A,LL B,LL &x,LL &y){    if(B==0) { x=1; y=0; return;}    LL x0,y0;    exgcd(B,A%B,x0,y0);    x=y0;    y=x0-(A/B)*y0;}int main(){    scanf("%lld%lld%lld",&p,&k,&a); phi=p-1;    if(a==0) { printf("1\n0\n"); return 0;}    g=getG(p);      inda=BSGS(g,a,p);           LL gcd=getgcd(k,p-1);    if(inda%gcd!=0) { printf("0\n"); return 0;}    LL A=k/gcd; LL B=(p-1)/gcd;    LL x,y; exgcd(A,B,x,y);     x=(x+B)%B; x=(x*(inda/gcd)+B)%B;    top=0;    for(;x<p-1;x+=B) stack[++top]=modpow(g,x,p);    sort(stack+1,stack+top+1);    printf("%lld\n",top);    for(int i=1;i<=top;i++)     printf("%lld\n",stack[i]);    return 0;}

---

补充：

为什么 a×indx≡indb (mod φ(p))解的个数为 gcd(a,φ(p))个？

一般化：　
a×x≡b (mod c)
a×x+c×y=b
x0+k×cgcd(a,c)为通解
因为 x<c
所以解的个数为 ccgcd(a,c)=gcd(a,c)