【简单数论题】【BZOJ2219】数论之神

题目地址：https://www.lydsy.com/JudgeOnline/problem.php?id=2219

Description

• 求 $x^N\equiv A(mod~P)$xN≡A(mod P) 在模 $P$P 意义下的解的数量，$P$P 是奇数。
• $N,A,P\le 10^9$N,A,P≤109。

Solution

• 神仙数论题。
• 我们知道，直接对这个方程求解有一定的困难。
• 我们设 $P=\prod_{i=1}^cp_i^{k_i}$P=∏i=1c​piki​​，$p_i$pi​ 是质数。
• 根据中国剩余定理，原方程解的数量即所有方程 $x^N\equiv A~\%~p_i^{k_i}(mod~p_i^{k_i})$xN≡A % piki​​(mod piki​​) 在模 $p_i^{k_i}$piki​​ 意义下的解的数量的乘积。
• 现在的问题即求形如 $x^N\equiv a(mod~p^k)$xN≡a(mod pk) 的方程在模 $p^k$pk 意义下的解的数量。
• 我们分类讨论一下。

1. 当 $(a,p^k)=a$(a,pk)=a 时

• 可以知道原方程等价于 $p^k|x^N$pk∣xN。
• 即 $x^N$xN 至少含有 $k$k 个质因子 $p$p。
• 所以 $p^{\lceil\frac{k}{N}\rceil}|x$p⌈Nk​⌉∣x。
• 所以在 $[0,p^k)$[0,pk) 中，$x$x 共有 $\frac{p^k}{p^{\lceil\frac{k}{N}\rceil}}=p^{k-\lceil\frac{k}{N}\rceil}$p⌈Nk​⌉pk​=pk−⌈Nk​⌉ 个。
• 该方程解即有 $p^{k-\lceil\frac{k}{N}\rceil}$pk−⌈Nk​⌉ 个。

2. 当 $(a,p^k)=1$(a,pk)=1 时

• 可以知道 $x\perp p$x⊥p 且 $x\in [0,p^k)$x∈[0,pk)。
• 我们找到 $p^k$pk 的一个原根 $g$g。
• 所以存在 $y=ind_gx$y=indg​x，$q=ind_ga$q=indg​a。
• 原方程等价于 $yN\equiv q(mod~\varphi(p^k))$yN≡q(mod φ(pk))。
• 这是一个线性同余方程。
• 当 $(N,\varphi(p^k))\nmid q$(N,φ(pk))∤q 时，原方程无解。
• 否则根据线性同余方程的性质，这个方程在 $y\in[0,\varphi(p^k))$y∈[0,φ(pk)) 的条件下，就有 $(N,\varphi(p^k))$(N,φ(pk)) 个解。
• 我们简单证明一下这个很常用的性质：
• 即证 $ax\equiv b(mod~p)$ax≡b(mod p) 若存在解，则有 $(a,p)$(a,p) 个模 $p$p 意义下的解。
• 证明：
  • 我们假设得到一个特解 $x_0$x0​，通解表示为 $x_0+k·\frac{p}{(a,p)}(k\in \mathbb Z)$x0​+k⋅(a,p)p​(k∈Z)。
  • 我们定义函数 $f(k)=x_0+\frac{p}{(a,p)}·k(k\in \mathbb Z)$f(k)=x0​+(a,p)p​⋅k(k∈Z)。
  • 即证 $f(k)\in [0,p)$f(k)∈[0,p) 恰有 $(a,p)$(a,p) 个整数解。
  • 显然 $f(k)$f(k) 单调递增。
  • 我们找到一个 $k_0$k0​ 使得 $f(k_0)\in[0,\frac{p}{(a,p)})$f(k0​)∈[0,(a,p)p​)。
  • 那么有 $f(k_0+(a,p))\in[p,p+\frac{p}{(a,p)})$f(k0​+(a,p))∈[p,p+(a,p)p​)，且 $f(k_0+(a,p)-1)\in[p-\frac{p}{(a,p)},p)$f(k0​+(a,p)−1)∈[p−(a,p)p​,p)。
  • 所以当 $k\in[k_0,k_0+(a,p))$k∈[k0​,k0​+(a,p)) 时，满足 $f(k)\in [0,p)$f(k)∈[0,p)，即不等式有 $(a,p)$(a,p) 个整数解。
  • 证毕。

3. 当 $1&lt;(a,p^k)&lt;a$1<(a,pk)<a 时

• 原方程可以写成 $x^N\equiv a&#x27;p^t(mod~p^k)$xN≡a′pt(mod pk)，$(a&#x27;,p)=1$(a′,p)=1。
• 当且仅当 $N|t$N∣t，原方程有解。
• 写成不定方程的形式，即 $x^N+p^kz=a&#x27;p^t$xN+pkz=a′pt。
• 两边同除 $p^t$pt，得 $(\frac{x}{p^{\frac{t}{N}}})^N+p^{k-t}z=a&#x27;$(pNt​x​)N+pk−tz=a′。
• 即 $(\frac{x}{p^{\frac{t}{N}}})^N\equiv a&#x27;(mod~p^{k-t})$(pNt​x​)N≡a′(mod pk−t)
• 即 $x&#x27;^N\equiv a&#x27;(mod~p^{k-t})$x′N≡a′(mod pk−t)。
• 对这个方程按照 2. 的方式求解即可得到 $[0,p^{k-t})$[0,pk−t) 内解的个数。
• 但实际上 $\frac{x}{p^{\frac{t}{N}}}\in [0,p^{k-\frac{t}{N}})$pNt​x​∈[0,pk−Nt​)。
• 因为我们解出来的 $x$x 都是模 $p^{k-t}$pk−t 的剩余系。
• 所以 $[0,p^{k-t})$[0,pk−t) 内的每一个解都对应 $[0,p^{k-\frac{t}{N}})$[0,pk−Nt​) 内的 $\frac{p^{k-\frac{t}{N}}}{p^{k-t}}=p^{{t-\frac{t}{N}}}$pk−tpk−Nt​​=pt−Nt​ 个解。
• 所以最后的答案还要乘以 $p^{{t-\frac{t}{N}}}$pt−Nt​。

#include <bits/stdc++.h>
 
const int MaxCnt = 1e5 + 5; 
 
int N, A, P; 
int cnt, p[MaxCnt], k[MaxCnt], pw[MaxCnt]; 
int a[MaxCnt], t[MaxCnt], tpw[MaxCnt]; 

inline int qpow(int a, int b)
{
	int res = 1; 
	for (; b; b >>= 1, a = a * a)
		if (b & 1)
			res = res * a; 
	return res; 
}

inline int qpow(int a, int b, const int &mod)
{
	int res = 1; 
	for (; b; b >>= 1, a = 1LL * a * a % mod)
		if (b & 1)
			res = 1LL * res * a % mod; 
	return res; 
}
 
inline int getphi(const int &p, const int &pw)
{
	return pw / p * (p - 1); 
}
 
inline int find_root(const int &mod, const int &phi)
{
	static int p[MaxCnt]; 
	static int cnt; 
	cnt = 0; 
	 
	int x = phi; 
	for (int i = 2; i * i <= phi; ++i)
		if (x % i == 0)
		{
			p[++cnt] = i; 
			while (x % i == 0)
				x /= i; 
		}
	if (x > 1)
		p[++cnt] = x; 
	 
	for (int g = 1; g <= mod; ++g)
	{
		if (std::__gcd(g, mod) > 1) continue; 
		bool flg = true; 
		for (int i = 1; i <= cnt; ++i)
			if (qpow(g, phi / p[i], mod) == 1)
			{
				flg = false; 
				break; 
			}
		if (flg)
			return g; 
	}
	return -1; 
}
 
inline int solve_BSGS(const int &a, const int &b, const int &mod)
{
	std::map<int, int> ha; 
	ha.clear(); 
	 
	int t = ceil(sqrt(mod)); 
	 
	int x = 1; 
	for (int i = 1; i <= t; ++i)
	{
		x = 1LL * x * a % mod; 
		ha[1LL * x * b % mod] = i; 
	}
	 
	int y = x; 
	for (int i = 1; i <= t; ++i)
	{
		if (ha.find(y) != ha.end())
			return i * t - ha[y]; 
		y = 1LL * y * x % mod; 
	}
	return -1; 
}
 
inline int exgcd(const int &a, const int &b, int &x, int &y)
{
	if (!b)
		return x = 1, y = 0, a; 
	int res = exgcd(b, a % b, y, x); 
	y -= a / b * x; 
	return res; 
}
 
inline int solve(const int &N, const int &a, const int &p_0, const int &p)
{
	int phi = getphi(p_0, p); 
	int g = find_root(p, phi); 
	int q = solve_BSGS(g, a, p); 
	 
	int d = std::__gcd(N, phi); 
	
	if (q % d)
		return 0; 
	 
	return d; 
}
 
int main()
{
	int T; 
	scanf("%d", &T); 
	while (T--)
	{
		scanf("%d%d%d", &N, &A, &P); 
		P = P << 1 | 1; 
		 
		cnt = 0; 
		int x = P; 
		for (int i = 2; i * i <= P; ++i)
			if (x % i == 0)
			{
				++cnt; 
				p[cnt] = i; 
				k[cnt] = 0; 
				pw[cnt] = 1; 
				while (x % i == 0)
					x /= i, ++k[cnt], pw[cnt] *= i; 
			}
		if (x > 1)
		{
			++cnt; 
			p[cnt] = pw[cnt] = x; 
			k[cnt] = 1; 
		}
		
		int ans = 1; 
		 
		bool flg = false; 
		for (int i = 1; i <= cnt; ++i)
		{
			a[i] = A % pw[i]; 
			t[i] = 0; 
			tpw[i] = 1; 
			while (a[i] && a[i] % p[i] == 0)
				a[i] /= p[i], ++t[i], tpw[i] *= p[i]; 
		 
			if (t[i] % N)
			{
				flg = true; 
				break; 
			}
			
			k[i] -= t[i]; 
			pw[i] /= tpw[i]; 
		}
		if (flg)
		{
			puts("0"); 
			continue; 
		}
		 
		for (int i = 1; i <= cnt; ++i)
		{
			if (!a[i])
				ans *= qpow(p[i], k[i] - ceil((double)k[i] / N)); 
			else if (!t[i])
				ans *= solve(N, a[i], p[i], pw[i]); 
			else
				ans *= solve(N, a[i], p[i], pw[i]) * qpow(p[i], t[i] - t[i] / N); 
			if (!ans)
				break; 
		}
		printf("%d\n", ans); 
	}
	return 0; 
}