《夜深人静写算法》数论篇 - (20) 中国剩余定理

前言

    当初学这个算法，纯粹是因为她叫 —— 中国剩余定理它的另一个名字 —— 孙子定理。有兴趣的朋友可以自行百度它原本的名字。

一、 物不知数

1、引例

【例题1】有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二。问物几何？

    这个问题是中国南北朝时期的数学著作《孙子算经》卷下第二十六题，叫做 “物不知数” 问题，这个问题翻译成现代文的意思是：一个整数除以三余二，除以五余三，除以七余二，求这个整数。表示如下:
$$(S):\{\begin{array}{l}x\equiv 2(mod3)\\ x\equiv 3(mod5)\\ x\equiv 2(mod7)\end{array}$$

    由于每一行都是一个同余方程，其中 $x$ 是公共的未知数。
    用这样的形式表示的我们称它为同余方程组，$S$ 就是这个同余方程组。

二、同余方程组

    我们将同余方程表示成更加一般的形式，如下：
$$(S):\{\begin{array}{ll}x& \equiv a_1(mod{m}_1)\\ x& \equiv a_2(mod{m}_2)\\ & ⋮\\ x& \equiv a_n(mod{m}_n)\end{array}$$

    我们知道同余方程是可以表示成等式的形式的，那么同余方程组自然也可以，我们可以这么来表示：
$$\begin{gathered} (S):\{\begin{array}{ll}x& =k_1{m}_1+a_1\\ x& =k_2{m}_2+a_2\\ & ⋮\\ x& =k_n{m}_n+a_n\end{array} \\ (k_i\in Z) \end{gathered}$$

    $n$ 个方程，$n+1$ 个未知数的情况下，解的情况可能存在，也可能不存在，存在的情只可能是无穷多解。

1、不存在解

    这种情况比较好理解，当存在 $i$ 和 $j$，同时满足 $m_i=m_j$ 且 $a_i\ne a_j$ 时，解必然是不存在的。
    因为任何一个数模另外一个数，得到的结果一定是确定的，不可能有两种结果。

2、存在无穷多解

    一旦存在解，那么对于任意两个 $i$ 和 $j$，只要 $m_i=m_j$ ，必然满足 $a_i=a_j$，出现这种情况时，变量数和方程数同时减少。
    换言之，未知数永远比方程数多 1，这样就一定有无穷多解了。

三、朴素算法

1、模数互素

    朴素算法就是穷举，以 “物不知数” 问题为例，我们可以从零开始枚举 $x$ 的值。然后不断试除，看什么时候能够满足三个同余方程同时成立。
$$(S):\{\begin{array}{l}x\equiv 2(mod3)\\ x\equiv 3(mod5)\\ x\equiv 2(mod7)\end{array}$$

    代码实现如下：

#include <iostream>
using namespace std;

const int MAXN = 3;
int m[MAXN] = {3, 5, 7};
int a[MAXN] = {2, 3, 2};

int main() {
    for(int x = 0; x < 400; ++x) {      // (1)
        bool bfind = true;
        for(int i = 0; i < MAXN; ++i) { 
            if( x % m[i] != a[i] ) {    // (2)
                bfind = false;
                break;
            }
        }
        if(bfind) {
            printf("%d\n", x);          // (3)
        }
    }	
    return 0;
} 

• $(1)$ 由于解有无数个，所以为了输出的时候能够观测到一些规律，我们只输出小于 400 的解。
• $(2)$ 对于每个 $(m_i,a_i)$ 的整数对，用枚举得到的 $x$ 进行试除，看能否整除，一旦不满足就跳出循环；
• $(3)$ 输出满足所有等式情况的解。

    输出结果为：

23
128
233
338

    这时我们发现，它其实是一个等差数列，并且公差是105，正好是3 * 5 * 7，也就是所有模数的乘积。那么，如果令最小的非负整数解为 $x^{\prime }$，我们猜测这个同余方程组的通解为：
$$\begin{gathered} x=x^{\prime }+k\prod _{i=1}^n{m}_i \\ (k=0,1,2,...) \end{gathered}$$

2、模数不互素

    事实上，上面这个公式只适合 模数互素 的情况，当模数不互素时，我们来看一个例子，如下：
$$(S):\{\begin{array}{l}x\equiv 2(mod4)\\ x\equiv 4(mod6)\\ x\equiv 7(mod15)\end{array}$$

    只需要将上面的代码稍微改一下，数据部分改成：

int m[MAXN] = {4, 6, 15};
int a[MAXN] = {2, 4, 7};

• 运行后得到输出结果如下：

22
82
142
202
262
322
382

    我们通过观察发现，它仍然是一个等差数列，并且公差是60，正好是 $lcm(4,6,15)$，也就是所有模数的最小公倍数。我们再联想刚才所有数互素的情况，它正好是所有数的乘积。
    那么，如果令最小的非负整数解为 $x^{\prime }$，我们猜测这个同余方程组的通解为：
$$\begin{gathered} x=x^{\prime }+k\times lcm(m_1,m_2,...,m_n) \\ (k=0,1,2,...) \end{gathered}$$

四、中国剩余定理

1、算法详解

    根据上述的猜测，我们能够发现 模数不互素 的情况，是完全包含 模数互素 的情况的。
    所以这里我直接介绍 模数不互素 的情况下，同余方程组通解的求解方式，网上有人叫它扩展中国剩余定理，我觉得没有 中国剩余定理 霸气，但是这不重要，知道算法原理就行啦！
    我们给同余方程组进行编号，如下：
$$(S):\{\begin{array}{llc}x& =k_1{m}_1+a_1& (1)\\ x& =k_2{m}_2+a_2& (2)\\ & ⋮& \\ x& =k_n{m}_n+a_n& (n)\end{array}$$

    然后尝试把第 $(1)$ 个方程和第 $(2)$ 个方程进行合并，得到：
$$k_1{m}_1+a_1=k_2{m}_2+a_2$$

    移项后得到：
$$k_1{m}_1-k_2{m}_2=a_2-a_1$$

    这是一个类似 $Ax+By=C$ 的同余方程。
    其中 $A$、$B$、$C$ 都是常数，$X$ 和 $Y$ 为未知数，所以我们可以得到如下的变量转换关系：
$$\{\begin{array}{l}A=m_1\\ B=-m_2\\ C=a_2-a_1\\ X=k_1\\ Y=k_2\end{array}$$

    如果 $Cmodgcd(A,B)\ne 0$ ，则方程无解，即整个方程组无解。
    否则，令 $A^{\prime }=\frac{A}{gcd(A,B)}$，$B^{\prime }=\frac{B}{gcd(A,B)}$，$C^{\prime }=\frac{C}{gcd(A,B)}$，得到新的等式：
$$A^{\prime }x+B^{\prime }y=C^{\prime }$$

    通过扩展欧几里德定理，必然能够求出一个最小的正整数解 $Y_0$，于是有 $Y$ 也就是 $k_2$ 的通解如下：
$$Y=Y_0+A^{\prime }t$$

    将 $Y$ 代入等式 $(2)$，得到：
$$\begin{array}{rl}x& =k_2{m}_2+a_2\\ & =Y{m}_2+a_2\\ & =(Y_0+A^{\prime }t)m_2+a_2\\ & =Y_0{m}_2+A^{\prime }{m}_{2}t+a_2\\ & =A^{\prime }{m}_{2}t+(Y_0{m}_2+a_2)\\ & =\frac{m_1{m}_2}{gcd(m_1,m_2)}t+(Y_0{m}_2+a_2)\\ & =lcm(m_1,m_2)t+(Y_0{m}_2+a_2)\end{array}$$

    其中 $lcm(m_1{m}_2)$ 和 $(Y_0{m}_2+a_2)$ 均为常数，$t$ 为未知数，于是我们得到了一个新的关于 $x$ 的方程，总的方程组数量减少了 1，用同样的方法，经过 $n-1$ 次迭代，必然可以将 $n$ 个方程转化成一个方程， 如下：
$$x=x^{\prime }+k\times lcm(m_1,m_2,...,m_n)$$

    而这个方程，正式我们上一章节猜测的那个方程，当 $k=0$ 时，$x^{\prime }$ 就是最小非负整数解了。

2、算法综述

算法流程如下：
    1）将所有方程组标准化，即模数和余数都变成正数；
    2）利用【1、算法详解】中的方法合并两个方程，使得方程数减一；合并过程中可能发生两种情况：
    2.a）无解的情况，直接返回；
    2.b）存在解，进入 3）；
    3）利用扩展欧几里得定理求出两个方程的通解，从而变成一个新的同余方程，如果只剩下最后一个同余方程了，那么它的余数就是同余方程组的一个最小非负整数解，否则回到 2）；

    通解如下：
$$x=x_{min}+k\times lcm(m_1,m_2,...,m_n)$$

    $x_{min}$ 即同余方程组的最小非负整数解，$k=0,1,2,3,...$；

3、算法实现

1）类设计

    首先，设计一个类，这个类只有两个成员m_(代表模数)和a_(代表余数)，用来代表上文提到的 方程：$x=km+a$，然后就可以用数组来代表方程组了。

#define ll long long

class ModPair {
private:
    ll m_;   // modulus
    ll a_;   // remainder
public:
    ModPair() {}
    ModPair(ll m, ll a) : m_(m), a_(a) {}
    ll getModulus() const { return m_;}
    void setModulus(ll m) {m_ = m;}
    ll getRemainder() const { return a_;}
    void setRemainder(ll a) {a_ = a;} 
    
    void standardization();               // (1)
    ll lcm (ll om) const;                 // (2)
};

    除了一些基本的set和get接口，额外提供两个接口：

• $(1)$ standardization()用来做标准化，即把m_和a_都标准化为正数，如下：

void ModPair::standardization() {
    if(m_ < 0) {
        m_ = -m_;
    }
    a_ = (a_ % m_ + m_) % m_;
}

    $(2)$ lcm (ll om)返回模数和传参的最小公倍数，实现如下：

ll ModPair::lcm (ll om) const {
    ll g = GCD(om, m_);
    return om / g * m_;
}

2）求解过程

    实际求解过程，先给出代码，通过需要来逐一解释：

ModPair chineseRemain(int n, ModPair* mp) {
    /********************(1)********************/
    for(int i = 0; i < n; ++i) {
        mp[i].standardization();
    } 
    ModPair ans = mp[0];
    ll A, B, C, X, Y;
    for(int i = 1; i < n; i++) {
        /********************(2)********************/
        A = ans.getModulus();         
        B = -mp[i].getModulus();
        C = mp[i].getRemainder() - ans.getRemainder();
        ll g = GCD(A, B);
        if( C % g ) {     
            ans.setRemainder(-1);             
            return ans;
        }
        /********************(3)********************/
        A /= g, B /= g, C /= g;
        /********************(4)********************/
        if(A < 0) {
            A = -A, B = -B, C = -C;
            B = ((B % A) + A) % A;
        }
        /********************(5)********************/
        ExpGcd(A, B, X, Y);
        Y = Product_Mod(Y, C, A);
        /********************(6)********************/
        ll tmpm = A * mp[i].getModulus();
        ll tmpr = (mp[i].getRemainder() + Product_Mod(mp[i].getModulus(), Y, tmpm)) % tmpm;
        
        ans.setModulus(tmpm);
        ans.setRemainder(tmpr);
    }
    return ans;
}

• $(1)$ 可以这么认为，任何一个方程 $x=km+a$，$k$ 可以是任意整数，所以 $m$ 为正为负是无关紧要的，为了标准，将 $m$ 转成整数；同样，$a=xmodm$，自然取值要落在 $[0,m)$ 才合理。于是，这一步就是对所有的方程都做了一次标准化，方便后续计算。
• 并且将第一个方程的参数的两个参数赋值给 ans。
• $(2)$ 这一步做的事情，就是抽象出一些变量 $A$、$B$、$C$，方便计算。并且判断 $gcd(A,B)$ 是否能整除 $C$，如果不能够整除，则方程组无解，直接将 ans的模数参数定为-1返回，调用方根据这个标记为去判断接下来要做的事情；
• $(3)$ 让等式两边同时除以 $gcd(A,B)$，等式是不变的，即 $A^{\prime }=\frac{A}{gcd(A,B)}$，$B^{\prime }=\frac{B}{gcd(A,B)}$，$C^{\prime }=\frac{C}{gcd(A,B)}$；
• $(4)$ 将 $A$ 和 $B$ 进行标准化，转化成正数；
• $(5)$ 利用扩展欧几里得，求出方程的最小正整数解，这里的Product_Mod是一个乘法求幂的接口，能够确保 64位 整数相乘但不溢出。
• $(6)$ 最后一步，就是得到的通解回代回ans的过程。
• 通过这样的方法，最后得到的ans的两个参数分别为所有模数的最小公倍数，以及方程组的最小非负整数解。