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

一、实现目标：

根据中国剩余定理，设正整数两两互素，那么对于任意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的最大公因数，他们的最大公因数赋值给CONif (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的模逆为Mjmultiply(Mi, Mj, Mm);  //Mm=Mi*Mjmultiply(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;}

实验结果展示：