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

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

Description

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

Solution

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

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

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

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

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

3. 当 $1<(a,p^k)<a$ 时

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

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