密码学基础算法(二)中国剩余定理

随便谷歌了一个图片做首图

原图地址: http://www.siwapu.com/etagid41968b0/

密码学基础系列:

• (一) 基于整数的欧几里得算法和扩展欧几里得算法
• (二) 中国剩余定理
• (三) 扩展域 $GF(2^n)$ 上的多项式运算
• (四) 扩展域 $GF(2^n)$ 上的欧几里得算法和扩展欧几里得算法

看书自学的一个比较大的问题就是，一旦对书中的问题有疑问，尤其是一些小的细节，很难找到人答疑。

学习中国剩余定理时(CRT: Chinese Reminder Theorem)，一直有点迷糊，没有完全明白中国剩余定理到底如何用。
陆续看了一些资料，

• 《深入浅出密码学》中没有对中国剩余定理的证明，只在 7.5.2 节中提到使用中国剩余定理快速加密
• 《密码编码学与网络安全》各版本(包括六/七/八版)中，对中国剩余定理的描述相当复杂(也可能是我太蠢了)
• 《数论概论》第三版第 11 章中，只说了中国剩余定理的表述，没有提供完整计算方式
• 《密码学原理与实践》第三版中，第 5.2.2 节对中国剩余定理有详细的证明和描述，是我见过最易懂的说明。

这里摘录一下《密码学原理与实践》中关于中国剩余定理的描述(第5.2.2节)：

(中国剩余定理) 假定$m_1,m_2,...,m_r$为两两互素的正整数，又假定$a_1,a_2,...,a_r$为整数。那么同余方程组:
$x\equiv a_i(mod{m}_i)(1\le i\le r)$ 有模 $M=m_1\times m_2\times ...\times m_r$的唯一解，此解由下式给出:
$$x=\sum _{i=1}^r{a}_i{M}_i{y}_{i}modM$$

其中, $M_i=M/m_i$，且$y_i=M_i^{-1}mod{m}_i,1\le i\le r$

从定理出发，明确了定义之后，就容易计算了。

这里以《孙子算经》中的原文为例:

今有物不知其数，三三数之余二，五五数之余三，七七数之余二，问物几何?

1. 问题

这个问题翻译过来就是:

已知 3 个同余方程组:

• $2\equiv xmod3,3\equiv xmod5,2\equiv xmod7$,

求$x$。

从条件得知:

1. $(m_1,m_2,m_3)\Rightarrow (3,5,7)$，即$m_1=3,m_2=5,m_3=7$
2. $(a_1,a_2,a_3)\Rightarrow (2,3,2)$，即$a_1=2,a_2=3,a_3=2$

2. 求解

第一步，计算同余方程组的模$M=m1\times m2\times ...\times m_r$

• $M=m_1\times m_2\times m_3=3\times 5\times 7=105$

第二步，求$M_i$及其基于$m_i$的逆元$M_i^{-1}$和$y_i$

这里的逆元$M_i^{-1}$使用《密码学基础算法(一)基于整数的欧几里得算法和扩展欧几里得算法》一文中求逆的函数int int_inv(int a, int b)来求取。

• $M_1=M/m_1=105/3=35$, $M_1^{-1}=int\_inv(M_1,m_1)=int\_inv(35,3)=2$, $y_1=M_1^{-1}mod{m}_1=2mod3=2$
• $M_2=M/m_2=105/5=21$, $M_2^{-1}=int\_inv(M_2,m_2)=int\_inv(21,5)=1$, $y_2=M_2^{-1}mod{m}_2=1mod5=1$
• $M_3=M/m_3=105/7=15$, $M_3^{-1}=int\_inv(M_3,m_3)=int\_inv(15,7)=1$, $y_3=M_3^{-1}mod{m}_3=1mod7=1$

第三步，根据第二步计算的$M/M_i/y_i$，使用中国剩余定理的公式求$x$

$$\begin{gathered} x=\sum _{i=1}^r{a}_i{M}_i{y}_{i}modM=(a_1{M}_1{y}_1+a_2{M}_2{y}_2+a_3{M}_3{y}_3)mod105 \\ =(2\times 35\times 2+3\times 21\times 1+2\times 15\times 1)mod105 \\ =(140+63+30)mod105 \\ =233mod105 \\ =23 \end{gathered}$$

所以最后得到要求的数为 23。

3. 代码

我们可以根据前面的计算过程将其翻译成代码

/**
 * @description: 中国剩余定理(CRT: Chinese Reminder Theorem)求同余方程解
 * @param {int} 指针 m, 指向整型数组, 包含每一组除数
 * @param {int} 指针 r, 指向整型数组, 包含每一组余数
 * @param {int} 整数 n, 指定传入指针 m 和 r 指向数组的大小
 * @return {*} 中国剩余定理的同余方程解
 */
int crt(int *m, int *r, int n)
{
    int i;
    int M, Mi, yi;
    int sum;

    // 计算同于方程组的模 M
    M = 1;
    for (i=0; i<n; i++)
    {
        M *= m[i];
    }

    sum = 0;
    for (i=0; i<n; i++)
    {
        // 计算每一个 Mi
        Mi = M / m[i];

        // 计算每一个 yi
        yi = int_inv(Mi, m[i]) % m[i];

        // 累积求和
        sum += r[i] * Mi * yi;
    }
    
    return sum % M;
}

其中调用的函数int_inv来自《密码学基础算法(一)基于整数的欧几里得算法和扩展欧几里得算法》，专门用于求整数的乘法逆元，如下:

/**
 * @description: 基于扩展欧几里得算法计算整数 a 关于 整数 b 的乘法逆元(展开 int_gcd_ex 后的优化版本)
 * @param {int} a 整数 a
 * @param {int} b 整数 b
 * @return {*} 返回整数 a 关于整数 b 的最小正整数乘法逆元, 乘法逆元不存在(gcd(a, b) != 1)则返回 0;
 */
int int_inv(int a, int b)
{
    int x, x0, x1;
    int q, r, t;

    /* 消除警告: "warning: ‘x’ may be used uninitialized in this function" */
    x = 0;

    /* 初始化最初两项系数 */
    x0 = 1;
    x1 = 0;

    /* 临时保存 b */
    t = b;

    q = a / b;
    r = a % b;

    while (r != 0)
    {
        /* 计算当前项 x */
        x = x0 - q * x1;

        /* 依次保存前两项到 x0, x1 */
        x0 = x1; x1 = x;

        a = b;
        b = r; /* 前一次的余数 */

        q = a / b;
        r = a % b;
    }

    /* gcd(a, b) != 1, 没有乘法逆元 */
    if (b != 1)
    {
        return 0;
    }

    /* 调整小于 0 的情况 */
    if (x < 0)
    {
        x += t;
    }

    return x;
}

4. 测试代码和示例

• 测试代码

int main(int argc, char *argv[])
{
    int m1[3] = { 3,  5,  7}, r1[3] = { 2,  3,  2}; // 示例1: 孙子算经
    int m2[3] = { 7, 11, 13}, r2[3] = { 5,  3, 10}; // 示例2: 《密码学原理与实践》第三版, p133, 例 5.3
    int m3[3] = {25, 26, 27}, r3[3] = {12,  9, 23}; // 示例3: 《密码学原理与实践》第三版, p178, 题 5.6
    int m4[2] = {37, 49},     r4[2] = {11, 42};     // 示例4: 《密码编码学与网络安全》第七版, p37, 示例
    int res;

    res = crt(m1, r1, 3);
    printf("同余方程组: 除数{ 7, 11, 13}, 余数{ 5,  3, 10}的解为 %d\n", res);

    res = crt(m2, r2, 3);
    printf("同余方程组: 除数{ 3,  5,  7}, 余数{ 2,  3,  2}的解为 %d\n", res);

    res = crt(m3, r3, 3);
    printf("同余方程组: 除数{25, 26, 27}, 余数{12,  9, 23}的解为 %d\n", res);

    res = crt(m4, r4, 2);
    printf("同余方程组: 除数{37, 49},     余数{11, 42}    的解为 %d\n", res);

    return 0;
}

• 编译和运行

$ gcc crt.c -o crt
$ ./crt
同余方程组: 除数{ 7, 11, 13}, 余数{ 5,  3, 10}的解为 23
同余方程组: 除数{ 3,  5,  7}, 余数{ 2,  3,  2}的解为 894
同余方程组: 除数{25, 26, 27}, 余数{12,  9, 23}的解为 14387
同余方程组: 除数{37, 49},     余数{11, 42}    的解为 973

5. 中国剩余定理的应用

根据《密码编码学与网络安全》中的说明，运用中国剩余定理，使得模 M 的大数运算可以转化为相对较小的数的运算。当 M 为 150 位或 150 位以上时，这种方法特别有效。

例如:

973 对 37 和 49 取模分别为: $973mod37=11,973mod49=42$，因此把 973 表示为(11, 42)的形式。

假如要计算 678 + 973 = 1651, 用处理 973 的办法同样处理 678: $678mod37=12,678mod49=41$，因此把 678 表示为(12, 41)

然后将上面得到的两个二元组的元素相加并化简:
$1651=973+678=((11+12)mod37,(42+41)mod49)=(23,34)$

从上面的例子可以看到，大数 1651 (其实很小) 对 37 和 49 取模转换成了小数 973 和 678 对 37 和 49 取模。

虽然看起来计算量差不多，但是对很大的数，前面提到的 150 位或更多这种，转换成两个数乘积的形式，新的数短一半，再进行操作是很有优势的。

6. 完整代码

#include <stdio.h>

/**
 * @description: 基于扩展欧几里得算法计算整数 a 关于 整数 b 的乘法逆元(展开 int_gcd_ex 后的优化版本)
 * @param {int} a 整数 a
 * @param {int} b 整数 b
 * @return {*} 返回整数 a 关于整数 b 的最小正整数乘法逆元, 乘法逆元不存在(gcd(a, b) != 1)则返回 0;
 */
int int_inv(int a, int b)
{
    int x, x0, x1;
    int q, r, t;

    /* 消除警告: "warning: ‘x’ may be used uninitialized in this function" */
    x = 0;

    /* 初始化最初两项系数 */
    x0 = 1;
    x1 = 0;

    /* 临时保存 b */
    t = b;

    q = a / b;
    r = a % b;

    while (r != 0)
    {
        /* 计算当前项 x */
        x = x0 - q * x1;

        /* 依次保存前两项到 x0, x1 */
        x0 = x1; x1 = x;

        a = b;
        b = r; /* 前一次的余数 */

        q = a / b;
        r = a % b;
    }

    /* gcd(a, b) != 1, 没有乘法逆元 */
    if (b != 1)
    {
        return 0;
    }

    /* 调整小于 0 的情况 */
    if (x < 0)
    {
        x += t;
    }

    return x;
}

/**
 * @description: 中国剩余定理(CRT: Chinese Reminder Theorem)求同余方程解
 * @param {int} 指针 m, 指向整型数组, 包含每一组除数
 * @param {int} 指针 r, 指向整型数组, 包含每一组余数
 * @param {int} 整数 n, 指定传入指针 m 和 r 指向数组的大小
 * @return {*} 中国剩余定理的同余方程解
 */
int crt(int *m, int *r, int n)
{
    int i;
    int M, Mi, yi;
    int sum;

    // 计算同于方程组的模 M
    M = 1;
    for (i=0; i<n; i++)
    {
        M *= m[i];
    }

    sum = 0;
    for (i=0; i<n; i++)
    {
        // 计算每一个 Mi
        Mi = M / m[i];

        // 计算每一个 yi
        yi = int_inv(Mi, m[i]) % m[i];

        // 累积求和
        sum += r[i] * Mi * yi;
    }

    return sum % M;
}

int main(int argc, char *argv[])
{
    int m1[3] = { 3,  5,  7}, r1[3] = { 2,  3,  2}; // 示例1: 孙子算经
    int m2[3] = { 7, 11, 13}, r2[3] = { 5,  3, 10}; // 示例2: 《密码学原理与实践》第三版, p133, 例 5.3
    int m3[3] = {25, 26, 27}, r3[3] = {12,  9, 23}; // 示例3: 《密码学原理与实践》第三版, p178, 题 5.6
    int m4[2] = {37, 49},     r4[2] = {11, 42};     // 示例4: 《密码编码学与网络安全》第七版, p37, 示例
    int res;

    res = crt(m1, r1, 3);
    printf("同余方程组: 除数{ 7, 11, 13}, 余数{ 5,  3, 10}的解为 %d\n", res);

    res = crt(m2, r2, 3);
    printf("同余方程组: 除数{ 3,  5,  7}, 余数{ 2,  3,  2}的解为 %d\n", res);

    res = crt(m3, r3, 3);
    printf("同余方程组: 除数{25, 26, 27}, 余数{12,  9, 23}的解为 %d\n", res);

    res = crt(m4, r4, 2);
    printf("同余方程组: 除数{37, 49},     余数{11, 42}    的解为 %d\n", res);

    return 0;
}

7. 其它

洛奇工作中常常会遇到自己不熟悉的问题，这些问题可能并不难，但因为不了解，找不到人帮忙而瞎折腾，往往导致浪费几天甚至更久的时间。

所以我组建了几个微信讨论群(记得微信我说加哪个群，如何加微信见后面)，欢迎一起讨论:

• 一个密码编码学讨论组，主要讨论各种加解密，签名校验等算法，请说明加密码学讨论群。
• 一个Android OTA的讨论组，请说明加Android OTA群。
• 一个git和repo的讨论组，请说明加git和repo群。

在工作之余，洛奇尽量写一些对大家有用的东西，如果洛奇的这篇文章让您有所收获，解决了您一直以来未能解决的问题，不妨赞赏一下洛奇，这也是对洛奇付出的最大鼓励。扫下面的二维码赞赏洛奇，金额随意：

洛奇自己维护了一个公众号“洛奇看世界”，一个很佛系的公众号，不定期瞎逼逼。公号也提供个人联系方式，一些资源，说不定会有意外的收获，详细内容见公号提示。扫下方二维码关注公众号：