【信息安全数学基础】中国剩余定理

  1 extern "C"
  2 {
  3 #include<miracl.h>
  4 #include<mirdef.h>
  5 }
  6 #include<stdio.h>
  7 #include<stdlib.h>
  8 
  9 #define NUM 3 //从文本读入方程组个数，数据以分行存储
 10 
 11 int main()
 12 {
 13     //中国剩余定理：x≡M1M1^1a1+M2M2^-1a2+……+MkMk^-1ak(mod m)其中m=m1m2m……mk,Mj=m/mj,MjMj^-1≡1(mod mj),j=1,2,……k
 14     int i, j;//设定循环变量
 15     FILE *fp;
 16     int flag = 0;
 17     big a[NUM], m[NUM], x[NUM], Mj[NUM], Mj1[NUM]; //设定大数数组，a存放同余号右边，m存放模值，x存放同余号左边，Mj与字面同义，Mj1存放Mj^-1
 18     big t0, t1, M, m1, X, y, W;                             
 19     miracl *mip = mirsys(4000, 10); //开辟空间
 20     mip->IOBASE = 10;                                     
 21     for (i = 0; i < NUM; i++)//初始化数组
 22     {
 23         a[i] = mirvar(0);
 24         m[i] = mirvar(0);
 25         x[i] = mirvar(0);
 26         Mj[i] = mirvar(0);
 27         Mj1[i] = mirvar(0);
 28     }
 29     t0 = mirvar(0);//初始化各变量
 30     t1 = mirvar(1);
 31     M = mirvar(1);
 32     m1 = mirvar(0);
 33     X = mirvar(0);
 34     y = mirvar(1);
 35     W = mirvar(0);
 36 
 37     fp = fopen("2.txt", "r");
 38     char tmp[4000];
 39     while (1)
 40     {
 41         fscanf(fp, "%s", tmp);
 42         if (flag<NUM)//从文件中读入a列并存入数组
 43         {
 44             cinstr(a[flag], tmp);
 45             cotnum(a[flag], stdout);
 46             flag++;
 47         }
 48         else//从文件中读入m列存入数组
 49         {
 50             cinstr(m[flag - NUM], tmp);
 51             cotnum(m[flag - NUM], stdout);
 52             flag++;
 53         }
 54         if (flag == 2 * NUM) break;
 55     }
 56     printf("结果为：");
 57     fclose(fp);
 58 
 59     for (i = 0; i<NUM; i++) //判断m列是否两两互素
 60     {
 61         for (j = i + 1; j<NUM; j++)
 62         {
 63             egcd(m[i], m[j], t0);//t0=存放公约数，若gcd(m[i], m[j])=1即m两两互素进行下一次判断，否则跳出
 64             if (mr_compare(t0, t1)) continue;
 65             else
 66             {
 67                 printf("不能直接应用中国剩余定理\n");
 68                 exit(0);
 69             }
 70         }
 71     }
 72 
 73     for (i = 0; i<NUM; i++)
 74     {
 75         multiply(m[i], M, M); //计算M=m1m2……mk
 76     }
 77 
 78     copy(M, W);//W=M
 79     for (i = 0; i<2; i++)
 80     {
 81         divide(M, m[i], Mj[i]); //计算Mj
 82         xgcd(Mj[i], m[i], Mj1[i], m1, t1); //计算Mj^1(mod mj)
 83         copy(W, M);//M=W，数值还原
 84     }
 85 
 86     for (i = 0; i<NUM; i++) //计算xj及其累加结果X
 87     {
 88         multiply(Mj[i], Mj1[i], t0);
 89         multiply(t0, a[i], x[i]);
 90         add(x[i], X, X);
 91     }
 92 
 93 
 94     powmod(X, y, M, X);//X mod M
 95 
 96     cotnum(X, stdout);
 97     mirexit();
 98 
 99     return 0;
100 
101 }