发布于 ,更新于 

同余方程组

同余方程组形式如下: \[ \begin{cases} x \equiv a_1 \pmod {n_1} \newline x \equiv a_2 \pmod {n_2} \newline \dots \newline x \equiv a_k \pmod {n_k} \end{cases} \]

中国剩余定理

原理

本算法用于所有 \(n_i\) 两两互质的情况。

实际上类似构造,先把 \(n_1\) 累积起来 \[ M = \prod_{i=1}^k n_i \] \[ M_i = \frac M{n_i} \]

先考虑如何构造出一个同余方程的特解,设 \(M_i ^{-1}\)\(M_i\)\(\bmod n_i\) 意义下的逆元,想到 \(M_iM_i^{-1} = 1\pmod {n_i}\),于是有 \(x = a_i M_i M_i^{-1}\) 满足方程 \(x \equiv a_i \pmod {n_i}\)
这实际上就是上面定义 \(M_i\) 为质数的原因:需要保证有逆元存在。

至于解的合并,把 \(x\) 相加即可。即 \(\sum_{i=1}^k a_iM_iM_i^{-1}\)。通解即为 \(nM+ \sum_{i=1}^k a_iM_iM_i^{-1}, n\in \mathbb Z\) 。为什么能这样合并?因为对于 \(i,j,i\neq j\) 来说,必然存在 \(n_j \mid M_i\),因为 \(n_i\) 两两互质,\(M_i\) 的因子包含 \(n_j\),因此不会对已经推出来的解造成任何影响。

通解方面:容易想到,加减 \(nM\) 也不会对解的可行性造成任何影响。

实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
llint ex_gcd(llint u, llint v, llint& x, llint& y) {
if (!v) {
x = 1, y = 0;
return u;
}

llint g = ex_gcd(v, u % v, x, y);
llint temp = x;
x = y;
y = temp - u / v * y;
return g;
}

llint crt() {
llint tmp = 1, ans = 0;
for (int i = 1; i <= n; ++ i) tmp *= b[i];
for (int i = 1; i <= n; ++ i) {
llint m = tmp / b[i], x, y;
ex_gcd(m, a[i], x, y);
ans = (ans + n[i] * m * x % tmp) % tmp;
}
return (ans % tmp + tmp) % tmp;
}

扩展中国剩余定理/同余方程合并

原理

\(n_i\) 不互质的时候,上面做法就不再适用(不能保证逆元存在),此时我们考虑逐一合并同余方程组。
方法是使用 \(\operatorname{ex\_gcd}\)

假设方程组只有两个方程

\[ \begin{cases} x \equiv a_1 \pmod{n_1} \newline x \equiv a_2 \pmod{n_2} \newline \end{cases} \]

可以写作

\[ \begin{cases} x = a_1 + k_1n_1 \newline x = a_2 + k_2n_2 \end{cases} \]\(a_1 + k_1n_1 = a_2 + k_2n_2\)
移项,得 \(k_1n_1 - k_2n_2 = a_2-a_1\) ,接下来按照扩展欧几里得的思路推导。
\(g = \gcd(n_1,n_2)\),若使用扩展欧几里得算法,可以求出 \(n_1k_1 + n_2(-k_2) = \gcd(n_1,n_2)\) 的一组特解 \(k_1, k_2\),即 \[ \begin{aligned} n_1k_1 + n_2(-k_2) &= g \newline (k_1 \cdot \frac{a_2-a_1}{g})n_1 + (-k_2 \cdot \frac{a_2-a_1}{g})n_2 &= a_2-a_1 \end{aligned} \] 抽出刚才表示 \(x\) 的方程组的一条,得到 \(x_0 = a_1 + k_1n_1\)\(k_1\) 的通解为 \(k_1 + \cfrac{n_2}gp, p\in \mathbb Z\),则 \[ \begin{aligned} x &= a_1 + (k_1+\cfrac{n_2}gp)n_1 \newline &= n_1k_1 + a_1 + \cfrac{n_1n_2}gp \newline &= n_1k_1 + a_1 + \operatorname{lcm}(n_1,n_2)\cdot p \end{aligned} \] 显然我们可以把这个等式写成同余式的形式(没必要)

\[ \begin{aligned} x_0 &\equiv x \pmod {\operatorname{lcm}(n_1,n_2)} \newline x&\equiv n_1k_1 + a_1 \pmod {\operatorname{lcm}(n_1,n_2)} \end{aligned} \] \(\operatorname{lcm}(n_1,n_2),n_1k_1+a_1\) 都是已知的,因此方程可解,合并完成。

实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#define N 100005

i128 ex_gcd(i128 a, i128 b, i128& x, i128& y) {
if (b) {
i128 tmp;
i128 g = ex_gcd(b, a % b, x, y);
tmp = x;
x = y;
y = tmp - a / b * y;
return g;
} else {
x = 1, y = 0;
return a;
}
}

i64 n;
i128 a[N], m[N]; // a === x (mod m)

i128 ex_crt() {
i128 k1, k2, a1, m1, a2, m2, c, gc, g;
a1 = a[1], m1 = m[1];
for (i64 i = 2; i <= n; ++i) {
a2 = a[i], m2 = m[i], c = ((a2 - a1) % m2 + m2) % m2;
g = ex_gcd(m1, m2, k1, k2);
gc = m2 / g;
k1 = k1 % gc * ((c / g) % gc) % gc;
a1 = k1 * m1 + a1;
m1 *= gc;
a1 = (a1 % m1 + m1) % m1;
}
return a1;
}

int main() {
read(n);
for (i64 i = 1; i <= n; ++i) {
read(m[i], a[i]);
}

write(ex_crt(), '\n');
return 0;
}