中国剩余定理（孙子定理）：数论中的一个关于一元线性同余方程组的定理，说明了一元线性同余方程组有解的准则以及求解方法。

需要预备知识：扩展欧几里得算法

中国剩余定理

适用问题

有以下一元线性同余方程组：

$$ (S) : \quad \left\{ \begin{matrix} x \equiv a_1 \pmod {m_1} \\ x \equiv a_2 \pmod {m_2} \\ \vdots \qquad\qquad\qquad \\ x \equiv a_n \pmod {m_n} \end{matrix} \right. $$

整数 $m_1,m_2,\dots,m_n$ 其中任两数互质，求解 $x$.

求解步骤

通解可以用如下方式构造得到：

1. 设 $M = m_1 \times m_2 \times \cdots \times m_n = \prod_{i=1}^n m_i$ 是整数 $m_1,m_2,\dots,m_n$ 的乘积，并设 $M_i = M/m_i, \; \; \forall i \in \{1, 2, \cdots , n\}$，即 $M_i$ 是除了 $m_i$ 以外的 $n-1$ 个整数的乘积。
2. 设 $t_i = M_i^{-1}$ 为 $M_i$ 模 $m_i$ 的数论倒数：$t_i M_i \equiv 1 \pmod {m_i}, \; \; \forall i \in \{1, 2, \cdots , n\}.$
3. 方程组 $(S)$ 的通解形式为：$x = a_1 t_1 M_1 + a_2 t_2 M_2 + \cdots + a_n t_n M_n + k M= k M + \sum_{i=1}^n a_i t_i M_i, \quad k \in \mathbb{Z}.$ 在模 $M$ 的意义下，方程组 $(S)$ 只有一个解：$x = \sum_{i=1}^n a_i t_i M_i.$

来源：Wikipedia (使用 CC BY-SA 3.0 协议)

其中，求数论倒数（即乘法逆元）的方式是扩展欧几里得算法，不可使用费马小定理。

简单证明

$$ \begin{align} x &= a_1 t_1 M_1 + a_2 t_2 M_2 + \cdots + a_n t_n M_n + k M\\ &= k M + \sum_{i=1}^n a_i t_i M_i, \quad k \in \mathbb{Z}. \end{align} $$

我们验证该解是否满足原一元线性同余方程组，首先考虑方程组中第一个方程：

$$ x\bmod m_1=(a_1 t_1 M_1) \bmod m_1+ (a_2 t_2 M_2) \bmod m_1+ \cdots + (a_n t_n M_n) \bmod m_1+ (k M)\bmod m_1 $$

首先，由于 $t_1 M_1 \equiv 1 \pmod {m_1}$，因此 $(a_1 t_1 M_1) \bmod m_1=a_1\bmod m_1$. 另外由于 $M_i,\; i\in\{2,3,\cdots,n\}$ 中均包含 $m_1$ 因子，因此 $(a_i t_i M_i) \bmod m_1=0,\; i\in\{2,3,\cdots,n\}$. 综上 $x\bmod m_1=a_1\bmod m_1$，即 $x \equiv a_1 \pmod {m_1}$.

同理可得，该解满足方程组中的所有方程，因此它是方程组的一个解。

扩展中国剩余定理

我们注意到，中国剩余定理的使用条件是整数 $m_1,m_2,\dots,m_n$ 其中任两数互质。如果不满足该条件，中国剩余定理失效。

针对这种情况，我们使用扩展中国剩余定理求解。（该情况有可能无解）

适用问题

有以下一元线性同余方程组：

$$ (S) : \quad \left\{ \begin{matrix} x \equiv a_1 \pmod {m_1} \\ x \equiv a_2 \pmod {m_2} \\ \vdots \qquad\qquad\qquad \\ x \equiv a_n \pmod {m_n} \end{matrix} \right. $$

不保证整数 $m_1,m_2,\dots,m_n$ 其中任两数互质，问：方程组是否有解？若有解，试求 $x$.

推导和求解步骤

联立其中两个方程

首先考虑 $(S)$ 中前两个方程：

$$ \left\{ \begin{matrix} x \equiv a_1 \pmod {m_1} \\ x \equiv a_2 \pmod {m_2} \end{matrix} \right.\tag{1} $$

$(1)$ 经过变形可得到：

$$ \left\{ \begin{matrix} x=a_1+m_1k_1\\ x=a_2+m_2k_2\\ \end{matrix} \right.\tag{2} $$

求解联立方程组

消去 $(2)$ 中 $x$ 并进行变形，得：

$$ m_1k_1-m_2k_2=a_2-a_1\tag{3} $$

该式对应的裴蜀等式为 $(4)$，其中 $d=\gcd(m_1,m_2)$：

$$ m_1k_1'-m_2k_2'=d\tag{4} $$

通过扩展欧几里得算法，求得 $k_1',k_2',d$，若 $a_2-a_1$ 是 $d$ 的倍数，即原式有解，一个特解为：

$$ \left\{ \begin{matrix} k_1^*=k_1'\cdot(a_2-a_1)/d\\ k_2^*=k_2'\cdot(a_2-a_1)/d \end{matrix} \right. $$

若 $a_2-a_1$ 不是 $d$ 的倍数，那么原一元线性同余方程组无解，可停止下面的计算。得到特解后，我们可以构造出通解（构造方式见扩展欧几里得算法的笔记）：

$$ \left\{ \begin{matrix} k_1=k_1^*+t\cdot (m_2/d)\\ k_2=k_2^*-t\cdot (m_1/d)\\ \end{matrix} \right.\; ,t\in\mathbb{Z} $$

将解回代

将解出的 $k_1$ 回代到 $(2)$ 中第一个方程，得到：

$$ \begin{align} x&=a_1+m_1k_1\\ &=a_1+m_1\cdot[k_1^*+t\cdot (m_2/d)]\\ &=(a_1+m_1k_1^*)+(m_1m_2/d)\cdot t\\ &=(a_1+m_1k_1^*)+\mathrm{lcm}(m_1,m_2)\cdot t\\ \end{align}\tag{5} $$

我们可以发现，$(5)$ 式与 $(2)$ 式的结构非常相似，结构均为：

$$ x=C_1+C_2\cdot m\;(C_1,C_2为常数,m为自变量)\tag{6} $$

将方程组合并为一个方程

我们上面的过程，即可看作将 $(2)$ 的两个式子合并为了 $(6)$ 这样的一个式子，$(6)$ 和 $(2)$ 同解。那么我们可以不断重复这个步骤，将式子不断合并。即再将 $(6)$ 与 $(S)$ 中第三个式子合并得到新式子，再将新式子与 $(S)$ 中第四个式子合并……

合并最后的结果便是一个式子，我们假设这个式子为：

$$ x=\alpha+b\cdot\beta\tag{7} $$

$(7)$ 式非常简单就能解出结果，而且原方程组 $(S)$ 与 $(7)$ 同解，那么我们解 $(7)$ 这个方程即可得到 $(S)$ 的解。

求解最终合并式

若我们要求得最小的正整数解，$(7)$ 式的解其实就是 $\alpha\bmod\beta$

模板

typedef long long ll;

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

// 如果中途输出-1,-1则需要立即停止
// first-a second-mod
pair<ll, ll> excrt(pair<ll, ll> a, pair<ll, ll> b)
{
    ll ya, yb;
    ll d = exgcd(a.second, b.second, ya, yb);
    if ((b.first - a.first) % d)
        return {-1, -1};
    ya *= (b.first - a.first) / d;
    ll tmp = b.second / d;
    ya = (ya % tmp + tmp) % tmp;
    a.first += a.second * ya;
    a.second = a.second / d * b.second;
    return a;
}

典型例题

原题链接：https://www.acwing.com/problem/content/206/

注意：

• 为防止乘法溢出，建议先除后乘。
• C 语言负数取模会成负数，因此取模需要写成 ($\alpha\bmod\beta+\beta)\bmod\beta$ 这种形式。
• 程序实现时并没有使用上述过程中的那么多变量，同一个变量被重复利用。

#include <bits/stdc++.h>

using namespace std;

typedef long long LL;
const int MAXN = 30;
int n;

LL exgcd(LL a, LL b, LL &x, LL &y)
{
    if (!b)
    {
        x = 1;
        y = 0;
        return a;
    }
    LL d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

int main()
{
    cin >> n;
    LL a1, m1;
    cin >> a1 >> m1;
    for (int i = 2; i <= n; i++)
    {
        LL a2, m2, y1, y2;
        cin >> a2 >> m2;
        LL d = exgcd(a1, a2, y1, y2);
        if ((m2 - m1) % d)
        {
            cout << -1 << endl;
            return 0;
        }
        y1 *= (m2 - m1) / d;
        LL tmp = a2 / d;
        y1 = (y1 % tmp + tmp) % tmp;
        m1 += a1 * y1;
        a1 = abs(a1 / d * a2);
    }
    cout << (m1 % a1 + a1) % a1 << endl;
    return 0;
}