高次剩餘學習筆記 / 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^a\equiv b\pmod{p^\alpha}(p\in\mathbb P)\) 怎麼解。注意到當 \(p\neq 2\) 時 \(p^{\alpha}\)
如果 \(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\)，那麼 \(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(\sqrt p)\)，而 \(|X|\) 上限是有保證的）。直接暴力列舉 \([0,p)\) 內所有 \(x\) 顯然會 T，我們可以剪枝呀！考慮對 \(\alpha\) 進行迭代，假設 \(x^a\equiv b\pmod{p^{\alpha-1}}\) 的解集已經算出為 \(X_{\alpha-1}\)，那麼 \(x\) 是 \(x^a\equiv b\pmod{p^{\alpha}}\) 的解的必要條件是滿足 \(\alpha-1\) 的方程，即 \(x\in X_{\alpha-1}\pmod{p^{\alpha-1}}\)，而對於每個 \(x_0\in X_{\alpha-1}\)，在 \(\bmod p^\alpha\) 意義下有兩個對應值，都快速冪判一下即可。這樣最壞複雜度依然是跑滿的，但是剪枝效果很好（可以通過模板題），畢竟 \(p^\beta\) 時減掉一個枝相當於直接砍掉 \(p^{\alpha-\beta}\) 個可能的答案。
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<<"!!\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;
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()<<"\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;
}