BZOJ2219 数论之神

题意：求 $x^a\equiv b(\mathrm{mod}p)$ 的解数，$a,b,p$ 给定。$a,b,p\le {10}^9+1$，$p$ 是奇数，$T(T\le 1000)$ 组数据。

容易发现单纯的 $p$ 不好做，于是显然根据中国剩余定理拆成 $p^k(p\in P)$ 做，然后直接把各部分的答案乘起来即可。

下面直接讨论 $x^a\equiv b(\mathrm{mod}{p}^k)$ 的解数。

Case 1：$b=0$

容易发现只要 $x$ 在 $p$ 的次幂至少是 $\lceil {\displaystyle \frac{k}{a}}\rceil$ 即可。于是此时的贡献就是 $p^{k-\lceil \frac{k}{a}\rceil }$。

Case 2：$bmodp=0\wedge b\ne 0$

此时假设 $b=c\times p^t$，那么此时 $x$ 在 $p$ 的次幂一定是 ${\displaystyle \frac{t}{a}}$（注意一定整除，否则无解），然后可以对两边同时除以 $p^t$，那么方程就变为了 $x^a\equiv c(\mathrm{mod}{p}^{k-t})$。注意到 $mod{p}^{k-t}$ 映射到 $mod{p}^k$ 时，由于需要乘上 $p^{\frac{t}{a}}$，所以剩下的 $k-(k-t)-{\displaystyle \frac{t}{a}}=t-{\displaystyle \frac{t}{a}}$ 次幂是完全定义域扩大的，所以需要乘上 $p^{t-\frac{t}{a}}$。剩下的是第三种情况，就可以解决这个问题了。

Case 3：$bmodp\ne 0$

这个可以使用离散对数（知道奇数的作用了吗？），假设 $b=g^t$，那么就是 $x^a\equiv g^t(\mathrm{mod}{p}^k)$，容易发现如果 $tmodgcd(\phi (p^k),a)\ne 0$，说明不可能凑出 $t$，这样的话就返回 0。

然后考虑有解的情况有多少个。假设 $x=g^y$，那么就是 $ay\equiv t(\mathrm{mod}\phi (p^k))$。容易发现此时当找到一个解时，$+{\displaystyle \frac{\text{lcm}(\phi (p^k),a)}{a}}$ 也是一个解，而且是最小周期。除一下，一共就有 $gcd(\phi (p^k),a)$ 个解。

把所有的乘起来即可。时间复杂度 $O(\sqrt{p}\log p)$ 或 $O(\sqrt{p})$，瓶颈在求离散对数的 BSGS 上。

int BSGS(int a, int b, int Mod)
{
	if (b == 1 || Mod == 1) return 0;
	std::unordered_map<int, int> H;
	int cur = b, K = std::sqrt(Mod) + 1;
	for (int i = 0; i < K; ++ i, cur = (LL) cur * a % Mod) H[cur] = i;
	int ak = cur = qpow(a, K, Mod);
	for (int i = 1; i <= K; ++ i, cur = (LL) cur * ak % Mod)
		if (H.count(cur)) return i * K - H[cur];
	return -1;
}

int findrt(int p, int k)
{
	int pk = qpow(p, k);
	std::vector<int> fac;
	int n = pk - 1;
	for (int i = 2; i <= n / i; ++ i)
		if (n % i == 0) {
			fac.push_back(i);
			while (n % i == 0) n /= i;
		}
	if (n ^ 1) fac.push_back(n);
	auto check = [&](int g) {
		for (int x : fac)
			if (qpow(g, pk / p * (p - 1) / x, pk) == 1) return false;
		return true;
	};
	for (int i = 2; i < pk; ++ i)
		if (check(i)) return i;
	return -1;
}

int solve(int a, int b, int p, int k)
{
	int pk = qpow(p, k);
	b %= pk;
	// std::cout << a << ' ' << b << ' ' << p << ' ' << k << '\n';
	if (b == 0) return qpow(p, k - (k + a - 1) / a);
	if (b % p == 0) {
		int t = 0;
		while (b % p == 0) b /= p, ++ t;
		if (t % a) return 0;
		return solve(a, b, p, k - t) * qpow(p, t - t / a);
	}
	int g = findrt(p, k), _b = BSGS(g, b, pk);
	// std::cout << pk << ' ' << g << ' ' << b << ' ' << _b << ' ' << Gcd(a, pk / p * (p - 1)) << std::endl;
	if (_b % Gcd(a, pk / p * (p - 1))) return 0;
	return Gcd(a, pk / p * (p - 1));
}

void work()
{
	int a, b, p, res = 1;
	scanf("%d %d %d", &a, &b, &p), p = 2 * p + 1;
	for (int i = 2; i <= p / i; ++ i)
		if (p % i == 0) {
			int t = 0;
			while (p % i == 0) ++ t, p /= i;
			// std::cout << "Solve " << a << ' ' << b << ' ' << i << ' ' << t << '\n';
			res *= solve(a, b, i, t);
		}
	if (p ^ 1) res *= solve(a, b, p, 1);
	printf("%d\n", res);
}