題解
K次剩余終極版！orz
寫一下，WA一年，bug不花一分錢
在很久以前，我還認為，數論是一個重在思維，代碼很短的東西
后來。。。我學了BSGS，學了EXBSGS，學了模質數的K次剩余……代碼一個比一個長……
直到今天，我寫了240行的數論代碼，我才發現數論這個東西= =太可怕了
好吧那么我們來說一下任意模數的K次剩余怎么搞
首先，如果模數是奇數，我們可以拆成很多個質數的指數冪，再用解同余方程的方法一個個合起來，細節之后探討
但是如果，模數有偶數呢
因為要輸出所有解，解的個數不多，我們可以倍增，也就是如果模數有一個\(2^k\)因子，那么它在模\(2^{k - 1}\)情況下的所有解x，如果還要在\(2^{k}\)成立，必定是\(x\)或\(x + 2^{k - 1}\)
我們對於每一個檢查一下就好了
這樣，我們就只要考慮模奇數的情況了
對於一個質數的指數冪\(p^{k}\) 有\(x^{A} \equiv C \pmod {p^{k}}\)
若\(C == 0\)
那么\(x\)中必然有\(p^{\lceil \frac{k}{A} \rceil}\)這個因子，之后從0枚舉，一個個乘起來就是\(x\)可能的取值
若\(C \% p == 0\)
也就是說，\(C\)可以寫成\(u * p^{e}\)的形式，有解必定有\(A|e\)
那么就是\(x^{A} \equiv u * p^{e} \pmod {p^{k}}\)
把同余式打開，可以有\(x^{A} = u * p^{e} + h * p^{k}\)
等式兩邊都除掉一個\(p^{e}\)就有
\((\frac{x}{p^{\frac{e}{A}}})^{A} = u + h * p^{k - e}\)
設\(t = \frac{x}{p^{\frac{e}{A}}}\)
我們要求的就是
\(t^{A} \equiv u \pmod {p^{k - e}}\)
這時候\(u\)必然和模數互質，可以套用模質數的K次剩余
此時求出來的指標要取模的數是\(\phi(p^{k - e})\)而不是\(\phi(p^k)\)
之后求出所有指標的上界是\(\phi(p^k)\) (就是不斷增加\(\frac{\phi(p^{k - e})}{gcd(\phi(p^{k - e},A))}\)的時候）
如果\(C\)和\(p\)互質
那么直接上模質數的K次剩余（雖然是質數的指數冪但是你不需要用到有質數的那些位置了）
最后求完了，和上一次的答案用同余方程合起來即可
（附贈51nod上tangjz的hack數據，我雖然ac了然而這個hack沒過，又調了一段時間才過掉）
輸入
1
3 531441 330750
輸出
264 19947 39630 59313 78996 98679 118362 138045 157728 177411 197094 216777 236460 256143 275826 295509 315192 334875 354558 374241 393924 413607 433290 452973 472656 492339 512022
代碼
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <ctime>
//#define ivorysi
#define MAXN 100005
#define eps 1e-7
#define mo 974711
using namespace std;
typedef long long int64;
typedef unsigned int u32;
typedef double db;
int64 A,B,C,G,eu;
int64 ans[MAXN],tmp[MAXN],R,L[MAXN],cntL;
int M;
struct node {
int next,num;
int64 hsh;
}E[MAXN];
int head[mo + 5],sumE;
int64 fpow(int64 x,int64 c,int64 MOD) {
int64 res = 1,t = x;
while(c) {
if(c & 1) res = res * t % MOD;
t = t * t % MOD;
c >>= 1;
}
return res;
}
int primitive_root(int64 P,int64 eu) {
static int64 factor[1005];
int cnt = 0;
int64 x = eu;
for(int64 i = 2 ; i <= x / i ; ++i) {
if(x % i == 0) {
factor[++cnt] = i;
while(x % i == 0) x /= i;
}
}
if(x > 1) factor[++cnt] = x;
for(int G = 2 ; ; ++G) {
for(int j = 1 ; j <= cnt ; ++j) {
if(fpow(G,eu / factor[j],P) == 1) goto fail;
}
return G;
fail:;
}
}
int64 gcd(int64 a,int64 b) {
return b == 0 ? a : gcd(b,a % b);
}
void ex_gcd(int64 a,int64 b,int64 &x,int64 &y) {
if(b == 0) {
x = 1,y = 0;
}
else {
ex_gcd(b,a % b,y,x);
y -= a / b * x;
}
}
int64 Inv(int64 num,int64 MOD) {
int64 x,y;
ex_gcd(num,MOD,x,y);
x %= MOD;x += MOD;
return x % MOD;
}
void add(int u,int64 val,int num) {
E[++sumE].hsh = val;
E[sumE].next = head[u];
E[sumE].num = num;
head[u] = sumE;
}
void Insert(int64 val,int num) {
int u = val % mo;
for(int i = head[u] ; i ; i = E[i].next) {
if(val == E[i].hsh) {
E[i].num = num;
return;
}
}
add(u,val,num);
}
int Query(int64 val) {
int u = val % mo;
for(int i = head[u] ; i ; i = E[i].next) {
if(val == E[i].hsh) {
return E[i].num;
}
}
return -1;
}
int BSGS(int64 A,int64 C,int64 P) {
memset(head,0,sizeof(head));sumE = 0;
int64 S = sqrt(P);
int64 t = 1,invt = 1,invA = Inv(A,P);
for(int i = 0 ; i < S ; ++i) {
if(t == C) return i;
Insert(invt * C % P,i);
t = t * A % P;
invt = invt * invA % P;
}
int64 tmp = t;
for(int i = 1 ; i * S < P ; ++i) {
int x = Query(tmp);
if(x != -1) {
return i * S + x;
}
tmp = tmp * t % P;
}
}
bool Process(int64 A,int64 C,int64 P,int k) {
int64 MOD = 1,g;
for(int i = 1 ; i <= k ; ++i) MOD *= P;
cntL = 0;
if(C % MOD == 0) {
int64 T = (k - 1) / A + 1;
L[++cntL] = 0;
if(T < k) {
int64 num = fpow(P,T,MOD);
for(int i = 1 ; i * num < MOD ; ++i) L[++cntL] = i * num;
}
}
else if(g = gcd(C % MOD,MOD) != 1){
int64 x = C % MOD;
int c = 0;
while(x % P == 0) ++c,x /= P;
if(c % A != 0) return false;
G = primitive_root(MOD / (C / x),eu / (C / x));
eu /= C / x;
int e = BSGS(G,x,MOD / (C / x));
g = gcd(A,eu);
if(e % g != 0) return false;
e /= g;
int64 s = Inv(A / g,eu / g) * e % (eu / g);
L[++cntL] = s;
while(1) {
if((L[cntL] + eu / g) % (eu * (C / x)) == L[1]) break;
L[cntL + 1] = L[cntL] + eu / g;
++cntL;
}
for(int i = 1 ; i <= cntL ; ++i) {
L[i] = fpow(G,L[i],MOD) * fpow(P,c / A,MOD) % MOD;
}
}
else {
int e = BSGS(G,C % MOD,MOD);
g = gcd(A,eu);
if(e % g != 0) return false;e /= g;
int s = Inv(A / g,eu / g) * e % (eu / g);
L[++cntL] = s;
while(1) {
if(L[cntL] + eu / g >= eu) break;
L[cntL + 1] = L[cntL] + eu / g;
++cntL;
}
for(int i = 1 ; i <= cntL ; ++i) L[i] = fpow(G,L[i],MOD);
}
if(!cntL) return false;
if(!M) {
M = cntL;
for(int i = 1 ; i <= M ; ++i) ans[i] = L[i];
sort(ans + 1,ans + M + 1);
M = unique(ans + 1,ans + M + 1) - ans - 1;
R = MOD;
return true;
}
int tot = 0;
for(int i = 1 ; i <= M ; ++i) {
for(int j = 1 ; j <= cntL ; ++j) {
tmp[++tot] = (R * Inv(R,MOD) % (R * MOD) * (L[j] - ans[i]) + ans[i]) % (R * MOD);
tmp[tot] = (tmp[tot] + R * MOD) % (R * MOD);
}
}
R *= MOD;
sort(tmp + 1,tmp + tot + 1);
tot = unique(tmp + 1,tmp + tot + 1) - tmp - 1;
for(int i = 1 ; i <= tot ; ++i) ans[i] = tmp[i];
M = tot;
return true;
}
void Solve() {
M = 0;
if(B % 2 == 0) {
int64 Now = 2;B /= 2;
if(C & 1) ans[++M] = 1;
else ans[++M] = 0;
while(B % 2 == 0) {
B /= 2;
Now *= 2;
int t = 0;
for(int i = 1 ; i <= M ;++i) {
if(fpow(ans[i],A,Now) == C % Now) tmp[++t] = ans[i];
if(fpow(ans[i] + Now / 2,A,Now) == C % Now) tmp[++t] = ans[i] + Now / 2;
}
for(int i = 1 ; i <= t ; ++i) ans[i] = tmp[i];
if(!t) goto fail;
M = t;
}
R = Now;
sort(ans + 1,ans + M + 1);
M = unique(ans + 1,ans + M + 1) - ans - 1;
}
for(int64 i = 3 ; i <= B / i ; ++i) {
if(B % i == 0) {
eu = (i - 1);
B /= i;
int num = i,cnt = 1;
while(B % i == 0) {
B /= i;eu *= i;num *= i;++cnt;
}
G = primitive_root(num,eu);
if(!Process(A,C,i,cnt)) goto fail;
}
}
if(B > 1) {
eu = B - 1;
G = primitive_root(B,eu);
if(!Process(A,C,B,1)) goto fail;
}
if(M == 0) goto fail;
sort(ans + 1,ans + M + 1);
for(int i = 1 ; i <= M ; ++i) {
printf("%d%c",ans[i]," \n"[i == M]);
}
return;
fail:
puts("No Solution");
}
int main() {
#ifdef ivorysi
freopen("f1.in","r",stdin);
#endif
int T;
scanf("%d",&T);
while(T--) {
scanf("%lld%lld%lld",&A,&B,&C);
Solve();
}
}