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

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

代码：

#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
#include<algorithm>
#define LL long long
using 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 \times indx \equiv indb \ (mod\ φ(p))$解的个数为 $gcd(a,φ(p))$个？

一般化：　
$a \times x \equiv b \ (mod\ c)$
$a \times x +c \times y = b$
$x_0+k\times \frac{c}{gcd(a,c)}$为通解
因为 $x<c$
所以解的个数为 $\frac{c}{\frac{c}{gcd(a,c)}}=gcd(a,c)$