51nod 1123 X^A Mod B 问题

先放在这好了。。。模2的次幂那一块还有问题，待改正

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<map>
using namespace std;
typedef long long LL;
inline int read()
{
int x=0;bool f=0;char c=getchar();
for (;c<'0'||c>'9';c=getchar()) f=c=='-'?1:0;
for (;c>='0'&&c<='9';c=getchar()) x=x*10+c-'0';
return f?-x:x;
}
const int N=100005;
int pr
,cnt,pm
,p
,pa
,alf
,n,A,B,C,r1
,
head
,nxt
,g,ans
,ansc,res
,resc,ndg=0;

inline int pw(int a,int k,int m){
int rec=1;
for (;k;a=1ll*a*a%m,k>>=1)
if (k&1) rec=1ll*rec*a%m;
return rec;
}

void init()
{
for (int i=2;i<N;i++)
{
if (!pm[i]) pr[++cnt]=i;
for (int j=1,k=i<<1;j<=cnt&&k<N;k=i*pr[++j])
{
pm[k]=pr[j];
if (i%pr[j]==0) break;
}
}
}

int gcd(int a,int b){
for (int c;b;c=a,a=b,b=c%b);
return a;
}

LL _x,_y,_b;
void _gcd(int a,int b){
if (!a) {_y=1;return;}
if (!b) {_x=1;return;}
_gcd(b,a%b);
LL t=_x;_x=_y;_y=t-_y*(a/b)%_b;
}
int rev(int a,int b){
_gcd(a,_b=b);
return (_x%b+b)%b;
}

int dv
,dc,ph,sq,_p,m,b,_B,b_p,_A;
map<int,int> cun;
map<int,int>::iterator it;
void solve(int id)
{
_p=p[id];m=pa[id];
dc=0;ph=pa[id]-pa[id]/p[id];
_B=B%m;_A=A%ph;b_p=0;

while (_B&&_B%_p==0) _B/=_p,b_p++;
if (b_p%A) {ndg=1;return;}
b_p/=A;

if (_B&&_p>2)
{
for (sq=2;sq*sq<=ph;sq++)
if (ph%sq==0)
{
dv[++dc]=sq;
if (sq*sq!=ph) dv[++dc]=ph/sq;
}
sort(dv+1,dv+dc+1);
for (g=2;;g++) if (gcd(g,m)==1)
{
bool ok=1;
for (int i=1;i<=dc&&ok;i++)
if (pw(g,dv[i],pa[id])==1) ok=0;
if (ok) break;
}
cun.clear();
for (int i=0,s=1,g_sq=pw(g,sq,m);i<ph;i+=sq,s=1ll*s*g_sq%m)
cun[s]=i;

for (int i=0,s=1,g_rv=pw(g,ph-1,m);i<sq;i++,s=1ll*s*g_rv%m)
if ((it=cun.find(1ll*_B*s%m))!=cun.end())
{b=i+it->second;break;}

int A_ph=gcd(A,ph);
if (b%A_ph) {ndg=1;return;}
_A/=A_ph;b/=A_ph;ph/=A_ph;
_A=1ll*rev(_A,ph)*b%ph;
int p_b_p=pw(_p,b_p,m);
for (int i=0;i<A_ph;i++)
ans[++ansc]=1ll*pw(g,_A,m)*p_b_p%m,_A+=ph,nxt[ansc]=head[id],head[id]=ansc;
}
else if (_B&&_p==2)
{
int cc=0;
for (int i=1;i<m;i+=2)
if (pw(i,_A,m)==_B)
ans[++ansc]=1ll*i*pw(2,b_p,m)%m,nxt[ansc]=head[id],head[id]=ansc,cc++;
if (!cc) {ndg=1;return;}
}
else
{
dc=0;
for (int i=1;i<m;i++)
if (gcd(i,m)==1) dv[++dc]=i;
for (int i=1;i<=dc;i++)
for (int j=(alf[id]-1)/A+1,val=pw(_p,j,m);j<alf[id]&&1ll*val*dv[i]<m;j++,val*=_p)
ans[++ansc]=1ll*val*dv[i]%m,nxt[ansc]=head[id],head[id]=ansc;
ans[++ansc]=0,nxt[ansc]=head[id],head[id]=ansc;
}
}

void dfs(int d,int a)
{
if (d>n) {res[++resc]=a;return;}
for (int i=head[d];i;i=nxt[i])
dfs(d+1,(a+1ll*ans[i]*r1[d])%C);
}

int main()
{
init();
for (int cas=read();cas--;)
{
memset(head,0,sizeof head);
A=read();C=read();B=read();n=0;ansc=resc=0;
int _C=C;ndg=0;
for (int i=1;i<=cnt&&pr[i]<=_C;i++)
if (_C%pr[i]==0)
{
p[++n]=pr[i];_C/=pr[i];pa
=pr[i];alf
=1;
while (_C%pr[i]==0) _C/=pr[i],pa
*=pr[i],alf
++;
}
if (_C>1) p[++n]=_C,pa
=_C,alf
=1;
for (int i=1;i<=n&&!ndg;i++)
solve(i);
if (ndg) {puts("No Solution");continue;}

for (int i=1;i<=n;i++)
{
int a=rev(C/pa[i],pa[i]),b=C/pa[i];
r1[i]=1ll*C/pa[i]*rev(C/pa[i],pa[i])%C;
}
dfs(1,0);
sort(res+1,res+resc+1);
for (int i=1;i<=resc;i++) printf("%d ",res[i]);
puts("");
}
return 0;
}