高次剩余学习笔记 / P5668 - 【模板】N次剩余 题解
考前学算法，就是这么浪。
求 \(x^a\equiv b\pmod p\) 在 \([0,p)\) 内的所有解。
核心思想是通过原根将幂变成指数。但是普通的模数 \(p\) 不一定有原根啊，那考虑使用分解质因数 \(p=\prod p_i^{\alpha_i}\)，奇质数的幂是有原根的。拆开的理论依据是什么？换个元，令 \(y=x^a\)，那么 \(y\equiv b\pmod p\) 可以由若干个 \(y\equiv b\pmod{p_i^{\alpha_i}}\) 联立而成的方程组等价表达，根据的是模数两两互质的 ExCRT（其实就是 CRT，可惜我不会写）。那么我们得到 \(y\equiv b\pmod{p_i^{\alpha_i}}\) 的解集 \(X_i\) 后，其实就是 \(x\) 是原方程的解当且仅当（方程和方程组的等价性已经用 CRT 证明）对所有 \(i\) 都有 \(x\in X_i\pmod{p_i^{\alpha_i}}\)，于是我们可以从每个 \(X_i\) 各挑一个值出来用 CRT 合并，得到 \(\prod|X_i|\) 个 \([0,p)\) 内的 \(x\) 便是原方程的所有解。这里使用了两次 CRT，正着一次反着一次，步步为等价变换。
下面讨论 \(x^a\equiv b\pmod{p^\alpha}(p\in\mathbb P)\) 怎么解。注意到当 \(p\neq 2\) 时 \(p^{\alpha}\) 必有原根，我们分成 \(p\overset{?}=2\) 两种情况讨论，先考虑 \(p\neq 2\)。首先如果 \(b=0\)，那么显然解集为 \(p^{\left\lceil\frac \alpha a\right\rceil}\mid x\)，下面 \(b\neq0\)。若 \(b\perp p^{\alpha}\)，求得 \(p^{\alpha}\) 的任意一个原根（可以用 1/4 次 log 的算法求最小原根）\(g\)，然后用 BSGS 求出 \(b=g^\beta\) 的 \(\left[0,\varphi=\varphi\!\left(p^\alpha\right)\right)\) 内唯一 \(\beta\)。此时显然 \(x\) 是解仅当 \(x\perp p^\alpha\)，于是令 \(x=g^y(y\in[0,\varphi))\)，求出 \(y\) 在 \(\bmod\varphi\) 意义下所有解，就相当于求出了 \(x\) 的所有解。原方程即 \(g^{ay}\equiv g^\beta\pmod{p^\alpha}\)，跟据原根的幂周期为 \(\varphi\)，这显然等价于 \(ay\equiv\beta\pmod{\varphi}\)。线性同余方程就没什么讲的必要了吧，exgcd 找特解然后枚举通解在 \([0,\varphi)\) 内所有取值即可。
如果 \(b\not\perp p^\alpha\) 呢？此时 \(b\) 在模意义下没有原根。设 \(b=p^\beta q\)，令 \(x=p^\gamma r\)，那么原方程等价于 \(p^{a\gamma}r^a\equiv p^\beta q\pmod{p^\alpha}\)。此时我们发现，由于 \(b\neq 0\)，也就是 \(\beta<+\infty\) 即 \(\beta<\alpha\)，那么 \(p^{\beta} q\) 所在剩余系中都恰好包含 \(\beta\) 个质因子 \(p\)，那么我们需要 \(a\gamma=\beta\)（并不是模意义下，而是严格相等！）！那么判一手 \(a\mid\beta\)，直接得到 \(\gamma=\dfrac\beta a\)，它现在是常数了。跟据同余式的性质将三边的 \(p^\beta\) 约掉，得到 \(r^a\equiv q\pmod{p^{\alpha-\beta}}\)，此时有 \(p\nmid q\)，于是归约到了 \(b\perp p^\alpha\) 的情况。求出 \(r\bmod p^{\alpha-\beta}\) 的所有可能值后，对于其中每一个，我们要找出对应的所有 \(\left[0,p^\alpha\right)\) 内的 \(p^\gamma r\)，那么找到所有 \(r\in\left[0,p^{\alpha-\gamma}\right)\) 即可，随便遍历一波，有 \(p^{(\alpha-\gamma)-(\alpha-\beta)}=p^{\beta-\gamma}\) 个值。
最后我们来考虑 \(p=2\) 的情形。由于没有原根，很可惜，这种情况只能暴力（有原根的情况下复杂度显然是 \(|X|+\mathrm O\!\left(\sqrt {p^\alpha}\right)\)，而 \(|X|\) 上限是有保证的）。直接暴力枚举 \([0,2^\alpha)\) 内所有 \(x\) 显然会 T，我们可以剪枝呀！考虑对 \(\alpha\) 进行迭代，假设 \(x^a\equiv b\pmod{2^{\alpha-1}}\) 的解集已经算出为 \(X_{\alpha-1}\)，那么 \(x\) 是 \(x^a\equiv b\pmod{2^{\alpha}}\) 的解的必要条件是满足 \(\alpha-1\) 的方程，即 \(x\in X_{\alpha-1}\pmod{2^{\alpha-1}}\)，而对于每个 \(x_0\in X_{\alpha-1}\)，在 \(\bmod 2^\alpha\) 意义下有两个对应值，都快速幂判一下即可。这样最坏复杂度依然是跑满的，但是剪枝效果很好（可以通过模板题），毕竟 \(2^\beta\) 时减掉一个枝相当于直接砍掉 \(2^{\alpha-\beta}\) 个可能的答案。
upd：被这个帖子 hack 掉了。还以为是什么大问题，没想到是个偷袭（雾）。设分解到的每个质因子的解数为 \(X_i\)，那么总共要做 \(\mathrm O(X=\prod X_i)\) 次 CRT，而 CRT 不会出现无法合并，也就是总解数也是 \(\mathrm O(X)\)，而题目保证 \(X\) 不会太大，所以复杂度是对的。但有一个漏洞：有 \(X_i=0\) 的话，「只要做 \(\mathrm O(X)\) 次 CRT」便不成立，这里在开始 CRT 之前特判一下就好了。
代码（已经换码风了，还要适应之前的码风好烦）：
#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<<"!!\n";
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<<":\n";
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;
vector<vector<int>> vv;
for(int i=0;i<dsr.size();i++){
vector<int> v=solve(a,b,dsr[i].X,dsr[i].Y);
vv.pb(v);
if(v.empty())return puts("0"),void();
}
for(int i=0;i<dsr.size();i++){
auto v=vv[i];
vector<int> 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()<<"\n";
for(int i=0;i<ans.size();i++)printf("%lld%c",ans[i]," \n"[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;
}
