同余方程组

同余方程组形式如下： $⎧ {⎨ {⎩ 𝑥 \equiv 𝑎 1 (m o d 𝑛 1) 𝑥 \equiv 𝑎 2 (m o d 𝑛 2) \dots 𝑥 \equiv 𝑎 𝑘 (m o d 𝑛 𝑘)$

中国剩余定理

原理

本算法用于所有 $𝑛 𝑖$ 两两互质的情况。

实际上类似构造，先把 $𝑛 1$ 累积起来 $𝑀 = 𝑘 \prod 𝑖 = 1 𝑛 𝑖$ $𝑀 𝑖 = 𝑀 𝑛 𝑖$

先考虑如何构造出一个同余方程的特解，设 $𝑀 - 1 𝑖$ 为 $𝑀 𝑖$ 在 $m o d 𝑛 𝑖$ 意义下的逆元,想到 $𝑀 𝑖 𝑀 - 1 𝑖 = 1 (m o d 𝑛 𝑖)$ ，于是有 $𝑥 = 𝑎 𝑖 𝑀 𝑖 𝑀 - 1 𝑖$ 满足方程 $𝑥 \equiv 𝑎 𝑖 (m o d 𝑛 𝑖)$ 。
这实际上就是上面定义 $𝑀 𝑖$ 为质数的原因：需要保证有逆元存在。

至于解的合并，把 $𝑥$ 相加即可。即 $\sum 𝑘 𝑖 = 1 𝑎 𝑖 𝑀 𝑖 𝑀 - 1 𝑖$ 。通解即为 $𝑛 𝑀 + \sum 𝑘 𝑖 = 1 𝑎 𝑖 𝑀 𝑖 𝑀 - 1 𝑖, 𝑛 \in ℤ$ 。为什么能这样合并？因为对于 $𝑖, 𝑗, 𝑖 \neq 𝑗$ 来说，必然存在 $𝑛 𝑗 ∣ 𝑀 𝑖$ ，因为 $𝑛 𝑖$ 两两互质， $𝑀 𝑖$ 的因子包含 $𝑛 𝑗$ ，因此不会对已经推出来的解造成任何影响。

通解方面：容易想到，加减 $𝑛 𝑀$ 也不会对解的可行性造成任何影响。

实现

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

复制

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

原理

当 $𝑛 𝑖$ 不互质的时候，上面做法就不再适用（不能保证逆元存在），此时我们考虑逐一合并同余方程组。
方法是使用 $e x_g c d$

假设方程组只有两个方程

${𝑥 \equiv 𝑎 1 (m o d 𝑛 1) 𝑥 \equiv 𝑎 2 (m o d 𝑛 2)$

可以写作

${𝑥 = 𝑎 1 + 𝑘 1 𝑛 1 𝑥 = 𝑎 2 + 𝑘 2 𝑛 2$ 则 $𝑎 1 + 𝑘 1 𝑛 1 = 𝑎 2 + 𝑘 2 𝑛 2$ 。
移项，得 $𝑘 1 𝑛 1 - 𝑘 2 𝑛 2 = 𝑎 2 - 𝑎 1$ ，接下来按照扩展欧几里得的思路推导。
设 $𝑔 = g c d (𝑛 1, 𝑛 2)$ ，若使用扩展欧几里得算法，可以求出 $𝑛 1 𝑘 1 + 𝑛 2 (- 𝑘 2) = g c d (𝑛 1, 𝑛 2)$ 的一组特解 $𝑘 1, 𝑘 2$ ，即 $𝑛 1 𝑘 1 + 𝑛 2 (- 𝑘 2) = 𝑔 (𝑘 1 \cdot 𝑎 2 - 𝑎 1 𝑔) 𝑛 1 + (- 𝑘 2 \cdot 𝑎 2 - 𝑎 1 𝑔) 𝑛 2 = 𝑎 2 - 𝑎 1$ 抽出刚才表示 $𝑥$ 的方程组的一条，得到 $𝑥 0 = 𝑎 1 + 𝑘 1 𝑛 1$ ， $𝑘 1$ 的通解为 $𝑘 1 + 𝑛 2 𝑔 𝑝, 𝑝 \in ℤ$ ，则 $𝑥 = 𝑎 1 + (𝑘 1 + 𝑛 2 𝑔 𝑝) 𝑛 1 = 𝑛 1 𝑘 1 + 𝑎 1 + 𝑛 1 𝑛 2 𝑔 𝑝 = 𝑛 1 𝑘 1 + 𝑎 1 + l c m (𝑛 1, 𝑛 2) \cdot 𝑝$ 显然我们可以把这个等式写成同余式的形式（没必要）

$𝑥 0 \equiv 𝑥 (m o d l c m (𝑛 1, 𝑛 2)) 𝑥 \equiv 𝑛 1 𝑘 1 + 𝑎 1 (m o d l c m (𝑛 1, 𝑛 2))$ $l c m (𝑛 1, 𝑛 2), 𝑛 1 𝑘 1 + 𝑎 1$ 都是已知的，因此方程可解，合并完成。

实现

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

复制

本文采用 GNU 自由文档许可证许可协议，转载请注明出处。