算法学习笔记(9): 中国剩余定理(CRT)以及其扩展(EXCRT)

合集 - 算法学习笔记(52)

1.算法学习笔记(∞)：杂项2023-08-17 2.算法学习笔记(1): 欧几里得算法及其扩展2023-01-12 3.算法学习笔记(2): 欧拉定理与逆元2023-01-13 4.算法学习笔记(3): 倍增2023-01-12 5.算法学习笔记(3.1): ST算法2023-10-13 6.算法学习笔记(4): 并查集及其优化2023-01-12 7.算法学习笔记(5): 最近公共祖先(LCA)2023-01-12 8.算法学习笔记(6): 树链剖分2023-01-12 9.算法学习笔记(7): 二分图2023-01-13 10.算法学习笔记(8): 网络流2023-01-13 11.算法学习笔记(8.0): 网络流前置知识2023-01-13 12.算法学习笔记(8.1): 网络最大流算法 EK, Dinic, ISAP2023-01-14 13.算法学习笔记(8.2): 上下界网络流2023-01-14 14.算法学习笔记(8.3): 网络最大流 - 模型篇2023-12-19

15.算法学习笔记(9): 中国剩余定理(CRT)以及其扩展(EXCRT)2023-01-15

扩展中国剩余定理

扩展中国剩余定理

讲解扩展之前，我们先叙述一下普通的中国剩余定理

“China Remain Theory” 也叫做孙子定理

难得是以中国命名的定理~~，谁敢不掌握~~

中国剩余定理

中国剩余定理通过一种非常精巧的构造求出了一个可行解

但是毕竟是构造，所以相对较复杂

{\begin{cases} x \equiv a_{1} (\mod m_{1}) \\ x \equiv a_{2} (\mod m_{2}) \\ x \equiv a_{3} (\mod m_{3}) \\ \dots \\ x \equiv a_{n} (\mod m_{n}) \end{cases}

$\begin{cases} x \equiv a_1 \pmod{m_1} \\ x \equiv a_2 \pmod{m_2} \\ x \equiv a_3 \pmod{m_3} \\ \dots \\ x \equiv a_n \pmod{m_n} \end{cases}$

对于上述同余方程组，首先需要满足 $m_1, m_2, m_3, \dots, m_n$ 互质

令 $m = \prod_{i=1}^{n} m_i, M_i = m / m_i,\ t_i$ 是线性同余方程 $M_it_i \equiv 1 \pmod {m_i}$ 的一个解

那么上述方程组有整数解，为 $x = \sum_{i=1}^n a_iM_it_i$

证明：

由于所有 $m_i$ 互质，所以 $M_i = m/m_i$ 是除了 $m_i$ 之外的所有模数的倍数，所以 $\forall j \ne i, a_iM_it_i \equiv 0 \pmod {m_k}$

所以代入 $x = \sum_{i=1}^n a_iM_it_i$ ，原方程组成立

中国剩余定理给出了方程的一个特解。通解可以表示为 $x + km (k \in \Z)$ 。如果需要求最小的非负整数解，只需要让 $x$ 对于 $m$ 取模就是。

参考代码

【模板】中国剩余定理（CRT）/ 曹冲养猪 - 洛谷

#include <iostream>

typedef long long ll;

ll m[11], a[11], Mm[11], y[11];

// 扩展欧几里得算法，由于求出线性同余方程 Mi * ti = 1 (mod p) 的一个特解
// 上述同余式可以转化为 Mi * ti + k * p = 1
// 所以跑一个 extgcd(Mi, p, ti, k) 即可 （k没有用的）
void extgcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) {
        x = 1, y = 0;
        return; 
    }
    extgcd(b, a % b, y, x);
    y -= (a / b) * x;
}

int main() {
    int n; scanf("%d", &n);

    ll M(1);
    // 这里读入，并求出m
    for (int i = 0; i < n; i++) {
        scanf("%lld%lld", m + i, a + i);
        M *= m[i];
    }

    // 求出Mi
    for (int i = 0; i < n; i++) {
        Mm[i] = M / m[i];
    }

    // 求出每一个不定方程的特解
    for (int i = 0; i < n; i++) {
        ll x, tmp, mi = m[i];
        extgcd(Mm[i], mi, x, tmp);
        y[i] = (x + mi) % mi;
    } 

    // 求和获得最小正整数解
    ll x(0);
    for (int i = 0; i < n; i++) {
        x = (x + a[i] * Mm[i] * y[i]) % M;
    }

    printf("%lld\n", x);
    return 0;
}

扩展中国剩余定理

还是考虑上述式子

{\begin{cases} x \equiv a_{1} (\mod m_{1}) \\ x \equiv a_{2} (\mod m_{2}) \\ x \equiv a_{3} (\mod m_{3}) \\ \dots \\ x \equiv a_{n} (\mod m_{n}) \end{cases}

$\begin{cases} x \equiv a_1 \pmod{m_1} \\ x \equiv a_2 \pmod{m_2} \\ x \equiv a_3 \pmod{m_3} \\ \dots \\ x \equiv a_n \pmod{m_n} \end{cases}$

此时，我们不保证 $m_i$ 不一定两两互质，所以中国剩余定理不再适用

所以扩展中国剩余定理就开始发挥作用了

虽然说是扩展……但是两者思路很不一样

扩展中国剩余定理采用数学归纳法，或者说两两合并法

假如我们需要合并方程

x \equiv a_{1} (\mod m_{1})

$x \equiv a_1 \pmod {m_1}$

x \equiv a_{2} (\mod m_{2})

$x \equiv a_2 \pmod {m_2}$

我们将同余式改写

k_{1} m_{1} + a_{1} = x

$k_1m_1 + a_1 = x$

k_{2} m_{2} + a_{2} = x

$k_2m_2 + a_2 = x$

将两者相减，整理后可得

k_{1} m_{1} - k_{2} m_{2} = a_{2} - a_{1}

$k_1m_1 - k_2m_2 = a_2 - a_1$

由于我们已知 $m_1, m_2, a_1, a_2$ ，也就是说我们只需要知道 $k_1, k_2$ 其中任意一个，就可以求出这时的 $x$ 是个什么东西了

所以，联想到了……扩展欧几里得算法

可以参考文章：算法学习笔记(1): 欧几里得算法及其扩展

我们令

\begin{aligned} d & = g c d (m_{1}, m_{2}) \\ c & = a_{2} - a_{1} \end{aligned}

$\begin{aligned} d &= gcd(m_1, m_2) \\ c &= a_2 - a_1 \end{aligned}$

那么化简一下有： $k_1 m_1 - k_2 m_2 = c$

所以我们求出 $k_1' \frac {m_1}d + k_2' \frac {m_2}d = 1$

那么换回去之后

\begin{aligned} k_{1} & = k_{1}^{'} \frac{c}{d} \\ k_{2} & = - k_{2}^{'} \frac{c}{d} \end{aligned}

$\begin{aligned} k_1 &= k_1' \frac cd \\ k_2 &= -k_2' \frac cd \end{aligned}$

由于 $exgcd(a, b, x, y)$ 求出的是 $xa + yb = \gcd(a, b)$ 的一组解……所以需要 $\frac cd$

令 $l = lcm(m_1, m_2) = m_1 m_2 / \gcd(m_1, m_2)$

那么此时 $x \equiv k_1' m_1\frac cd + a_1 \pmod l$

也就是可以得到 $x$ 的一个解集 $x \in X = \{kl + x|k \in \R\}$

也就是说，我们将上述两个式子合并成了

x \equiv k_{1}^{'} \frac{c}{d} m_{1} + a_{1} (\mod l c m (m_{1}, m_{2}))

$x \equiv k_1' \frac cd m_1 + a_1 \pmod {lcm(m_1, m_2)}$

那么考虑代码怎么写

typedef long long lint;

struct Equation {
    lint a, p;
};

inline lint gcd(lint x, lint y) {
    return y ? gcd(y, x % y) : x;
}

inline lint exgcd(lint x, lint y, lint &s, lint &t) {
    if (!y) {
        s = 1, t = 0;
        return x;
    }
    lint r = exgcd(y, x % y, t, s);
    t -= (x / y) * s;
    return r;
}

bool merge(const Equation &a, const Equation &b, Equation &res) {
    // 第一部分：预先求出需要的值
    lint d = gcd(a.p, b.p), s, t;
    lint lcm = a.p / d * b.p;

    // 第二部分：求出上文中的 k1', k2'
    lint c = (b.a - a.a) % lcm;
    if (c < 0) c += lcm;
    if (c % d) return false; // 判断有无解，返回false则无解
    exgcd(a.p / d, b.p / d, s, t);

    // 第三部分：求出上文中的 k1
    s = s * (c / d) % lcm;
    if (s < 0) s += lcm;

    // 第四部分：求出一个特殊的解
    lint x = (a.a + s * a.p % lcm) % lcm;
    if (x < 0) x += lcm;
    res = {x, lcm};
    return true;
}