中国剩余定理

文章转载自算法学习笔记(10): 中国剩余定理 - 知乎 (zhihu.com)，如有侵权，请联系作者删除。

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

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

实际上就是求解下面的同余方程组： \[ \left\{ \begin{aligned} x & \equiv b_1 & (mod\ a_1)\\ x & \equiv b_2 & (mod\ a_2)\\ x & \equiv b_n & (mod\ a_n) \end{aligned} \right. \] 这个方程组有解的一个充分条件是：$a_1、a_2······a_n$两两互质，可以使用构造法得到下面方程的通解，一下面这道模板题为例：

（洛谷P1495 曹冲养猪）

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

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

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

要想直接找到一个$n$使得方程组$\left\{ \begin{aligned} x & \equiv 2 &(mod\ 3)\\ y & \equiv 3 &(mod\ 5)\\ z & \equiv 2 &(mod\ 7) \end{aligned} \right.$ 成立时不容易的，但是要找到

$n_1,n_2,n_3$使得$\left\{ \begin{aligned} n_1 & \equiv 2& (mod\ 3) \\ n_2 & \equiv 3& (mod\ 5) \\ n_3 & \equiv 2& (mod\ 7) \end{aligned} \right.$是相对容易的。

直接令$x=n_1+n_2+n_3$可以吗？恐怕未必。在什么情况下$n_1\equiv 2(mod\ 3)$可以推出$n_1+n_2\equiv2(mod\ 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的倍数。也就是说，现在我们需要求三个同余方程：

$\left\{ \begin{aligned} 35m_1 & \equiv 2 & (mod\ 3) \\ 21m_2 & \equiv 3 & (mod\ 5) \\ 15m_3 & \equiv 2 & (mod\ 7) \end{aligned} \right.$

注意到模数两两互质，则$gcd(35,3)=gcd(21,5)=gcd(15,7)=1$，所以我们可以用拓展欧几里得的方法解（我愿称之为最妙的一步qwq，这里用的是求逆元的方法）：

$ { \[\begin{aligned} 35w_1 & \equiv 1& (mod\ 3) \\ 21w_2 & \equiv 1& (mod\ 5)\\ 15w_3 & \equiv 1& (mod\ 7) \end{aligned}\]

. $

解得$w_1=2,w_2=1,w_3=1$，然后可得：$\left\{ \begin{aligned} m_1 & =2w_1=4 \\ m_2 & =3w_2=3 \\ m_3 & =2w_3=2 \end{aligned} \right.$

于是：$\left\{ \begin{aligned} n_1 & = 35m_1=140 \\ n_2 & = 21m_2=63 \\ n_3 & = 15m_3=30 \end{aligned} \right.$。

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

接下来将这个过程一般化（个人感觉严重劝退qwq）。设$p=\prod_{i=1}^{n}a_i$（即所有模数的乘积），并设$r_i=\frac{p}{a_i}$（在“物不知数”中即为35、21和15）。于是$w_i={inv(r_i)|}_{a_i}$（表示$r_i$在模$a_i$意义下的逆元），$m_i=b_iw_i$，而$n_i=r_im_i$，所有$n_i$相加得到$x$。

以上这些综合起来就是：

$x\equiv \sum_{i=1}^nb_ir_i{[r_i]^{-1}|}_{a_i}\ (mod\ 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;
}

自己的实现（其实就是照抄作者的qwq）：

#include<bits/stdc++.h>
using namespace std;
using ll=long long;
const int N=1e5+10;
ll a[10],b[10];
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;
}
ll inv(ll a,ll p){
    ll x,y;
    exgcd(a,p,x,y);
    return (x%p+p)%p;
}
ll CRT(ll a[],ll b[],ll n){
    ll p=1,x=0;
    for(int i=1;i<=n;i++){
        p*=a[i];
    }
    for(int i=1;i<=n;i++){
        ll r=p/a[i];
        x+=(b[i]*r*inv(r,a[i]))%p;
    }
    return x%p;
}
int main(){
    int n;
    cin>>n;
    for(int i=1;i<=n;i++){
        cin>>a[i]>>b[i];
    }
    cout<<CRT(a,b,n);
    return 0;
}