算法学习笔记(10): 中国剩余定理

中国剩余定理，也叫孙子定理，之所以叫这个名字，是因为《孙子算经》中有这样一个问题：

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

这个被叫做“物不知数”的问题本质上是解下面的同余方程组：

$\begin{cases}x\equiv b_1 \pmod{a_1}\\x\equiv b_2 \pmod{a_2}\\...\\x\equiv b_n \pmod{a_n}\end{cases}$

后来的数学家在研究中发现，这一方程组有解的一个充分条件是 $a_1, \ a_2, \ ... \ a_n$ 两两互质，并用构造法给出了这种情况下方程的通解。而这种方法在算法竞赛中也常常会用到，例如这道模板题：

（洛谷P1495 曹冲养猪）

题目描述
自从曹冲搞定了大象以后，曹操就开始捉摸让儿子干些事业，于是派他到中原养猪场养猪，可是曹冲满不高兴，于是在工作中马马虎虎，有一次曹操想知道母猪的数量，于是曹冲想狠狠耍曹操一把。举个例子，假如有16头母猪，如果建了3个猪圈，剩下1头猪就没有地方安家了。如果建造了5个猪圈，但是仍然有1头猪没有地方去，然后如果建造了7个猪圈，还有2头没有地方去。你作为曹总的私人秘书理所当然要将准确的猪数报给曹总，你该怎么办？
输入格式
第一行包含一个整数n (n <= 10) – 建立猪圈的次数，解下来n行，每行两个整数ai, bi( bi <= ai <= 1000), 表示建立了ai个猪圈，有bi头猪没有去处。你可以假定ai,aj互质.
输出格式
输出包含一个正整数，即为曹冲至少养母猪的数目。

（本质上就是给“物不知数”套了个背景。）

我们从“物不知数”这个具体问题出发。（接下来一大波数学公式，请做好心理再看）

要想直接找到一个 $x$ 使得方程组 $\begin{cases}x\equiv 2\pmod{3}\\x\equiv 3\pmod{5}\\x\equiv 2\pmod{7}\end{cases}$ 成立当然是不容易的，但是要找到 $n_1, \ n_2,\ n_3$ 使得 $\begin{cases}n_1\equiv 2\pmod{3}\\n_2\equiv 3\pmod{5}\\n_3\equiv 2\pmod{7}\end{cases}$ 是相对容易的。

那么令 $x=n_1+n_2+n_3$ 可以吗？那恐怕未必。在什么情况下 $n_1\equiv 2\pmod{3}$ 可以推出 $n_1+n_2\equiv 2\pmod{3}$ 呢？显然，那只有当 $n_2$ 是3的倍数时成立。同理，要使 $n_1+n_2+n_3$ 也符合前式，需要 $n_2$ 和 $n_3$ 都是3的倍数。

这样推下去， $x=n_1+n_2+n_3$ 符合方程组的条件是 $n_1$ 是35的倍数， $n_2$ 是21的倍数， $n_3$ 是15的倍数。也就是说，现在我们只需要解三个同余方程：

$\begin{cases}35m_1\equiv 2\pmod{3}\\21m_2\equiv 3\pmod{5}\\15m_3\equiv 2\pmod{7}\end{cases}$ 。

注意到模数两两互质，则 $\gcd(35,3)=\gcd(21,5)=\gcd(15,7)=1$ ，所以我们可以用拓展欧几里得的方法解

$\begin{cases}35w_1\equiv 1\pmod{3}\\21w_2\equiv 1\pmod{5}\\15w_3\equiv 1\pmod{7}\end{cases}$ 。（其实相当于求逆元）

解得 $w_1=2,\ w_2=1,\ w_3=1$ ，然后可得 $\begin{cases}m_1=2w_1=4\\ m_2=3w_2=3\\ m_3=2w_3=2\end{cases}$ ，于是 $\begin{cases}n_1=35m_1=140\\ n_2=21m_2=63\\ n_3=15m_3=30\end{cases}$ 。

三者相加，即得一特解233（这里的233不是网络意义下的233，但我算出来不禁233了）。所有与233在模105意义下同余的数都是这个方程组的解，要求最小正数解只需对105取模即可，这里得出来是23。

现在我们把刚刚这个过程一般化。我们设 $p=\prod_{i=1}^{n}a_i$ （即所有模数的乘积），并设 $r_i=\frac{p}{a_i}$ （在“物不知数”中即为35、21和15）。于是 $w_i=\left.\text{inv}\left(r_i\right)\right|_{a_i}$ （表示 $r_i$ 在模 $a_i$ 意义下的逆元）， $m_i=b_iw_i$ ，而 $n_i=r_im_i$ ，所有 $n_i$ 相加即得 $x$ 。

我们把以上这些综合成一个（看起来可能有点劝退的）公式即是：

$x \equiv \sum_{i=1}^{n}b_i{r_i\left.\left[r_i\right]^{-1}\right|_{a_i}}\pmod{p}$

现在我们来看（可能相对没那么劝退的）代码吧：

inline ll CRT(ll a[], ll b[], ll n)  // a是模数数组，b是余数数组，n是数组长度
{
    ll p = 1, x = 0;
    for (int i = 0; i < n; ++i)
        p *= a[i];
    for (int i = 0; i < n; ++i)
    {
        ll r = p / a[i];
        x += (b[i] * r * inv(r, a[i])) % p;   // 逆元的求法参见上篇文章，或者下面有完整代码
    }
    return x % p;
}

这个函数返回的是符合方程组的最小正数解，一般要求的正是这个。

再附上曹冲养猪的完整AC代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
ll exgcd(ll a, ll b, ll &x, ll &y)
{
    if (b == 0)
    {
        x = 1;
        y = 0;
        return a;
    }
    ll d = exgcd(b, a % b, y, x);
    y -= (a / b) * x;
    return d;
}
inline ll inv(ll a, ll p)
{
    ll x, y;
    exgcd(a, p, x, y);
    return (x % p + p) % p;
}
inline ll CRT(ll a[], ll b[], ll n)
{
    ll p = 1, x = 0;
    for (int i = 0; i < n; ++i)
        p *= a[i];
    for (int i = 0; i < n; ++i)
    {
        ll r = p / a[i];
        x += (b[i] * r * inv(r, a[i])) % p;
    }
    return x % p;
}
int main()
{
    ll n, a[10], b[10];
    scanf("%lld", &n);
    for (int i = 0; i < n; ++i)
        scanf("%lld%lld", a + i, b + i);
    printf("%lld\n", CRT(a, b, n));
    return 0;
}