230527 // 中国剩余定理

真 __,复 __ 周 __ !


给定下列关于 \(x\) 的一元同余方程组:

\[\begin {cases} x \equiv a_1 \pmod {m_1} \\ x \equiv a_2 \pmod {m_2} \\ \quad \quad \vdots \\ x \equiv a_k \pmod {m_k} \end {cases} \]

其中对于 \(\forall \, 1\le i, j\le k\, (i\ne j)\),满足 \(m_i \, \bot\, m_j\)。下面是求解该方程组的具体方法:

\(M = \prod \limits _{i = 1} ^ k m_i, \, M_i=\dfrac M {m_i}\),则 \(x = \sum \limits _{i = 1} ^k a_i\times M_i\times M_i ^{-1} \pmod M\)

其中,\(M_i^{-1}\)\(M_i\) 在模 \(m_i\) 意义下的逆元(所以 \(M_i\times M_i^{-1}\) 的值并不是视觉上的 \(-1\))。

正确性证明:令 \(A=\sum\limits_{i=1}^k a_i\dot M_i \dot {M_i}^{-1}\)。对于 \(\forall \, A'\in \Z\)\(A'\equiv A\pmod M\),下证 \(A'\) 是原同余方程的解:

对于任意 \(i\in \{1, 2, \cdots, k\}\)

没写完

时间复杂度 \(\mathcal O(n\log n)\)。注意求解逆元时要用到 exgcd。

inline int CRT(int k, int *a, int *m) {
	int M = 1, x = 0;
	static int M1[maxn], Mi[maxn]; 
	for (int i = 1; i <= k; ++i)
		M *= m[i];
	for (int i = 1; i <= k; ++i) {
		M1[i] = M / m[i];
		Mi[i] = getinv(M1[i], m[i]);
		(x += a[i] * M1[i] * Mi[i]) %= M;
	}
	return x;
}

A. 曹冲养猪

http://222.180.160.110:1024/contest/3642/problem/1

板。

#define int long long
namespace XSC062 {
using namespace fastIO;
const int maxn = 15;
int n;
int a[maxn], b[maxn];
int exgcd(int a, int &x, int b, int &y) {
	if (b == 0) {
		x = 1, y = 0;
		return a;
	}
	int u = exgcd(b, x, a % b, y);
	int t = x;
	x = y, y = t - (a / b) * y;
	return u;
}
inline int getinv(int a, int p) {
	int k1, k2;
	int g = exgcd(a, k1, p, k2);
	if (g != 1)
		return -1;
	return (k1 % p + p) % p;
}
inline int CRT(int k, int *a, int *m) {
	int M = 1, x = 0;
	static int M1[maxn], Mi[maxn]; 
	for (int i = 1; i <= k; ++i)
		M *= m[i];
	for (int i = 1; i <= k; ++i) {
		M1[i] = M / m[i];
		Mi[i] = getinv(M1[i], m[i]);
		(x += a[i] * M1[i] * Mi[i]) %= M;
	}
	return x;
}
int main() {
	read(n);
	for (int i = 1; i <= n; ++i)
		read(a[i]), read(b[i]);
	print(CRT(n, b, a));
	return 0;
}
} // namespace XSC062
#undef int

B. 阿九大战朱最学

http://222.180.160.110:1024/contest/3642/problem/2

它甚至连样例都不愿意改一下,就是上一道题啊。


C. 韩信点兵

http://222.180.160.110:1024/contest/3642/problem/3

还是板。

#define int long long
namespace XSC062 {
using namespace fastIO;
const int maxn = 15;
int n, res;
int a[maxn], b[maxn];
int exgcd(int a, int &x, int b, int &y) {
	if (b == 0) {
		x = 1, y = 0;
		return a;
	}
	int u = exgcd(b, x, a % b, y);
	int t = x;
	x = y, y = t - (a / b) * y;
	return u;
}
inline int getinv(int a, int p) {
	int k1, k2;
	int g = exgcd(a, k1, p, k2);
	if (g != 1)
		return -1;
	return (k1 % p + p) % p;
}
inline int CRT(int k, int *a, int *m) {
	int M = 1, x = 0;
	static int M1[maxn], Mi[maxn]; 
	for (int i = 1; i <= k; ++i)
		M *= m[i];
	for (int i = 1; i <= k; ++i) {
		M1[i] = M / m[i];
		Mi[i] = getinv(M1[i], m[i]);
		(x += a[i] * M1[i] * Mi[i]) %= M;
	}
	return x;
}
int main() {
	n = 3;
	a[1] = 3, a[2] = 5, a[3] = 7;
	read(b[1]), read(b[2]), read(b[3]);
	res = CRT(n, b, a);
	if (res < 10 || res > 100)
		puts("no answer");
	else print(res);
	return 0;
}
} // namespace XSC062
#undef int

扩展中国剩余定理

设有如下同余方程组:

\[\begin {cases} x \equiv a_1 \pmod {m_1} \\ x \equiv a_2 \pmod {m_2} \end {cases} \]

不保证 \(m_1 \, \bot \, m_2\)

\(x = m_1 \times p_1 + a_1 = m_2 \times p_2 + a_2\),其中 \(p_i \in \Z\),则有 \(p_1 \times m_1 - p_2 \times m_2 = a_2 - a_1\)

\(p_1\) 的值可以用 exgcd 求解,则原方程的解满足 \(x \equiv m_1\times p_1 + a_1 \pmod {M}\),其中 \(M=\operatorname{lcm}(m_1, m_2)\)

这样我们就得到了一个新的同余方程。对于 \(k>2\) 的情况,我们不断合并两个同余方程即可得到最终同余方程。

此时根据题目要求(一般求最小值即 \(m_k\times p_k + a_k\))求解答案即可。

时间复杂度 \(\mathcal O(n\log n)\),其中 \(\log\) 来自 exgcd。

inline int calc(int m1, int a1,
					int m2, int a2) {
	int p1, p2;
	if (a2 - a1 < 0)
		swap(a1, a2), swap(m1, m2);
	int g = exgcd(m1, p1, m2, p2);
	if ((a2 - a1) % g)
		return -1;
	p1 *= (a2 - a1) / g, m2 /= g;
	p1 = (p1 % m2 + m2) % m2;
	return p1 * m1 + a1;
}
inline int exCRT(int k, int *a, int *m) {
	for (int i = 2; i <= k; ++i) {
		a[i] = calc(m[i - 1],
				a[i - 1], m[i], a[i]);
		if (a[i] == -1)
			return -1;
		m[i] = lcm(m[i - 1], m[i]);
	}
	return a[k] % m[k];
}

D. Strange Way to Express Integers

http://222.180.160.110:1024/contest/3642/problem/4

板。但是要开 __int128。不知道为什么智力只用开 long long 就能跑过,我和揭哥就不行。

#define int __int128
namespace XSC062 {
using namespace fastIO;
const int maxn = 1e5 + 5;
int n;
int a[maxn], m[maxn];
int gcd(int x, int y) {
	return y ? gcd(y, x % y) : x;
}
inline int lcm(int x, int y) {
	return x / gcd(x, y) * y;
}
inline void swap(int &x, int &y) {
	x ^= y ^= x ^= y;
	return;
}
int exgcd(int a, int &x, int b, int &y) {
	if (b == 0) {
		x = 1, y = 0;
		return a;
	}
	int u = exgcd(b, x, a % b, y);
	int t = x;
	x = y, y = t - (a / b) * y;
	return u;
}
inline int calc(int m1, int a1,
					int m2, int a2) {
	int p1, p2;
	if (a2 - a1 < 0)
		swap(a1, a2), swap(m1, m2);
	int g = exgcd(m1, p1, m2, p2);
	if ((a2 - a1) % g)
		return -1;
	p1 *= (a2 - a1) / g, m2 /= g;
	p1 = (p1 % m2 + m2) % m2;
	return p1 * m1 + a1;
}
inline int exCRT(int k, int *a, int *m) {
	for (int i = 2; i <= k; ++i) {
		a[i] = calc(m[i - 1],
				a[i - 1], m[i], a[i]);
		if (a[i] == -1)
			return -1;
		m[i] = lcm(m[i - 1], m[i]);
	}
	return a[k] % m[k];
}
int main() {
	while (read(n)) {
		for (int i = 1; i <= n; ++i)
			read(m[i]), read(a[i]);
		print(exCRT(n, a, m), '\n');
	}
	return 0;
}
} // namespace XSC062
#undef int

__ 你,__ 去 __ 吧!

posted @ 2023-06-02 19:53  XSC062  阅读(26)  评论(0编辑  收藏  举报