中国剩余定理(CRT)及其扩展(EXCRT)详解

博客园食用效果更佳

问题背景

孙子定理是中国古代求解一次同余式方程组的方法。是数论中一个重要定理。又称中国余数定理。一元线性同余方程组问题最早可见于中国南北朝时期（公元5世纪）的数学著作《孙子算经》卷下第二十六题，叫做“物不知数”问题，原文如下：
有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二。问物几何？即，一个整数除以三余二，除以五余三，除以七余二，求这个整数。《孙子算经》中首次提到了同余方程组问题，以及以上具体问题的解法，因此在中文数学文献中也会将中国剩余定理称为孙子定理。

用现代数学的语言来分析这个问题，中国剩余定理给出了以下的一元线性同余方程组：

$\left\{\begin{array}{l}x \equiv a_{1}\left(\bmod m_{1}\right) \\ x \equiv a_{2}\left(\bmod m_{2}\right) \\ \cdot \\ \cdot\\ \cdot \\ x \equiv a_{n}\left(\bmod m_{n}\right)\end{array}\right.$

在中国剩余定理（以下简称 $CRT$ ) 中 $m_1,m_2,...,m_n$ 为两两互质的整数，在扩展中国剩余定理（以下简称 $EXCRT$ ）中 $m_1,m_2,...,m_n$ 并不满足互质条件，后者相对于前者更难解决。而两者都要运用一个算法来解决：扩展欧几里得算法。

扩展欧几里得算法

有啥用呢：

我们在中学就知道，计算两个正整数的最大公因数有两个比较常用的方法：更相减损术和辗转相除法，其中辗转相除法也叫欧几里得算法： $gcd(a,b) = gcd(b,a\%b)$ （一些题外话：个人认为以人名来命名是一种非常原始和不科学的做法，我们无法从人名中提取到关于这个定理的信息。中文在这个时候就显得博大精深了，在定义一个定理的同时十分精炼地概括了定理的精华部分，实在是妙不可言）。而扩展欧几里得算法是欧几里得算法的扩展（废话），广泛应用于 $RSA$ 加密等领域。在 $CRT$ 问题中，我们用于求得逆元和求解裴蜀等式（后面会讲的），在 $EXCRT$ 中我们用来求解最小公因数和裴蜀等式。

裴蜀等式：

裴蜀定理 $ax+by = gcd(a,b)$ （或贝祖定理,同样是看名字一脸懵的定理）得名于法国数学家艾蒂安·裴蜀，说明了对任何整数 $a、b$ 和它们的最大公约数 $gcd(a,b)$ ，关于未知数 $x$ 和 $y$ 的线性二元一次不定方程（称为裴蜀等式）：一定存在整数 $x,y$ ，使 $ax+by=gcd(a,b)$ 成立。它的一个重要推论是： $a,b$ 互质的充要条件是存在整数 $x,y$ 使 $ax+by=1$ 。证明我就略去了，来讲一下扩展欧几里得算法怎么得到裴蜀等式的一个解（有多个解，求出一个解可以写出通解）。

裴蜀等式求解过程：

应用欧几里得算法是一个递归的过程，将规模较大的问题转化为等价的小问题去解决，核心代码就是 $gcd(a,b) = gcd(b,a\%b)$ ,证明也略去了（笔者数学一般，见谅）。它的边界情况出现在 $gcd(k,0)$ 的时候，这时 $k$ 就是 $gcd(a,b)$ 啦。而扩展欧几里得算法则多了几个额外的步骤，在递归返回的时候做了一些小操作，使得我们可以得到 $ax+by=gcd(a,b)$ 的一组解。思路如下：

第一层： $ax+by=gcd(a,b)$ ，递归到下一层， $gcd(a,b) = gcd(b,a\%b)$ , 同时向下一层传 $x$ 和 $y$ ，设为 $x'$ 和 $y'$
第二层中： $bx'+a\%b*y'=gcd(b,a\%b)$ ,我们可以把 $a\%b$ 表示为 $a-\lfloor \frac{a}{b} \rfloor~b$ ,代入整理可得

$ay'+b(x'-\lfloor \frac{a}{b} \rfloor~y') = gcd(b,a\%b) = gcd(a,b)$

比较可得，第一层中的 $x$ 和 $y$ 和第二层中 $x'$ 和 $y'$ 的关系为
$x=y^{\prime} ,\quad y=x^{\prime}-\left\lfloor\frac{a}{b}\right\rfloor y^{\prime}$

第二层的 $x'$ 和 $y'$ 又可以由第三层的 $x''$ 和 $y''$ 代入同样的求得，在递归终点的时，我们只需让 $x = 1~y = 0$ 即可。注意因为上一层需要知道下一层的 $x$ 和 $y$ ，所以我们可以用传引用的方法，在递归返回时 $x$ 和 $y$ 就是下一层的 $x'$ 和 $y’$ 了。代码很短，但是不要写错细节哦（第二种写法更简洁）。

//求解 ax+by = gcd(a,b)，注意是传引用
void exGcd(ll a,ll b,ll &x,ll &y){  
    if(b == 0) { x = 1,y = 0; return;}         // b = 0时,a = gcd(原a,原b)，返回
    exGcd(b,a%b,x,y);
    ll temX = x;       //保留下一层的x'
    x = y,y = temX - a/b * y;   //x = y', y = x' - a/b * y' (x'和y'是递归下一层返回后的x和y)
}
/*void exGcd(ll a,ll b,ll &x,ll &y){  //更简洁的写法
    if(b == 0) { x = 1,y = 0; return;}
    exGcd(b,a%b,y,x);           //x和y换位
    y = y- a/b*x;
}*/

洛谷有道扩展欧几里得算法的模板题，充分理解后做出这题以后扩欧都难不住你

同余方程求解过程：

同余方程 $ax \equiv c(mod~b)$ 同样可以化成一个二元一次不定方程 $ax+by=c$ ,注意二元一次不定方程是不一定有整数解的，有整数解的充要条件是 $gcd(a,b)|c$ ，表示 $gcd(a,b)$ 是 $c$ 的一个因数。我们设 $gcd(a,b) = d$ ，则我们先用扩展欧几里得算法解出翡蜀等式 $ax+by = d$ 的解，再乘上 $c/d$ 即可得到同余方程的一个特解 (设为 $x_0,y_0$ )。要求通解也很简单，如下

$\left\{\begin{array}{l}x=x_{0}+b_{1} t \\ y=y_{0}-a_{1} t\end{array}\right.$

其中， $a_1 = \frac{a}{gcd(a,b)},b_1=\frac{b}{gcd(a,b)}$ ，~~证明当然也没有啦~~，可以手推一下的。

CRT问题的解决方法

构造出解

$CRT$ 问题主要利用了 $m_1,m_2,...,m_n$ 为两两互质的整数这一美好的性质，我们可以让

$M=m_{1} \times m_{2} \times \cdots \times m_{n}=\prod_{i=1}^{n} m_{i}$ ,同时定义 $M_{i}=M / m_{i}, \quad \forall i \in\{1,2, \cdots, n\}$

我们知道 $M_i$ 是和 $m_i$ 互质的（因为有两两互质这一性质保证）。再定义逆元： $t_{i}=M_{i}^{-1}$ 为 $M_i$ 模 $m_i$ 意义下的数论倒数（ $t_i$ 为 $M_i$ 模 $m_i$ 意义下的逆元，并不是真正的我们以前学的倒数），满足 $M_{i} t_{i} \equiv 1 \quad\left(\bmod m_{i}\right)$ ，然后我们可以构造出一个解 $x_0=\sum_{i=1}^{n} a_{i} M_{i} t_{i}$ ,因为对于 $\forall j \in[1, n]$ ：

$\left\{\begin{aligned} a_jM_jt_j\%m_i &= 0~~~&,j\ne i\\a_jM_jt_j\%m_i &= a_i~~~&,j = i \end{aligned} \right.$

所以它满足了 $x_0\%m_i=a_i$ ， $\forall i \in[1, n]$ ，解是成立的。

通解就是 $x = x_0 +kM$ , 而 $CRT$ 问题所求的最小整数解其实就是 $ans = x \% M$ 。

逆元求法

求解逆元的方法主要有如下三种，本文只介绍使用扩展欧几里得的方法。

1、费马小定理，限制模数必须为质数。 $CRT$ 问题中模数只是互质，不一定是质数，所以不可用。

2、欧拉定理，~~蒟蒻还不会~~。

3、扩展欧几里得，进入正题。

逆元构造的方程 $M_{i} t_{i} \equiv 1 \quad\left(\bmod m_{i}\right)$ 其实就是一个同余方程嘛（ $M_i$ 和 $m_i$ 互质， $gcd(M_i,m_i) = 1$ ，方程有整数解） ,代入扩展欧几里得模板就可以求出来了。

CRT问题完整代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll ;
ll mod[15],yu[15],M = 1,ans;//mod[i]即为mi,yu[i]存放模后余数

void exGcd(ll a,ll b,ll &x,ll &y){   //求解 ax+by = gcd(a,b)，注意是传引用
    if(b == 0) { x = 1,y = 0; return;}         // b = 0时,a = gcd(原a,原b)
    exGcd(b,a%b,x,y);
    ll temX = x;
    x = y,y = temX - a/b * y;   //x = y', y = x' - a/b * y' (x'和y'是递归下一层返回后的x和y)
}
/*void exGcd(ll a,ll b,ll &x,ll &y){  //更简洁的写法
    if(b == 0) { x = 1,y = 0; return;}
    exGcd(b,a%b,y,x);           //x和y换位
    y = y- a/b*x;
}*/
int main() {
    int n;
    cin>>n;       //方程组数
    for (int i = 1; i <= n ; ++i) {
        scanf("%ld %ld",&mod[i],&yu[i]);  //模数和余数，模数互质
        M*=mod[i];
    }
    for (int i = 1; i <= n ; ++i) {
        ll Mi = M / mod[i],inv,y; // Mi为所有模数乘积除以第i个模数，inv为Mi在模mi意义下的逆元
        exGcd(Mi, mod[i], inv, y);
        inv = inv % mod[i];
        ans = (ans + yu[i] * Mi * inv) % M;
    }
    cout<< (ans + M) % M;          //保证结果不出现负数
    return 0;
}

题目链接洛谷P1495 【模板】中国剩余定理(CRT)/曹冲养猪

题目链接洛谷P3868 [TJOI2009]猜数字

前者是裸题，后者要稍作变换，并且要使用快速乘。

EXCRT问题的解决方法

为啥不能继续用CRT的代码解决问题了呢？

~~显然不能~~。因为同余方程组不再满足 $m_1,m_2,...,m_n$ 为两两互质的整数这一美好的性质了，上面求出的 $M_i$ 与 $m_i$ 不一定互质，从而导致了

$\left\{\begin{aligned} a_jM_jt_j\%m_i &= 0~~~&,j\ne i\\a_jM_jt_j\%m_i &= a_i~~~&,j = i \end{aligned} \right.$

这一条件失效，逆元也求不出来。所以我们不能再用逆元来解决 $EXCRT$ 的问题，必须换一种思路。

合并同余方程组

观察一下几个简单的式子（摘自一位让我学会 $EXCRT$ 的大佬阮行止的博客，里面的数学证明表示更加严谨）

$\left\{\begin{array}{ll}x \equiv 2 & (\bmod 4) \\ x \equiv 4 & (\bmod 6)\end{array} \Rightarrow x \equiv 10 \quad(\bmod 12)\right.\\ \left\{\begin{array}{ll}x \equiv 4 &(\bmod 6) \\ x \equiv 3 &(\bmod 5)\end{array} \Rightarrow x \equiv 28 \quad(\bmod 30)\right.\\\left\{\begin{array}{ll}x \equiv 2 &(\bmod 4) \\ x \equiv 3 &(\bmod 6) \end{array} \Rightarrow \varnothing\right.$

我们可以发现，同余方程在一定条件下是可以合并的，但是也会出现无解的情况。合并后的方程仍然可以表示为同余方程的形式，并且模数是原来模数的最小公倍数（上述这些规律其实是要证明的，但是蒟蒻不太会，感兴趣的同学可以自己研究一下上面大佬的博客和其他资料）。那么，解决 $EXCRT$ 问题的关键，就在于如何合并这 $n$ 个同余方程组，并判断是否有解。

合并流程

假设前 $k-1$ 个同余方程合并得到的新方程为： $x = r_1(mod~M),$ $M$ 是前 $k-1$ 个同余方程模数的最小公倍数，现在考虑合并第 $k$ 个方程： $x = r_2(mod~m_k)$ 。

对于前 $k-1$ 个同余方程，其通解为 $x = r_1+i*M,i\in Z$ , 其中 $r_1$ 是我代码里面的 $ans$ 。我们在这个通解里找到一个 $t$ ，使得 $(r_1+t*M) \%m_k = r_2$ ,即可以满足第 $k$ 个方程。那么合并后前 $k$ 个同余方程的通解就是 $x = r_1+t*M+i*lcm(M,m_k),i \in Z$ , 再对通解模 $lcm(M,m_k)$ 即可得到新的 $ans$ 作为前 $k$ 个同余方程的最小整数解。具体实现的时候还要注意可能会出现负数和爆long long的情况，要用一些小技巧避免。

找 $t$ 的过程如下：方程整理为 $M*t+m_k*y=r_2-r_1$ ,其中 $M$ 就是翡蜀等式的 $a$ , $t$ 就是裴蜀等式的 $x$ , $r_2-r_1$ 要满足 $gcd(M,m_k)|(r_2-r_1)$ ,否则无解。我们先用老朋友扩欧算法来解出 $M*t+m_k*y=gcd(M,m_k)$ ，顺便求出 $gcd(M,m_k)$ ,再将得到的解 $t = t*\frac{r_2-r_1}{gcd(M,m_k)}$ ,即可求出 $t$ ,注意这里得用快速（龟速）乘，会爆 $long~ long$ ,而且 $r_1-r_1$ 不能是负数（否则快速乘会 $TLE$ ) , 所以这里都要使用一些取模的技巧。之后要更新 $ans$ 和 $M$ ,让 $M$ 成为前 $k$ 个同余方程模数的最下公倍数。这里有个公式是 $gcd(a,b)*lcm(a,b)=a*b$ ,可以用它来更新 $M$ 。

EXCRT 问题的完整代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
ll n,mod[100009],yu[100009];

//要用快速乘（龟速乘），防止爆long long
ll qMul(ll a,ll b,ll mo){
    ll an = 0;
    while(b) {
        if(b&1) an =(an+a) % mo;
        a = (a+a)%mo;
        b>>=1;
    }return an%mo;
}

//扩展欧几里得算法，返回gcd(a,b),并计算出ax+by = gcd(a,b)中的x和y
ll exGcd(ll a,ll b,ll &x,ll &y){
    if( b == 0 ) { x = 1;y = 0; return a;}
    ll gcd = exGcd(b,a%b,y,x);  //注意x和y的顺序
    y = y - a/b*x;
    return gcd;
}

int main() {
    ios::sync_with_stdio(false);//加速cin和cout
    cin>>n;
    for(int i = 1;i <= n;i++) cin>>mod[i]>>yu[i];
    ll ans = yu[1],M = mod[1] ,t,y;  //ans表示前i-1个方程式的特解（余数），M为前i-1个方程式的模数的最小公倍数(i从2开始)
    for(int i = 2;i <= n;i++){
        ll mi = mod[i],res = ((yu[i] - ans)%mi + mi)%mi;  //res是等式的右边部分，不能出现负数
        ll gcd = exGcd(M,mi,t,y);        //求出gcd(mi,M)
        if(res % gcd != 0) { cout<<-1<<endl;exit(0); }   //如果等式右边不能整除gcd，方程组无解
        t = qMul(t,res/gcd,mi);            //求出t还要乘上倍数，注意是快速乘取模mi (对谁取模要分清)
        ans += t * M;                               //得到前i个方程的特解（余数）
        M = mi /gcd * M;                         //M等于lcm(M,mi)，注意乘法要在除法后面做，否则会爆long long
        ans = (ans%M+M)%M;                //让特解范围限定在0~(M-1)内，防止会出现负数
    }
    cout<<ans;
    return 0;
}

裸题链接洛谷P4777 【模板】扩展中国剩余定理（EXCRT）

牛客NC15068 一个小问题

两道题的唯一区别是：前者保证一定有解，后者不一定有解。

蒟蒻初学扩展欧几里得算法和 $CRT$ 和 $EXCRT$ 问题，可能会出现一些笔误和逻辑错误，希望能得到指正。全文很多内容整理自一些大佬的博客和百度百科，旨在帮助从没有学过这些算法的同学在理解上提供尽可能多的帮助（我也是昨天才学的QAQ），如果有没有理解的地方可以私信我哦。如果有帮助,希望给我一个赞鼓励我(markdown使用过度，页面渲染好像出了点问题)~~(@^_^@)~~