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

问题背景

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

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

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

扩展欧几里得算法

有啥用呢：

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

裴蜀等式：

  裴蜀定理\(ax+by = gcd(a,b)\)（同样是看名字一脸懵的定理，还有这个字念"pei",不知打错了多少次“翡”）得名于法国数学家艾蒂安·裴蜀，说明了对任何整数\(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)\) 的一组解。思路如下：

1. 第一层： \(ax+by=gcd(a,b)\)，递归到下一层，\(gcd(a,b) = gcd(b,a\%b)\), 同时向下一层传 \(x\) 和 \(y\) ，设为\(x'\)和\(y'\)

2. 第二层中：\(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'\) 的关系为

3. 第二层的 \(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\))。要求通解也很简单，如下

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

\(\mathrm{CRT}\) 问题的解决方法

构造出解

  \(\mathrm{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]\)：

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

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

逆元求法

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

  1、费马小定理，限制模数必须为质数。 \(\mathrm{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\)，方程有整数解） ,代入扩展欧几里得模板就可以求出来了。

\(\mathrm{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]猜数字

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

\(\mathrm{EXCRT}\) 问题的解决方法

为啥不能继续用 \(\mathrm{CRT}\) 的代码解决问题了呢？

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

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

合并同余方程组

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

  我们可以发现，同余方程在一定条件下是可以合并的，但是也会出现无解的情况。合并后的方程仍然可以表示为同余方程的形式，并且模数是原来模数的最小公倍数（上述这些规律其实是要证明的，但是蒟蒻不太会，感兴趣的同学可以自己研究一下上面大佬的博客和其他资料）。那么，解决 \(\mathrm{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_2-r_1\) 不能是负数（否则快速乘会 \(TLE\) ) , 所以这里都要使用一些取模的技巧。之后要更新 \(ans\) 和 \(M\) ,让 \(M\) 成为前 \(k\) 个同余方程模数的最下公倍数。这里有个公式是 \(gcd(a,b)*lcm(a,b)=a*b\) ,可以用它来更新 \(M\) 。

\(\mathrm{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 一个小问题

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

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