题解
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();
}
}原文地址：https://www.cnblogs.com/ivorysi/p/9050357.html
时间： 2024-10-07 21:32:21