1. 程式人生 > >miracl庫c語言實現中國剩餘定理

miracl庫c語言實現中國剩餘定理

 

#include"miracl.h"
#include"mirdef.h"
#include<stdio.h>

#define NUM 3 //從文字讀入方程組個數,資料以分行儲存

FILE *fp;

int flag = 0; 

int main()
{
	int i, j;
	char tmp[100];
	big a[NUM], m[NUM], x[NUM], mm1[NUM], mm2[NUM];        /*
	big t0, t1, M, m1, X, y, W;                             *這裡mm1[] 為 Mj
	miracl *mip = mirsys(500, 10);                          *mm2[] 為 M^-1 j
	mip->IOBASE = 10;                                       *M 為 m[] 的連乘
	for(i = 0; i < NUM; i++)                               */
	{
		a[i] = mirvar(0);
		m[i] = mirvar(0);
		x[i] = mirvar(0);
		mm1[i] = mirvar(0);
		mm2[i] = mirvar(0);
	}
	t0 = mirvar(0);
	t1 = mirvar(1);
	M = mirvar(1);
	m1 = mirvar(0); 
	X = mirvar(0);
	y = mirvar(1);
	W = mirvar(0);

	fp = fopen("C:\\Users\\asus\\Desktop\\02.txt", "r");
	while( !feof(fp) )
	{
		fgets(tmp, 100, fp); // 按行讀入資料
		//printf("%s", tmp);
		if(flag<NUM)
		{
			cinstr(a[flag], tmp);
			cotnum(a[flag], stdout);
			flag++;
		}
		else
		{
			cinstr(m[flag-NUM], tmp);
			cotnum(m[flag-NUM], stdout);
			flag++;
		}
	}
	printf("\n");
	fclose(fp);
	
	for(i=0; i<NUM; i++) //判斷mi是否兩兩互素
	{
		for(j=i+1; j<NUM; j++)
		{
			egcd(m[i], m[j], t0);
			if(divisible(t0, t1)) continue;
			else 
			{
				printf("不能直接應用中國剩餘定理\n");
				exit(0);
			}
		}
	}

	for(i=0; i<NUM; i++)
	{
		multiply(m[i], M, M); //計算M=m1 x m2 x ... x mn
	}

	copy(M, W);
	for(i=0; i<2; i++)
	{
		divide(M, m[i], mm1[i]); //計算mm1 = Mj ???運算一次後M的值就清0了!!!釋放空間了???
		xgcd(mm1[i], m[i], mm2[i], m1, t1 ); //計算mm2 = M1~j
		copy(W, M);
	}

	for(i=0; i<NUM; i++) //計算xj
	{
		multiply(mm1[i], mm2[i], t0);
		multiply(t0, a[i], x[i]);
	}
	
	for(i=0; i<NUM; i++) //計算xj和
	{
		add(x[i], X, X);
	}

	
	powmod(X, y, M, X);//X mod M

	cotnum(X, stdout);
	mirexit();
	
	return 0;
	
}