SGU 261 Discrete Roots N次余下

SGU 261 Discrete Roots N次剩余

链接：vjudge

题意：给出两个素数P，K (2 <= P <= 10^9, 2 <= K <= 100000)和一个整数A (0 <= A < P) ，找出 x ^ K ≡ A mod P的解。

思路：N次剩余，模板题，复杂度O(sqrt(p))。

代码：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <map>
#include <cstdlib>
#include <queue>
#include <stack>
#include <vector>
#include <ctype.h>
#include <algorithm>
#include <string>
#include <set>
#define PI acos(-1.0)
#define maxn 1000005
#define INF 0x7fffffff
#define eps 1e-8
typedef long long LL;
typedef unsigned long long ULL;
using namespace std;
LL extend_gcd(LL a, LL b, LL &x, LL &y)
{
    if(b==0)
    {
        x=1;
        y=0;
        return a;
    }
    LL r=extend_gcd(b,a%b,x,y);
    LL t=x;
    x=y;
    y=t-a/b*y;
    return r;
}
LL pow_mod(LL aa,LL ii,LL nn)
{
    if(ii==0)
        return 1%nn;
    LL temp=pow_mod(aa,ii>>1,nn);
    temp=temp*temp%nn;
    if(ii&1)
        temp=temp*aa%nn;
    return temp;
}
vector <LL> a;
bool g_test(LL g,LL p)
{
    for(LL i=0; i<a.size(); i++)
        if(pow_mod(g,(p-1)/a[i],p)==1)
            return 0;
        return 1;
}
LL primitive_root(LL p)
{
    a.clear();
    LL tmp=p-1;
    for(LL i=2; i<=tmp/i; i++)
        if(tmp%i==0)
        {
            a.push_back(i);
            while(tmp%i==0)
                tmp/=i;
        }
    if(tmp!=1)
        a.push_back(tmp);
    LL g=1;
    while(true)
    {
        if(g_test(g,p))
            return g;
        g++;
    }
}
struct b_step
{
    int i,m;
} bb[maxn];
LL inval(LL a,LL b,LL n)
{
    LL e,x,y;
    extend_gcd(a,n,x,y);
    e=((LL)x*b)%n;
    return e<0?e+n:e;
}
bool cmp(b_step a,b_step b)
{
    return a.m==b.m?a.i<b.i:a.m<b.m;
}
int BiSearch(int m,LL num)
{
    int low=0,high=m-1,mid;
    while(low<=high)
    {
        mid=(low+high)>>1;
        if(bb[mid].m==num)
            return bb[mid].i;
        if(bb[mid].m<num)
            low=mid+1;
        else
            high=mid-1;
    }
    return -1;
}
LL giant_step_baby_step(LL b,LL n,LL p)
{
    LL tt=1%p;
    for(int i=0; i<100; i++)
    {
        if(tt%p==n)
            return i;
        tt=((LL)tt*b%p);
    }
    LL D=1%p;
    int d=0,temp;
    while((temp=__gcd(b,p))!=1)
    {
        if(n%temp)
            return -1;
        d++;
        n/=temp;
        p/=temp;
        D=((b/temp)*D)%p;
    }
    int m=(int)ceil(sqrt((double)(p)));
    bb[0].i=0,bb[0].m=1%p;
    for(int i=1; i<=m; i++)
    {
        bb[i].i=i;
        bb[i].m=bb[i-1].m*b%p;
    }
    sort(bb,bb+m+1,cmp);
    int top=1;
    for(int i=1; i<=m; i++)
        if(bb[i].m!=bb[top-1].m)
            bb[top++]=bb[i];
    int bm=pow_mod(b,m,p);
    for(int i=0; i<=m; i++)
    {
        int tmp=inval(D,n,p);
        if(tmp>=0)
        {
            int pos=BiSearch(top,tmp);
            if(pos!=-1)
                return i*m+pos+d;
        }
        D=((LL)(D*bm))%p;
    }
    return -1;
}
vector <LL>residue (LL p,LL N,LL a)
{
    LL g=primitive_root(p);
    LL m=giant_step_baby_step(g,a,p);
    vector <LL>ret;
    if(a==0)
    {
        ret.push_back(0);
        return ret;
    }
    if(m==-1)
        return ret;
    LL A=N,B=p-1,C=m,x,y;
    LL d=extend_gcd(A,B,x,y);
    if(C%d!=0)
        return ret;
    x=x*(C/d)%B;
    LL delta =B/d;
    for(int i=0;i<d;i++)
    {
        x=((x+delta)%B+B)%B;
        ret.push_back((LL)pow_mod(g,x,p));
    }
    sort(ret.begin(),ret.end());
    ret.erase(unique(ret.begin(),ret.end()),ret.end());
    return ret;
}
int main()
{
    LL p,N,a;
    scanf("%I64d%I64d%I64d",&p,&N,&a);
    vector<LL>ans=residue(p,N,a);
    printf("%d\n",ans.size());
    for(int i=0;i<ans.size();i++)
        printf("%I64d\n",ans[i]);
        return 0;
}