中国剩余定理（孙子定理）的算法实现（基于miracl大数运算库）中国剩余定理（孙子定理）的算法实现：

53 阅读 0 评论 35 点赞

我是靠谱客的博主机智豌豆，这篇文章主要介绍中国剩余定理（孙子定理）的算法实现（基于miracl大数运算库）中国剩余定理（孙子定理）的算法实现：，现在分享给大家，希望可以做个参考。

中国剩余定理（孙子定理）的算法实现：

一、实现目标：

根据中国剩余定理，设正整数两两互素，那么对于任意k个整数，同余方程组：

必有解，模的解数为1。方程组元素的传入是通过文本文件读入的，顺序是，，每个数字之间是通过换行来分割的，数字大小最大设值为500位。判断正整数是否两两互素；是，则通过中国剩余定理算出同余方程组的解；否则跳出，输出“不能直接利用中国剩余定理”。

二、方案设计

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

2.1关于中国剩余定理的数学定义

设正整数m1,m2,⋯,mk两两互素，那么对于任意k个整数a1,a2,⋯,ak，同余方程组

必有解，模i=1kmi的解数为1。该解是x≡M1M1-1a1+⋯+MkMk-1ak(modm)

其中m=j=1kmj，∀1≤j≤k,Mj=mmj，Mj-1是满足MjMj-1≡1(modmj)的一个整数。

2.2 中国剩余定理的算法步骤

根据中国剩余定理，要实现其算法实现，算法步骤可以大致分为以下的四步：

Step1：判断正整数m1,m2,⋯,mk是否两两互素：是，则继续，转到Step2；否则跳出，输出“不能直接利用中国剩余定理”。

Step2：计算

Step3：计算

Step4：计算

三、方案实现

3.1 主要函数的介绍

3.1.1 multiply

函数原型: void multiply(big x, big y, big z);

功能说明: 两个大数相乘，z=x*y。

3.1.2 fdiv

函数原型：void fdiv(x,y,z)；

功能说明：将两个大数相除，z=x/y。

3.1.3 xgcd

函数原型: int xgcd(bigx, big y, big xd, big yd, big z);

功能说明: 计算两个大数的扩展最大公约数，也可以用来计算模逆，这个函数比mad 函数运算速度稍慢。z=gcd(x,y)=x*xd+y*yd

例子: xgcd(x,p,x,x,x); //计算x^-1 mod p

/* x = 1/x mod p (p is prime) */

在miracl库中，要求得一个数它的逆元的话，有一个函数invers()可以求得逆元，但是invers()只能够求得int型数据的逆元。所以在miracl库中搜索可以求得逆元的函数时就找到了这个函数，跟着这个函数所举出的例子，在函数中确实可以求出逆元，但是一开始时我对于这个函数的运行机理产生了疑惑。就拿这个例子来说，按照函数说明，第二个x是第一个x的逆元，第三个x是p的逆元，最后一个x是第一个x和p的最大公因数（在此处应该是1才对），既然有了这么多x，函数是怎么知道应该用哪个x呢？在函数的操作手册中我找到了答案。其限制条件是这么说的：If xd and yd are not distinct, only xd is returned. The GCD is only returned if z distinct from both xd and yd.也就是说，例子中的传入参数xd、yd以及z都是一样的，所以说函数只会返回x在模p下的逆元。

3.1.4 add

函数原型: void add(big x, big y, big z);

功能说明: 两个大数相加，z=x+y。

Example: add(x,x,x); // This doubles the value of x.

3.1.5 powmod

函数原型: void powmod(big x, big y,big z, big w);

功能说明: 模幂运算，w=xy(modz)。

3.2 实验代码

复制代码

#include "miracl.h"
#define z 50   //用作下面的big型数组初始化，这里宏定义为数组里面有五十个big型的数据，如果从文件中导入的数据超过了五十个，可以修改；
int main() {
    FILE *fp;
    big a[z], m[z], x[z], t[z];
    big Xul, Mul, CON, Mi, Mj, Mm, re;
    int i, k, b, j;
    miracl *mip = mirsys(5000, 10);   //初始化系统，分配五千位的十进制空间
    mip->IOBASE = 10;    //以十进制的形式输入输出数据
    Xul = mirvar(0);
    Mul = mirvar(1);
    CON = mirvar(0);
    Mi = mirvar(0);
    Mj = mirvar(0);
    Mm = mirvar(0);
    re = mirvar(0);
    i = 0;
    k = 0;
    b = 0;
    j = 0;
    if ((fp = fopen("2.txt", "r")) == NULL)  //打开文件
    {
       printf("Open File Failue.....n");
       return -1;
    }
    //从文件中读入所有的数据放入数组t[]中，也就将所有的a和m先读入进来，下面的操作进行分开，这里b的值加上1就是数组t[]中的元素个数（在实验中，元素个数一定为奇数）
    for (b = 0; !feof(fp); b++)
    {
       t[b] = mirvar(0);    //数组中的每一个数据在被读入之前，都必须初始化
       cinnum(t[b], fp);    //读入数据
    }
   
    k = (b + 1) / 2;   //k的值大小代表着 数组中元素个数的一半
    i = 0;//在VC6.0中，不能有for (int i = 0; i < k; i++)这样的形式，负责计数的i，必须在for循环开始之前就已经被声明并且赋上初始值，后面的每个for循环前都会有这样的操作
    for (i = 0; i < k; i++) {
       a[i] = mirvar(0);
       a[i] = t[i];
    }//在实验中传入的文本文件的格式为前一半的数据为a1，a2，a3，....这样的，将矩阵t[]中前一半的元素传入数组a[]即可得到操作数a所在的数组
   
    i = k;
    for (i = k; i <=b; i++) {
       m[i - k] = mirvar(0);
       m[i - k] = t[i];
    }//将数组t[]后半部分的数据m1，m2，m3，....传入数组m[]得到操作数m所在的数组
   
    i = 0;
    for (i = 0; i < k; i++) {
       x[i] = mirvar(1);
    }
    //以上的工作是初始化，并且将数据从文本文件中转移到相对应的数组中
    i = 0;
    j = 0;
    for (i = 0; i < k; i++)   //两两比较mi和mj，以达到中国剩余定理的条件：两两互素
    {
       for (j = 0; j < k; j++)
       {
           if (i == j)
           {
              continue;  //自己和自己不需要比较
           }
           else
           {
              egcd(m[i], m[j], CON);    //计算两个mi和mj的最大公因数，他们的最大公因数赋值给CON
              if (compare(CON, mirvar(1))) {
                  break;      //最大公因数不为1的话，不满足条件
              }
           }
       }
       if (compare(CON, mirvar(1)))
       {
           break;
       }
    }
    if (compare(CON, mirvar(1)))
    {
       printf("Input positive integer m（s）don not match the requirement of Chinese remainder theorem!");     //上述部分判断是否满足中国剩余定理的条件
    }
    else    //满足中国剩余定理条件
    {
       i = 0;
       for (i = 0; i < k; i++)    //将所有的m乘起来，这里的Mul相当于书上的m（所有的mj的乘积）
       {
           multiply(Mul, m[i], Mul);   //multiply函数作用：Mul=Mul*m[i]
       }
       i = 0;
       for (i = 0; i < k; i++)
       {
           fdiv(Mul, m[i], Mi);//fdiv函数作用：Mi = Mul/ m[i];
           xgcd(Mi, m[i], Mj, Mj, Mj);//Mj = invers(Mi, m[i]);   计算两个大数的扩展最大公约数，也可以用来计算模逆，z=gcd(x,y)=x.xd+y.yd，在本例中是求得Mi的模逆为Mj
           multiply(Mi, Mj, Mm);  //Mm=Mi*Mj
           multiply(Mm, a[i], x[i]);//x[i] = (Mi*Mj*a[i]);
       }
       i = 0;
       for (i = 0; i < k; i++)
       {
           add(Xul, x[i], Xul);//X = X + x[i];
       }
       powmod(Xul, mirvar(1), Mul, re);//re = X % M;
       printf("the resault of x is :n");
       cotnum(re, stdout);
    }
    mirkill(Xul);
    mirkill(Mul);
    mirkill(CON);
    mirkill(Mi);
    mirkill(Mj);
    mirkill(Mm);
    mirkill(re);
    return 0;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
#include "miracl.h"
#define z 50   //用作下面的big型数组初始化，这里宏定义为数组里面有五十个big型的数据，如果从文件中导入的数据超过了五十个，可以修改；
int main() {
    FILE *fp;
    big a[z], m[z], x[z], t[z];
    big Xul, Mul, CON, Mi, Mj, Mm, re;
    int i, k, b, j;
    miracl *mip = mirsys(5000, 10);   //初始化系统，分配五千位的十进制空间
    mip->IOBASE = 10;    //以十进制的形式输入输出数据
    Xul = mirvar(0);
    Mul = mirvar(1);
    CON = mirvar(0);
    Mi = mirvar(0);
    Mj = mirvar(0);
    Mm = mirvar(0);
    re = mirvar(0);
    i = 0;
    k = 0;
    b = 0;
    j = 0;
    if ((fp = fopen("2.txt", "r")) == NULL)  //打开文件
    {
       printf("Open File Failue.....n");
       return -1;
    }
    //从文件中读入所有的数据放入数组t[]中，也就将所有的a和m先读入进来，下面的操作进行分开，这里b的值加上1就是数组t[]中的元素个数（在实验中，元素个数一定为奇数）
    for (b = 0; !feof(fp); b++)
    {
       t[b] = mirvar(0);    //数组中的每一个数据在被读入之前，都必须初始化
       cinnum(t[b], fp);    //读入数据
    }
   
    k = (b + 1) / 2;   //k的值大小代表着 数组中元素个数的一半
    i = 0;//在VC6.0中，不能有for (int i = 0; i < k; i++)这样的形式，负责计数的i，必须在for循环开始之前就已经被声明并且赋上初始值，后面的每个for循环前都会有这样的操作
    for (i = 0; i < k; i++) {
       a[i] = mirvar(0);
       a[i] = t[i];
    }//在实验中传入的文本文件的格式为前一半的数据为a1，a2，a3，....这样的，将矩阵t[]中前一半的元素传入数组a[]即可得到操作数a所在的数组
   
    i = k;
    for (i = k; i <=b; i++) {
       m[i - k] = mirvar(0);
       m[i - k] = t[i];
    }//将数组t[]后半部分的数据m1，m2，m3，....传入数组m[]得到操作数m所在的数组
   
    i = 0;
    for (i = 0; i < k; i++) {
       x[i] = mirvar(1);
    }
    //以上的工作是初始化，并且将数据从文本文件中转移到相对应的数组中
    i = 0;
    j = 0;
    for (i = 0; i < k; i++)   //两两比较mi和mj，以达到中国剩余定理的条件：两两互素
    {
       for (j = 0; j < k; j++)
       {
           if (i == j)
           {
              continue;  //自己和自己不需要比较
           }
           else
           {
              egcd(m[i], m[j], CON);    //计算两个mi和mj的最大公因数，他们的最大公因数赋值给CON
              if (compare(CON, mirvar(1))) {
                  break;      //最大公因数不为1的话，不满足条件
              }
           }
       }
       if (compare(CON, mirvar(1)))
       {
           break;
       }
    }
    if (compare(CON, mirvar(1)))
    {
       printf("Input positive integer m（s）don not match the requirement of Chinese remainder theorem!");     //上述部分判断是否满足中国剩余定理的条件
    }
    else    //满足中国剩余定理条件
    {
       i = 0;
       for (i = 0; i < k; i++)    //将所有的m乘起来，这里的Mul相当于书上的m（所有的mj的乘积）
       {
           multiply(Mul, m[i], Mul);   //multiply函数作用：Mul=Mul*m[i]
       }
       i = 0;
       for (i = 0; i < k; i++)
       {
           fdiv(Mul, m[i], Mi);//fdiv函数作用：Mi = Mul/ m[i];
           xgcd(Mi, m[i], Mj, Mj, Mj);//Mj = invers(Mi, m[i]);   计算两个大数的扩展最大公约数，也可以用来计算模逆，z=gcd(x,y)=x.xd+y.yd，在本例中是求得Mi的模逆为Mj
           multiply(Mi, Mj, Mm);  //Mm=Mi*Mj
           multiply(Mm, a[i], x[i]);//x[i] = (Mi*Mj*a[i]);
       }
       i = 0;
       for (i = 0; i < k; i++)
       {
           add(Xul, x[i], Xul);//X = X + x[i];
       }
       powmod(Xul, mirvar(1), Mul, re);//re = X % M;
       printf("the resault of x is :n");
       cotnum(re, stdout);
    }
    mirkill(Xul);
    mirkill(Mul);
    mirkill(CON);
    mirkill(Mi);
    mirkill(Mj);
    mirkill(Mm);
    mirkill(re);
    return 0;
}