hdu 3930 X^N=a(mod) p 求X

传送门

/*hdu 3930
  题意：
  给定newx, k, m, 方程 (x^k)%m=newx, 求在模m意义下的所有解x。
  限制：
  0 <= newx, m, k <= 1.5*10^15; m是素数。
  思路：
  N次剩余
 */
#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <vector>
#include <algorithm>
using namespace std;
#define LL __int64
#define PB push_back
LL mul(LL a,LL b,LL m)
{
	LL ret = 0;
	a %= m;
	while(b)
	{
		if(b & 1) ret = (ret + a) % m;
		a = (a + a) % m;
		b >>= 1;
	}
	return ret;
}
LL a_b_MOD_c(LL a,LL b,LL m)
{
	LL ret = 1;
	a %= m;
	while(b)
	{
		if(b&1) ret = mul(ret,a,m);
		a = mul(a,a,m);
		b >>= 1;
	}
	return ret;
}

LL ext_gcd(LL a,LL b,LL &x,LL &y)
{
	if(b==0)
	{
		x=1, y=0;
		return a;
	}
	LL ret= ext_gcd(b,a%b,y,x);
	y-= a/b*x;
	return ret;
}
vector<LL> a;
bool g_test(LL g,LL p)
{
	for(LL i=0; i<a.size(); ++i)
		if(a_b_MOD_c(g,(p-1)/a[i],p)==1)
			return 0;
	return 1;
}
LL pri_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 Node  
{  
    LL idx;  
    LL val;  
} baby[1000000];  
bool cmp(Node n1,Node n2)  
{  
    return n1.val!=n2.val?n1.val<n2.val:n1.idx<n2.idx;  
}  
LL gcd(LL a,LL b)  
{  
    return b==0?a:gcd(b,a%b);  
}  
LL inval(LL a,LL b,LL n)  
{  
    LL e,x,y;  
    ext_gcd(a,n,x,y);  
    e=((LL)x*b)%n;  
    return e<0?e+n:e;  
}  
LL PowMod(LL a,LL b,LL MOD)  
{  
    LL ret=1%MOD,t=a%MOD;  
    while(b)  
    {  
        if(b&1)  
            ret=((LL)ret*t)%MOD;  
        t=((LL)t*t)%MOD;  
        b>>=1;  
    }  
    return ret;  
}  
LL BinSearch(int num,int m)  
{  
    int low=0,high=m-1,mid;  
    while(low<=high)  
    {  
        mid=(low+high)>>1;  
        if(baby[mid].val==num)  
            return baby[mid].idx;  
        if(baby[mid].val<num)  
            low=mid+1;  
        else  
            high=mid-1;  
    }  
    return -1;  
}  
LL BabyStep(LL A,LL B,LL C)  
{  
    LL tmp,D=1%C;  
    LL temp;  
    for(LL i=0,tmp=1%C; i<100; i++,tmp=((LL)tmp*A)%C)  
        if(tmp==B)  
            return i;  
    LL d=0;  
    while((temp=gcd(A,C))!=1)  
    {  
        if(B%temp) return -1;  
        d++;  
        C/=temp;  
        B/=temp;  
        D=((A/temp)*D)%C;  
    }  
    LL m=(LL)ceil(sqrt((long double)C));  
    for(LL i=0,tmp=1%C; i<=m; i++,tmp=((LL)tmp*A)%C)  
    {  
        baby[i].idx=i;  
        baby[i].val=tmp;  
    }  
    sort(baby,baby+m+1,cmp);  
    int cnt=1;  
    for(int i=1; i<=m; i++)  
        if(baby[i].val!=baby[cnt-1].val)  
            baby[cnt++]=baby[i];  
    LL am=PowMod(A,m,C);  
    for(LL i=0; i<=m; i++,D=((LL)(D*am))%C)  
    {  
        LL tmp=inval(D,B,C);  
        if(tmp>=0)  
        {  
            LL pos=BinSearch(tmp,cnt);  
            if(pos!=-1)  
                return i*m+pos+d;  
        }  
    }  
    return -1;  
}  
/*n次剩余
  任务：
  给定N, a, p, 求出(x^N)%p=a 在模p意义下的所有解x。
  说明：
  令g为p的原根，因为p为素数，所以phi(p)=p-1。
  由原根的性质得：
  如果g为p的原根，则：g^i mod p != g^j mod p (p为素数), 其中i != j且i, j介於1至(p-1)之间
  所以，可以设g^y=x, g^t=a，则有：
  g^(y*N)%p=g^t
  又由原根的性质：
  g^(y*N)%p=g^t -> (y*N)%(p-1)=t (此方程可以由拓展欧几里得解)
  另外g^t=a可以由离散对数求出
 */
vector<LL> residue(LL p, LL N, LL a)
{
	LL g = pri_root(p);
	g %= p;
	LL m = BabyStep(g, a, p);
	vector<LL> ret;
	if(a == 0)
	{
		ret.PB(0);
		return ret;
	}
	if(m == -1)
		return ret;
	LL A = N, B = p - 1, C = m, x, y;
	LL d = ext_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.PB(a_b_MOD_c(g, x, p));
	}
	sort(ret.begin(), ret.end());
	ret.erase(unique(ret.begin(), ret.end()), ret.end());
	return ret;
}
int main()
{
	int cas = 0;
	LL k,m,newx;
	while(scanf("%I64d%I64d%I64d",&k, &m, &newx)!=EOF)
	{
		vector<LL> ans;
		ans = residue(m,k,newx);
		printf("case%d:\n",++cas);
		if(ans.size()==0) puts("-1");
		for(int i = 0; i < ans.size(); ++i)
			printf("%I64d\n",ans[i]);
	}
	return 0;
}