考前学算法，就是这么浪。
求 (x^aequiv bpmod p) 在 ([0,p)) 内的所有解。
核心思想是通过原根将幂变成指数。但是普通的模数 (p) 不一定有原根啊，那考虑使用分解质因数 (p=prod p_i^{alpha_i})，奇质数的幂是有原根的。拆开的理论依据是什么？换个元，令 (y=x^a)，那么 (yequiv bpmod p) 可以由若干个 (yequiv bpmod{p_i^{alpha_i}}) 联立而成的方程组等价表达，根据的是模数两两互质的 ExCRT（其实就是 CRT，可惜我不会写）。那么我们得到 (yequiv bpmod{p_i^{alpha_i}}) 的解集 (X_i) 后，其实就是 (x) 是原方程的解当且仅当（方程和方程组的等价性已经用 CRT 证明）对所有 (i) 都有 (xin X_ipmod{p_i^{alpha_i}})，于是我们可以从每个 (X_i) 各挑一个值出来用 CRT 合并，得到 (prod|X_i|) 个 ([0,p)) 内的 (x) 便是原方程的所有解。这里使用了两次 CRT，正着一次反着一次，步步为等价变换。
下面讨论 (x^aequiv bpmod{p^alpha}(pinmathbb P)) 怎么解。注意到当 (p eq 2) 时 (p^{alpha}) 必有原根，我们分成 (poverset{?}=2) 两种情况讨论，先考虑 (p eq 2)。首先如果 (b=0)，那么显然解集为 (p^{leftlceilfrac alpha a ight ceil}mid x)，下面 (b eq0)。若 (bperp p^{alpha})，求得 (p^{alpha}) 的任意一个原根（可以用 1/4 次 log 的算法求最小原根）(g)，然后用 BSGS 求出 (b=g^eta) 的 (left[0,varphi=varphi!left(p^alpha ight) ight)) 内唯一 (eta)。此时显然 (x) 是解仅当 (xperp p^alpha)，于是令 (x=g^y(yin[0,varphi)))，求出 (y) 在 (modvarphi) 意义下所有解，就相当于求出了 (x) 的所有解。原方程即 (g^{ay}equiv g^etapmod{p^alpha})，跟据原根的幂周期为 (varphi)，这显然等价于 (ayequivetapmod{varphi})。线性同余方程就没什么讲的必要了吧，exgcd 找特解然后枚举通解在 ([0,varphi)) 内所有取值即可。
如果 (b otperp p^alpha) 呢？此时 (b) 在模意义下没有原根。设 (b=p^eta q)，令 (x=p^gamma r)，那么原方程等价于 (p^{agamma}r^aequiv p^eta qpmod{p^alpha})。此时我们发现，由于 (b eq 0)，也就是 (eta<+infty) 即 (eta<alpha)，那么 (p^{eta} q) 所在剩余系中都恰好包含 (eta) 个质因子 (p)，那么我们需要 (agamma=eta)（并不是模意义下，而是严格相等！）！那么判一手 (amideta)，直接得到 (gamma=dfraceta a)，它现在是常数了。跟据同余式的性质将三边的 (p^eta) 约掉，得到 (r^aequiv qpmod{p^{alpha-eta}})，此时有 (p mid q)，于是归约到了 (bperp p^alpha) 的情况。求出 (rmod p^{alpha-eta}) 的所有可能值后，对于其中每一个，我们要找出对应的所有 (left[0,p^alpha ight)) 内的 (p^gamma r)，那么找到所有 (rinleft[0,p^{alpha-gamma} ight)) 即可，随便遍历一波，有 (p^{(alpha-gamma)-(alpha-eta)}=p^{eta-gamma}) 个值。
最后我们来考虑 (p=2) 的情形。由于没有原根，很可惜，这种情况只能暴力（有原根的情况下复杂度显然是 (|X|+mathrm O!left(sqrt {p^alpha} ight))，而 (|X|) 上限是有保证的）。直接暴力枚举 ([0,2^alpha)) 内所有 (x) 显然会 T，我们可以剪枝呀！考虑对 (alpha) 进行迭代，假设 (x^aequiv bpmod{2^{alpha-1}}) 的解集已经算出为 (X_{alpha-1})，那么 (x) 是 (x^aequiv bpmod{2^{alpha}}) 的解的必要条件是满足 (alpha-1) 的方程，即 (xin X_{alpha-1}pmod{2^{alpha-1}})，而对于每个 (x_0in X_{alpha-1})，在 (mod 2^alpha) 意义下有两个对应值，都快速幂判一下即可。这样最坏复杂度依然是跑满的，但是剪枝效果很好（可以通过模板题），毕竟 (2^eta) 时减掉一个枝相当于直接砍掉 (2^{alpha-eta}) 个可能的答案。
code（写起来非常爽啊）
#include<bits/stdc++.h>
using namespace std;
#define int long long
#define pb push_back
#define mp make_pair
#define X first
#define Y second
int qpow(int x,int y,int p){
int res=1;
while(y){
if(y&1)res=res*x%p;
x=x*x%p;
y>>=1;
}
return res;
}
int exgcd(int a,int b,int &x,int &y){
if(!b)return x=1,y=0,a;
int d=exgcd(b,a%b,y,x);
return y-=a/b*x,d;
}
int inv(int a,int p){
int x,y,d=exgcd(a,p,x,y);
if(d==1)return 0;
return (x%p+p)%p;
}
int excrt(int a1,int a2,int p1,int p2){
int k1,k2,d=exgcd(p1,p2,k1,k2);
assert(d==1);
k1=(k1%(p1*p2)*(a2-a1)%(p1*p2)+p1*p2)%(p1*p2);
return (a1+p1*k1)%(p1*p2);
}
vector<int> solve2(int a,int b,int alpha){
b%=1<<alpha;
if(alpha==0)return vector<int>({0});
vector<int> v=solve2(a,b,alpha-1),ret;
for(int i=0;i<v.size();i++){
if(qpow(v[i],a,1<<alpha)==b)ret.pb(v[i]);
if(qpow(v[i]+(1<<alpha-1),a,1<<alpha)==b)ret.pb(v[i]+(1<<alpha-1));
}
return ret;
}
int proot(int p,int alpha,int palpha){
int phi=palpha/p*(p-1);
vector<int> dsr;
for(int i=2;i*i<=phi;i++)if(phi%i==0){
dsr.pb(i);
while(phi%i==0)phi/=i;
}
if(phi>1)dsr.pb(phi);
phi=palpha/p*(p-1);
for(int i=2;;i++)if(i%p){
bool flg=true;
for(int j=0;j<dsr.size();j++)flg&=qpow(i,phi/dsr[j],palpha)!=1;
if(flg)return i;
}
}
int bsgs(int g,int b,int p){
unordered_map<int,int> mmp;
int B=sqrt(p)+1,now=1,now0=1;
for(int i=1;i<=B;i++)now=now*g%p,mmp[now*b%p]=i;
for(int i=1;i<=B;i++){
now0=now0*now%p;
int j=mmp[now0];
if(j)return i*B-j;
}
}
vector<int> solve(int a,int b,int p,int alpha){
int palpha=qpow(p,alpha,1e12);
// cout<<palpha<<"!!
";
b%=palpha;
if(p==2)return solve2(a,b,alpha);
if(b==0){
int div=(alpha+a-1)/a,fst=qpow(p,div,1e12);
vector<int> ret;
for(int i=0;i<palpha;i+=fst)ret.pb(i);
return ret;
}
if(b%p==0){
int beta=0,q=b;
while(q%p==0)beta++,q/=p;
if(beta%a)return vector<int>();
vector<int> v=solve(a,q,p,alpha-beta);
int gamma=beta/a;
vector<int> ret;
int lim=qpow(p,beta-gamma,1e12),mul=qpow(p,alpha-beta,1e12),pgamma=qpow(p,gamma,1e12);
for(int i=0;i<v.size();i++)for(int j=0;j<lim;j++)ret.pb((v[i]+j*mul)%palpha*pgamma%palpha);
return ret;
}
int g=proot(p,alpha,palpha),beta=bsgs(g,b,palpha),phi=palpha/p*(p-1);
int x,y,d=exgcd(a,phi,y,x);
if(beta%d)return vector<int>();
int mod=phi/d,gap=qpow(g,mod,palpha);y=(y*(beta/d)%mod+mod)%mod;
int now=qpow(g,y,palpha);
vector<int> ret;
for(int i=0;y+i*mod<phi;i++)ret.pb(now),now=now*gap%palpha;
return ret;
}
void mian(){
int a,p,b;
cin>>a>>p>>b;
// cerr<<a<<" "<<b<<" "<<p<<":
";
vector<pair<int,int> > dsr;
for(int i=2;i*i<=p;i++)if(p%i==0){
int cnt=0;
while(p%i==0)cnt++,p/=i;
dsr.pb(mp(i,cnt));
}
if(p>1)dsr.pb(mp(p,1));
vector<int> ans;ans.pb(0);
int now=1;
for(int i=0;i<dsr.size();i++){
vector<int> v=solve(a,b,dsr[i].X,dsr[i].Y),nw;
for(int j=0;j<ans.size();j++)for(int k=0;k<v.size();k++)nw.pb(excrt(ans[j],v[k],now,qpow(dsr[i].X,dsr[i].Y,1e12)));
ans.swap(nw),now*=qpow(dsr[i].X,dsr[i].Y,1e12);
}
sort(ans.begin(),ans.end());
cout<<ans.size()<<"
";
for(int i=0;i<ans.size();i++)printf("%lld%c",ans[i],"
"[i+1==ans.size()]);
}
signed main(){
// freopen("P5668_1.in","r",stdin);freopen("out.out","w",stdout);
int t;
cin>>t;
while(t--)mian();
return 0;
}