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

时间 2019-11-24

标签 51nod mod 任意模数剩余繁體版

原文原文链接

题解

K次剩余终极版！orz
写一下，WA一年，bug不花一分钱node

在好久之前，我还认为，数论是一个重在思惟，代码很短的东西
后来。。。我学了BSGS，学了EXBSGS，学了模质数的K次剩余……代码一个比一个长……
直到今天，我写了240行的数论代码，我才发现数论这个东西= =太可怕了ios

好吧那么咱们来讲一下任意模数的K次剩余怎么搞
首先，若是模数是奇数，咱们能够拆成不少个质数的指数幂，再用解同余方程的方法一个个合起来，细节以后探讨ui

可是若是，模数有偶数呢
由于要输出全部解，解的个数很少，咱们能够倍增，也就是若是模数有一个\(2^k\)因子，那么它在模\(2^{k - 1}\)状况下的全部解x，若是还要在\(2^{k}\)成立，一定是\(x\)或\(x + 2^{k - 1}\)
咱们对于每个检查一下就行了spa

这样，咱们就只要考虑模奇数的状况了code

对于一个质数的指数幂\(p^{k}\) 有\(x^{A} \equiv C \pmod {p^{k}}\)
若\(C == 0\)
那么\(x\)中必然有\(p^{\lceil \frac{k}{A} \rceil}\)这个因子，以后从0枚举，一个个乘起来就是\(x\)可能的取值string

若\(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))}\)的时候）it

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

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

（附赠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 512022stream

代码

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