【51nod】1123 X^A Mod B （任意模數的K次剩餘）

題解

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();
    }
}