51nod1123 X^A Mod B （任意模数的K次剩余）

Posted 2020-11-06 sigongzi

题解

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(k)\)
之后求出所有指标的上界是\(\phi(k)\) (就是不断增加\(\frac{\phi(p^{k - e})}{gcd(\phi(p^{k - e},A))}\)的时候）

如果\(C\)和\(p\)互质
那么直接上模质数的K次剩余（虽然是质数的指数幂但是你不需要用到有质数的那些位置了）

最后求完了，和上一次的答案用同余方程合起来即可

（附赠51nod上tangjz的hack数据，我虽然ac了然而这个hack没过，又调了一段时间才过掉）
1
3 531441 330750

代码

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