BZOJ2219 数论之神 题解

\(\text{BZOJ2219 数论之神 题解}\)

记模数 \(P=2\times K+1\)。我们希望像上一个题一样转化为指标来求解，然而这里并没有保证 \(P\) 为质数，换句话说 \(P\) 甚至不一定有原根。于是考虑对 \(P\) 进行质因数分解，将 \(P\) 分解为 \(P=\prod p_i^{a_i}\) 的形式，这样一来对于 \(P\) 的同余转化为了一系列形如 \(x^A\equiv b_i\pmod {p_i^{a_i}}\) 的形式。考虑到这些个模数是互质的，由中国剩余定理可以知道在 \([0,P)\) 区间内只有唯一解，于是这些个线性同余方程组是互不相关的，于是求出每一个方程组在 \([0,p_i^{a_i})\) 范围内解的个数，答案就是其乘积。

现在考虑计算单个方程组解的个数。对于方程 \(x^A\equiv b_i\pmod {p_i^{a_i}}\)：

• 当 \(\gcd(b_i,p_i^{a_i})=1\) 时，形式和上一题是相同的。
• 当 \(\gcd(b_i,p_i^{a_i})=p_i^{a_i}\) 时：

显然可以转化为 \(x^A\equiv 0\pmod {p_i^{a_i}}\)。于是设 \(x=a\times p_i^w(\gcd(a,p_i)=1)\)，那么 \(a^Ap^{wA}\equiv 0\pmod{p_i^{a_i}}\)。这样一来 \(wA\ge a_i\)，于是有 \(w\ge \lceil \dfrac {a_i}{A}\rceil\)。由于要求 \(\lceil \dfrac {a_i}{A}\rceil\mid x\)，因此 \(x\) 的取值事实上是 \(\dfrac{p^{a_i}}{\lceil \frac {a_i}{A}\rceil}\) 种。

• 当 \(\gcd(b_i,p_i^{a_i})=p_i^w\) 时：

考虑给整个式子同除以一个 \(p_i^w\)。于是转化为了 \(\dfrac{x^A}{p_i^w}\equiv \pmod {p_i^{a_i-w}}\)。转化一下事实上左边的形式是 \((\dfrac x{p_i^{\frac wA}})^A\)。这样一来和第一种形式是相通的。但需要留意的是这样求出的 \(x\) 的范围实际上是 \([0,p_i^{a_i-w+\frac wA})\)，于是答案 \(\times p_i^{w-\frac wA}\)。

代码写的有些复杂了。

#include <bits/stdc++.h>
#define int long long
#define N 500005
using namespace std;
int T;
int a, b, k;
int qpow(int x, int y, int p) {
	int ans = 1;
	while (y) {
		if (y & 1) ans = ans * x % p;
		x = x * x % p;
		y >>= 1;
	}
	return ans;
}
int BSGS(int a, int b, int p) {
	unordered_map<int, int>mp;
	a %= p, b %= p;
	int s = ceil(sqrt(p)), t = b;
	mp[t] = 0;
	for (int i = 1; i <= s; i++) {
		t = t * a % p;
		mp[t] = i;
	}
	int x = qpow(a, s, p);
	t = 1;
	for (int i = 1; i <= s; i++) {
		t = t * x % p;
		if (mp.count(t)) return i * s - mp[t];
	}
	return -1;
}
int p[N], t[N], cnt;
void divd(int x) {
	cnt = 0;
	for (int i = 2; i * i <= x; i++)
		if (x % i == 0) {
			p[++cnt] = i;
			while (x % i == 0)
				x /= i, ++t[cnt];
		}
	if (x > 1) p[++cnt] = x, t[cnt] = 1;
}
int phi(int x) {
	int ans = x;
	for (int i = 2; i * i <= x; i++)
		if (x % i == 0) {
			ans = ans / i * (i - 1);
			while (x % i == 0) x /= i;
		}
	if (x > 1) ans = ans / x * (x - 1);
	return ans;
}
int chk(int g, int m, int ph) {
	if (qpow(g, ph, m) != 1) return 0;
	for (int i = 1; i <= cnt; i++)
		if (qpow(g, ph / p[i], m) == 1) return 0;
	return 1;
}
int fnd(int m, int ph) {
	for (int i = 1; i < m; i++)
		if (chk(i, m, ph)) return i;
	return 0;
}
int sex(int a, int b, int p) {
	int g = __gcd(a, p);
	if (b % g) return 0;
	return g;
}

int fp[N], ft[N], fcnt;
void solve() {
	memset(p, 0, sizeof p);
	memset(t, 0, sizeof t);
	memset(fp, 0, sizeof fp);
	memset(ft, 0, sizeof ft);
	cnt = fcnt = 0;
	cnt = fcnt = 0;
	int a, b, k, P, ph;
	cin >> a >> b >> k;
	P = 2 * k + 1;
	divd(P);
	int ans = 1;
	for (int i = 1; i <= cnt; i++) fp[i] = p[i], ft[i] = t[i];
	fcnt = cnt;
	for (int i = 1; i <= cnt; i++) {
		int np = qpow(p[i], t[i], P + 1);
		int B = b % np;
		int gc = __gcd(B, np);
		if (gc == 1) {
			ph = phi(np);
			divd(ph);
			int g = fnd(np, ph);
			int t = BSGS(g, B, np);
			ans = ans * sex(a, t, ph);
		}
		else if (gc == np) ans = ans * qpow(p[i], t[i] - (int)ceil(t[i] * 1.0 / a), P + 1);
		else {
			int w = 0;
			for (int j = 1; j <= 50; j++)
				if (qpow(p[i], j, P + 1) == gc) w = j;
			np /= qpow(p[i], w, P + 1);
			ph = phi(np);
			B /= qpow(p[i], w, P + 1);
			int bp = p[i];
			divd(ph);
			int g = fnd(np, ph);
			int t = BSGS(g, B, np);
			ans = ans * sex(a, t, ph);
			ans = ans * qpow(bp, w - (int)ceil(w * 1.0 / a), P + 1);
		}
		cnt = fcnt;
		for (int j = 1; j <= cnt; j++) p[j] = fp[j], t[j] = ft[j];
	}
	cout << ans << '\n';
}

signed main() {
	ios::sync_with_stdio(0);
	cin.tie(0);
	cin >> T;
	while (T--) solve();
	return 0;
}