洛谷 P5345: 【XR-1】快乐肥宅

题目传送门:洛谷 P5345

很荣幸为 X Round 1 贡献了自己的一题。

题意简述:

给定 \(n\)\(k_i,g_i,r_i\)\(0\le k_i,r_i<g_i\le 10^7\))。

求关于 \(x\) 的方程组 \(\left\{\begin{matrix}k_1^x\equiv r_1\pmod{g_1}\\k_2^x\equiv r_2\pmod{g_2}\\\vdots\\k_n^x\equiv r_n\pmod{g_n}\end{matrix}\right.\)(定义 \(0^0=1\))在 \([0,10^9]\) 内的最小整数解,或判断在这个范围内无解。

题解:

分为两个部分解决。

第一部分:分别求出每个方程的解。

首先观察 \(k^j\bmod g\) 的规律,以 \(g=495616\)\(k=124\) 为例:

定义 \(f(i)=ki\bmod g\),则把每个在 \([0,g)\) 之间的整数看作一个节点后,形成了一个基环内向森林。

\(1\bmod g\) 开始走唯一的出边,必然形成一个 \(\rho\) 形结构,让我们模拟这个过程:

可以看到,出现了长度为 \(5\) 的循环。

  1. 如果 \(r=245760\),则 \(x\equiv 7\pmod{5}\),且 \(x\ge 7\)

  2. 如果 \(r=12544\),则 \(x=4\)

  3. 如果 \(r=2\),则无解。

以上就是每个方程的 \(3\) 种解,无限解,唯一解或无解。

那么,问题就是,如何区分这 \(3\) 种解,以及如何求出解。

首先,需要将 \(\rho\) 形的“尾巴”和循环分开考虑,不难证明“尾巴”的长度不超过 \(\log_2g\)

而且可以发现,若在“尾巴”上找到了解,那么就只有一组解。所以需要判断一个值是否在尾巴上。

其实这是很简单的,如果 \(\gcd(x,m)\neq\gcd(f(x),m)\),那么 \(x\) 就在尾巴上,反之就在循环内。

据此,我们区分了在“尾巴”上的解,也就是只有一组解的情况,接下来考虑不在“尾巴”上的情况。

因为解不在“尾巴”上,则如果解不在循环内就是无解了,所以首先判断解是否在循环内。

假设第一个在循环上的数为 \(c\),尾巴长度为 \(o\),则有 \(c\equiv k^o\pmod{g}\)

\(d=\gcd(c,g)\),则原方程 \(k^x\equiv r\pmod{g}\) 可化为 \(c\cdot k^{x-o}\equiv r\pmod{g}\)\(x\ge o\),这是因为我们假定解一定在循环内。

可以进一步化为 \(k^{x-o}\equiv\frac{r}{d}\left(\frac{c}{d}\right)^{-1}\pmod{\frac{g}{d}}\),若 \(r\) 不是 \(d\) 的倍数则无解。

这是因为循环内的值均是 \(d\) 的倍数,可以让方程同除 \(d\),又因为 \(\frac{c}{d}\perp\frac{g}{d}\),可以取逆元。

\(a=k,b=\frac{r}{d}\left(\frac{c}{d}\right)^{-1},m=\frac{g}{d}\),则有 \(a\perp m\)

使用 BSGS 求出 \(a^y\equiv b\pmod{m}\) 的最小自然数解 \(y\) 后,有原方程的最小自然数解为 \(y+o\)

再使用 BSGS 求出 \(a^z\equiv 1\pmod{m}\) 的最小正整数解 \(z\),即 \(a\)\(m\) 的阶。

则原方程的解为 \(x\equiv y+o\pmod{z}\)\(x\ge y+o\)

至此,第一部分解决。

第二部分:合并每个方程的解。

首先,若前面的方程出现了无解的情况,则方程组也无解。

若前面的方程出现了只有唯一解的情况,只需要检查此唯一解是否满足所有方程即可。

接下来讨论以上每个方程的解均形如 \(x\equiv r\pmod{q}\)\(x\ge r\) 的形式。

分开考虑前半部分和后半部分,对于前半部分,显然是 ExCRT 的形式。

对于后半部分,可以化为 \(x\ge\max r_i\),令 \(\max r_i=x_0\),放到最后考虑。

考虑使用 ExCRT 解决前半部分,令前 \(i\) 个方程组的解为 \(x\equiv P_i\pmod{Q_i}\),有 \(P_0=0,Q_0=1\)

但是这里出现一个问题,那就是 \(P_i,Q_i\) 可能很大,但是这里其实没必要使用高精度,要怎么处理这种情况呢?

因为以上方程的解 \(x\equiv r\pmod{q}\) 中,\(q\) 均小于 \(10^7\),所以不需要担心这部分。

考虑这样一种情况:\(Q_{i-1}>10^9\),而再合并进 \(x\equiv r_i\pmod{q_i}\) 可能会让 \(Q_i\) 超出 long long 能够表示的范围。

这时其实不需要再合并了,只需要判断 \(P_{i-1}\) 是否满足第 \(i\) 个方程即可。因为若不满足,合并后 \(x\) 必然会至少增加一倍的 \(Q_{i-1}\),这就超出了 \(10^9\) 的范围,从而可以直接输出无解。

最后,当方程成功合并完之后,\(Q_n\) 可能不是真实值,但是当 \(Q_n\le 10^9\) 时必然是真实值。

此时可以再考虑 \(x_0\),真实的解应该为 \(x=\left\lceil\frac{x_0-P_n}{Q_n}\right\rceil\cdot Q_n+P_n\)

\(Q_n> 10^9\),则必然要满足 \(P_n\ge x_0\),否则无解。

据此写出代码,复杂度 \(\mathcal{O}\left(\sum_{i=1}^{n}\left(\sqrt{g_i}+\log g_i\right)\right)\)

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>

#define Fail return puts("Impossible"), 0
#define mp std::make_pair

typedef long long LL;
typedef std::pair<int, int> pii;
const int MG = 10000005, MS = 3175;
const int U = 1e9;

int gcd(int a, int b) { return b ? gcd(b, a % b) : a; }

template<typename T>
T exgcd(T a, T b, T &x, T &y) {
	if (!b) return x = 1, y = 0, a;
	int d = exgcd(b, a % b, y, x);
	return y -= a / b * x, d;
}

int buk[MG], stk[MS];
inline int BSGS(int a, int b, int m) {
	int S = sqrt(m - 1) + 1;
	int A = 1, f = -1, t = 0;
	for (int i = 0; i < S; ++i) {
		buk[stk[++t] = (LL)b * A % m] = i;
		A = (LL)A * a % m;
	}
	int C = 1;
	for (int i = 1; !~f &&  i <= S; ++i)
		if (~buk[C = (LL)C * A % m])
			f = i * S - buk[C];
	while (t) buk[stk[t--]] = -1;
	return f;
}

inline pii ExBSGS(int a, int b, int m) {
	int o = 0, A = 1 % m, d = 1, nd, x, y;
	while (1) {
		if (d == (nd = gcd((LL)A * a % m, m))) break;
		if (A == b) return mp(o, -1);
		++o, A = (LL)A * a % m, d = nd;
	}
	if (b % d) return mp(-1, -1);
	m /= d, b /= d, A /= d;
	exgcd(A, m, x, y);
	b = (LL)b * (x + m) % m;
	x = BSGS(a, b, m);
	if (!~x) return mp(-1, -1);
	y = BSGS(a, 1 % m, m);
	return mp(x % y + o, y);
}

inline bool Combine(LL &a1, LL &m1, LL a2, LL m2) {
	LL k1, k2, g = exgcd(m1, m2, k1, k2);
	if ((a2 - a1) % g) return 0;
	a1 += (k1 * ((a2 - a1) / g) % m2 + m2) * m1 % (m1 / g * m2);
	return a1 %= m1 *= m2 / g, 1;
}

const int MN = 1005;

int N;
int k[MN], r[MN], g[MN];
int x[MN], q[MN], X, MaxX;

int main() {
	memset(buk, -1, sizeof buk), X = -1;
	scanf("%d", &N);
	for (int i = 1; i <= N; ++i) {
		scanf("%d%d%d", &k[i], &g[i], &r[i]);
		pii ans = ExBSGS(k[i] % g[i], r[i] % g[i], g[i]);
		x[i] = ans.first, q[i] = ans.second;
		if (!~x[i]) Fail;
		if (!~q[i]) X = x[i];
		MaxX = std::max(MaxX, x[i]);
	}
	if (~X) {
		for (int i = 1; i <= N; ++i) {
			if (~q[i] && (X < x[i] || (X - x[i]) % q[i])) Fail;
			if (!~q[i] && X != x[i]) Fail;
		}
		return printf("%d\n", X), 0;
	}
	LL P = 0, Q = 1;
	for (int i = 1; i <= N; ++i) {
		if (Q > U && (P - x[i]) % q[i]) Fail;
		if (Q <= U && !Combine(P, Q, x[i] % q[i], q[i])) Fail;
	}
	if (P < MaxX) P += ((MaxX - P - 1) / Q + 1) * Q;
	if (P > U) Fail;
	printf("%lld\n", P);
	return 0;
}
posted @ 2019-05-05 18:45  粉兔  阅读(1975)  评论(0编辑  收藏  举报