同余方程组

同余方程组形式如下： \[ \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;
}