6.2.3 中国剩余定理（CRT & exCRT）

中国剩余定理

中国剩余定理（CRT）

孙子算经

今有物不知其数，三三数之余二；五五数之余三；七七数之余二。问物几何？

像这样的问题，我们可以将它简化为下面的一个同余方程组

\[\left\{ \begin{align*} x \equiv 2 \ (\bmod3) \\ x \equiv 3 \ (\bmod5) \\ x \equiv 2 \ (\bmod7) \\ \end{align*} \right. \]

这个方程组满足条件，所有模数两两互质

我们并不好直接求解这个同余方程组，所以把它转化为下面的形式

\[\left\{ \begin{align*} x_1 \equiv 1 \ (\bmod3) \\ x_1 \equiv 0 \ (\bmod5) \\ x_1 \equiv 0 \ (\bmod7) \\ \end{align*} \right. \]

\[\left\{ \begin{align*} x_2 \equiv 0 \ (\bmod3) \\ x_2 \equiv 1 \ (\bmod5) \\ x_2 \equiv 0 \ (\bmod7) \\ \end{align*} \right. \]

\[\left\{ \begin{align*} x_3 \equiv 0 \ (\bmod3) \\ x_3 \equiv 0 \ (\bmod5) \\ x_3 \equiv 1 \ (\bmod7) \\ \end{align*} \right. \]

那么对于这三个方程组的解，如果已经解出来了，那么最后的解就是 \(x_1r_1 + x_2r_2 + x_3r_3 (\bmod m_1m_2m_3)\)

那么对于上面的每个方程组，我们又该如何求解呢？

我们可以先看最后一个，其他两个是一样的，相当于解这样一个同余方程

\[15x_3 \equiv 1 \ (\bmod 7) \\ 15x_3 = 1 + 7y_3 \\ 15x_3 - 7y_3 = 1 \]

到这里，我们就可以用 \(exgcd\) 来求解 \(x_3\) 具体是多少

写成一般形式

\[\left\{ \begin{align*} x \equiv r_1 \ (\bmod m_1) \\ x \equiv r_2 \ (\bmod m_2) \\ x \equiv r_3 \ (\bmod m_3) \\ \end{align*} \right. \]

设 \(M = \prod_{i = 1}^nm_i\) ，对每个式子我们都可以写出一个这样的同余方程

\[\frac M {m_i} x_i - m_i y_i = 1 \]

因为模数两两互质，所以 \(\gcd (\frac M{m_i},m_i) = 1\) ，所以 \(\gcd (\frac M{m_i},m_i) \mid 1\) ，所以保证一定有解

那么我们最终的答案就是 \(\displaystyle \sum_{i = 1}^n x_ir_i\)

code

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
using namespace std;
const int N = 1e5 + 10;
namespace michaele {
    typedef __int128 i8;
    typedef long long ll;
    const int N = 1e5 + 10;
    
    int n, r[N], m[N];
    void exgcd (ll a, ll b, ll &x, ll &y) {
        if (b == 0) {
            x = 1;
            y = 0;
            return;
        }
        exgcd (b, a % b, x, y);
        ll t = x;
        x = y;
        y = t - a / b * y;
    }
    void main () {
        cin >> n;
        ll M = 1;
        for (int i = 1; i <= n; i++) {
            cin >> m[i] >> r[i];
            M *= m[i];
        }
        ll ans = 0, x, y;
        for (int i = 1; i <= n; i++) {
            ll c = M / m[i];
            exgcd (c, m[i], x, y);
            ans += (i8)x * c % M * r[i] % M;
            ans %= M;
        }
        cout << (ans % M + M) % M << endl;
    }
}

int main () {
    michaele :: main ();
    return 0;
}

中国剩余定理使用的前提是所有 \(m_i\) 两两互质，这个条件其实是非常强的，如果不满足这个条件，又该如何求解呢？

扩展中国剩余定理（EXCRT）

EXCRT 就是用来解决不保证 \(m_i\) 两两互质的情况的

它的主要思想是将方程组中方程选两个进行合并，合并成一个，不断进行合并，合并 \(n - 1\) 次，最后得到一个方程，就是我们的解系

考虑只有两个方程的情况

\[\left\{ \begin{align*} x \equiv r_1 \ (\bmod m_1) \\ x \equiv r_2 \ (\bmod m_2) \end{align*} \right. \]

对于第一个方程而言，解的形式为 \(x = r_1 + k_1m_1\)

带入第二个方程，得到 \(r_1 + k_1m_1 = r_2 + k_2m_2\) ，这里只有 \(k_1,k_2\) 未知

那这个东西其实就很像是 \(ax + by = c\) 的形式，给它变个形 \(k_1m_1 - k_2m_2 = r_2 - r_1\)

如果要用 \(exgcd\) 求解，我们还需要先判断一下这个方程有没有解，设 \(d = \gcd (m_1,m_2)\)

• 如果 \(d\) 不整除 \((r_2 - r_1)\) 那么直接退出返回无解即可
• 否则有解

设 \(p_1 = m1 / d,p_2 = m_2 / d\) ，刚才的方程转化为 \(k_1p_1 - k_2p_2 = \frac {r_2 - r_1} d\)

这时 \(\gcd (p_1, p_1) = 1\) 我们用 \(exgcd\) 求解 \(t_1p_1 - t_2p_2 = 1\) 的方程的解 \(t_1\)，得到 \(k_1 = t_1 \frac{r_2 - r_1} d\) ，设 \(K_1\) 为 \(k_1\) 的通解，那么 \(K_1 = k_1 + z\frac {m_2} d\)

那么 \(x = K_1m_1 + r_1 = m_1k_1 + z \frac {m_1m_2} d + r_1\) ，等价于 \(x \equiv m_1k_1 + r_1 \ (\bmod \frac {m_1m_2} d)\)

这样，我们就成功将两个方程合并成为了一个方程，最后的解其实就是合并成一个方程以后的 \(r\)

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>

using namespace std;

namespace michaele {

    typedef long long ll;
    typedef __int128 i8;

    const int N = 1e5 + 10;

    ll n, r[N], m[N];

    ll gcd (ll a, ll b) {
        int az = __builtin_ctzll (a);
        int bz = __builtin_ctzll (b);
        int z = min (az, bz);
        ll diff;
        b >>= bz;
        while (a) {
            a >>= az;
            diff = a - b;
            if (!diff) break;
            az = __builtin_ctzll (diff);
            b = a < b ? a : b;
            a = diff < 0 ? -diff : diff;
        }
        return b << z;
    }

    void exgcd (ll a, ll b, ll &x, ll &y) {
        if (b == 0) {
            x = 1;
            y = 0;
            return;
        }
        exgcd (b, a % b, x, y);
        ll t = x;
        x = y;
        y = t - a / b * y;
    }

    ll excrt () {
        cin >> n;
        for (int i = 1; i <= n; i++) {
            scanf ("%lld%lld", &m[i], &r[i]);
        }
        ll m1, r1, m2, r2, d, k1, k2, p1, p2;
        m1 = m[1], r1 = r[1];
        for (int i = 2; i <= n; i++) {
            m2 = m[i], r2 = r[i];
            d = gcd (m1, m2);
            if ((r2 - r1) % d) return -1;
            p1 = m1 / d;
            p2 = m2 / d;
            exgcd (p1, p2, k1, k2);
            k1 = (i8) k1 * (r2 - r1) / d % (m2 / d);
            k1 = (k1 % (m2 / d) + (m2 / d)) % (m2 / d);
            r1 = r1 + k1 * m1;
            m1 = m1 / d * m2;
        }
        return (r1 % m1 + m1) % m1;
    }

    void main () {
        printf ("%lld\n", excrt ());
    }
}

int main () {
    michaele :: main ();
    return 0;
}