中國剩余定理(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. \]

  在中國剩余定理（以下簡稱 \(\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'\) 的關系為

\[x=y^{\prime} ,\quad y=x^{\prime}-\left\lfloor\frac{a}{b}\right\rfloor y^{\prime} \]

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\))。要求通解也很簡單，如下

\[\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)}\) ，證明當然也沒有啦，可以手推一下的。

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

\[\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\) , 而 \(\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\) 不一定互質，從而導致了

\[\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. \]

這一條件失效，逆元也求不出來。所以我們不能再用逆元來解決 \(\mathrm{EXCRT}\) 的問題，必須換一種思路。

合並同余方程組

  觀察一下幾個簡單的式子（摘自一位讓我學會 \(\mathrm{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.\]

  我們可以發現，同余方程在一定條件下是可以合並的，但是也會出現無解的情況。合並后的方程仍然可以表示為同余方程的形式，並且模數是原來模數的最小公倍數（上述這些規律其實是要證明的，但是蒟蒻不太會，感興趣的同學可以自己研究一下上面大佬的博客和其他資料）。那么，解決 \(\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使用過度，頁面渲染好像出了點問題😂)