nfunc.c

calc大数库

	printf(" The matrix entered is A = \n\n");
	PRINTMATI(0,MATI1->R-1,0,MATI1->C-1,MATI1);
	z = MIN((int) MATI1->R, (int) MATI1->C);
	printf("enter MAX_ENTRY:");
	MAX_ENTRY = INPUTI(&u);
	Flush();
	printf("The test matrix is %u x %u\n\n", MATI1->R, MATI1->C);
	printf("Enter alpha=m1/n1: select m1 and n1 (1/4 < alpha <= 1) :");
	scanf("%u %u", &m1, &n1);
	Flush();
	printf("MAX_ENTRY CUTOFF= "); PRINTI(MAX_ENTRY);printf("\n\n");
	MLLLVERBOSE = 0;
	printf("VERBOSE? (Y/N)\n\n");
	answer = GetYN();
	if (answer)
		MLLLVERBOSE = 1;
	printf("Do you want P and Q? (Y/N)\n\n");
	answer = GetYN();
	if (answer)
	{
		M = SMITHI(MATI1, &MATI2, &MATI3, &r, MAX_ENTRY, m1, n1);
		FREEMPI(MAX_ENTRY);
		printf("The row unit matrix is P = \n");
		PRINTMATI(0,MATI2->R-1,0,MATI2->C-1,MATI2);
		strcpy(buff, "smithp.out");
		outfile = fopen(buff, "w");
		FPRINTMATI(outfile,0,MATI2->R-1,0,MATI2->C-1,MATI2);
		fclose(outfile);
		printf("The column unit matrix is Q = \n");
		PRINTMATI(0,MATI3->R-1,0,MATI3->C-1,MATI3);
		strcpy(buff, "smithq.out");
		outfile = fopen(buff, "w");
		FPRINTMATI(outfile,0,MATI3->R-1,0,MATI3->C-1,MATI3);
		fclose(outfile);
		printf("Do you want to check that PAQ is actually equal to the quoted SNF?(Y/N)\n\n");
		answer = GetYN();
		if (answer){
			Temp = MULTMATI(MATI2,MATI1);
			Temp1 = MULTMATI(Temp,MATI3);
			for (i = 0; i < Temp1->R; i++)
			{
				for (j = 0; j < Temp1->C; j++)
				{
					if (i !=j && (Temp1->V[i][j])->S )
					{
						fprintf(stderr, "PAQ is not diagonal!\n");
						exit (1);
					}
				}
			}
			for(i = 0; i < r; i++)
			{
				if (!EQUALI(Temp1->V[i][i], M[i]))
				{
					fprintf(stderr, "PAQ is diagonal, but the diagonals are not right!\n");
					exit (1);
				}
				
			}
			FREEMATI(Temp);
			FREEMATI(Temp1);
		}
		FREEMATI(MATI2);
		FREEMATI(MATI3);
		printf("validated!\n\n");
	}
	else
	{
		M = SMITHI1(MATI1, &r, MAX_ENTRY, m1, n1);
		FREEMPI(MAX_ENTRY);
	}
	MLLLVERBOSE = 0;
	FREEMATI(MATI1);
	strcpy(buff, "smith.out");
	outfile = fopen(buff, "w");
	fprintf(outfile, "%u %u\n", r, 1);
	for(i = 0; i < r; i++)
	{
		printf("invariant factor d[%u] = ", i + 1);
		PRINTI(M[i]);
		FPRINTI(outfile, M[i]);
		fprintf(outfile, "\n");
		printf("\n");
		FREEMPI(M[i]);
	}
	printf("rank A = %u\n", r);
	fclose(outfile);
	ffree((char *)M, z * sizeof(MPI *));
	return;
}

void DECODEX(MPI *Eptr, MPI *Pptr, MPI *Qptr)
/*
 * *Eptr is the encrytion key, *Pptr and *Qptr are the RSA primes.
 * The deciphering key *Dptr is computed and file "encoded" is decoded.
 * The decoded message is sent to the screen and also to the file "decoded".
 */
{
	MPI *Rptr, *Dptr;
	MPI *Tmp1, *Tmp2, *Tmp3;

	Rptr = MULTI(Pptr, Qptr);
	Tmp1 = SUB0_I(Pptr, 1);
	Tmp2 = SUB0_I(Qptr, 1);
	Tmp3 = MULTI(Tmp1, Tmp2);
	printf("(p-1)(q-1) =  :" );PRINTI(Tmp3);printf("\n");
	Dptr = INVERSEM(Eptr, Tmp3);
	printf("d =  :" );PRINTI(Dptr);printf("\n");
	FREEMPI(Tmp1);
	FREEMPI(Tmp2);
	FREEMPI(Tmp3);
	DECODE(Dptr, Rptr);
	FREEMPI(Dptr);
	FREEMPI(Rptr);
	return;
}
/*
MPI *RSAEX()
{
	unsigned int u;
	MPI *Pptr, *Qptr, *Eptr;

	printf("enter a prime p > 355142:" );
	Pptr = INPUTI(&u);
	printf("enter another prime q > 355142:" );
	Qptr = INPUTI(&u);
	Flush(); 
	Eptr = RSAE(Pptr, Qptr);
	FREEMPI(Pptr);
	FREEMPI(Qptr);
	return (Eptr);
}
*/

MPI *GCD_ARRAYV(MPIA M, MPIA *Y)
/*
 * n is length of array
 * Returns d=gcd(M[0],...,M[n-1]) and an array Y[] of MPI's such that
 * d = M[0]*Y[0]+...+M[n-1]*Y[n-1]. Here n > 1.
 */
{
	MPI *D;
	int i;
	unsigned long n = M->size;
	MPMATI *Temp, *Temp1, *Q;
	
	Temp = BUILDMATI(n, 1);
	for (i = 0; i < n; i++)
	  Temp->V[i][0] = COPYI(M->A[i]);
	Temp1 = EXTGCD(Temp, &D, &Q, 3, 4);
	printf("The multipliers are\n");
	PRINTMATI(0,Temp1->R-1,0,Temp1->C-1,Temp1);
	*Y = BUILDMPIA();
	for (i = 0; i < n; i++)
	  ADD_TO_MPIA(*Y, Temp1->V[0][i], i);
	FREEMATI(Temp);
	FREEMATI(Temp1);
	FREEMATI(Q);
	return (D);

/*
	B = (MPI **)mmalloc((USL)(n * sizeof(MPI *)));
	for (i = 0; i < n; i++)
		B[i] = COPYI(M[i]);	
	for (i = 1; i < n; i++)
	{
		tmp = B[i];
		B[i] = GCD(B[i], B[i - 1]);
		FREEMPI(tmp);
	}
	D = COPYI(B[n - 1]);
	tmp = B[n - 1];
	B[n - 1] = COPYI(M[n - 1]);
	FREEMPI(tmp);
	k = ONEI();
	for (i = n - 1; i >= 1; i--)
	{
		G = EUCLIDI(B[i], B[i - 1], &U, &V);
		FREEMPI(G);
		tmp = B[i];
		B[i] = MULTI(U, k);
		FREEMPI(U);
		FREEMPI(tmp);
		if ( i == 1)
			break;
		tmp = k;
		k = MULTI(k, V);	
		FREEMPI(V);
		FREEMPI(tmp);
		tmp = B[i - 1];
		B[i - 1] = COPYI(M[i - 1]);
		FREEMPI(tmp);
	}
	tmp = B[0];
	B[0] = MULTI(k, V);
	FREEMPI(tmp);
	FREEMPI(k);
	FREEMPI(V);
	*Y = B;
	return (D);
*/
}

void EXTGCDX()
/* This is algorthm 1 of Havas, Majewski, Matthews. */
{
	USI answer, m1, n1, m, n, i, j, p;
	MPMATI *MATI1, *MATI2, *Q, *M, *MATI3;
	char buff[20];
	FILE *outfile;
	MPI *A, *T, **XX, **X, *Temp;
	int e;

	HAVASFLAG = 1;
	printf("Do you wish to enter your sequence of numbers from an existing column matrix? (Y/N)\n");
	answer = GetYN();
	if (answer)
	{
		MATI1 = FINPUTMATFILEI_I();
		if (MATI1 == NULL)
			exit (0);
	}
	else
	{
		printf("WARNING: Make sure the first integer in the sequence is non-zero) :\n");
		MATI1 = INPUTMATI();
	}
	printf("The matrix entered is\n\n");
	PRINTMATI(0,MATI1->R-1,0,MATI1->C-1,MATI1);
	MLLLVERBOSE = 0;
	printf("VERBOSE? (Y/N)\n");
	answer = GetYN();
	printf("answer = %d\n", answer);
	if (answer)
		MLLLVERBOSE = 1;
	printf("MLLLVERBOSE = %u\n", MLLLVERBOSE);
	printf("to enter alpha=m1/n1: select m1 and n1 (normally 1 and 1) :");
	scanf("%u %u", &m1, &n1);
	Flush();
	MATI3 = EXTGCD(MATI1, &A, &MATI2, m1, n1);
	FREEMATI(MATI1);
	printf("The multiplier matrix found by LLL is\n");
	PRINTMATI(0,MATI2->R-1,0,MATI2->C-1,MATI2);
	printf("gcd = "); PRINTI(A); printf("\n");
	FREEMPI(A);
	printf("The multipliers found by LLL are\n");
	PRINTMATI(0,MATI3->R-1,0,MATI3->C-1,MATI3);
	m = MATI2->C;
	T = LENGTHSQRI(MATI3, 0);
	printf("The Euclidean norm squared = ");
	PRINTI(T);
	FREEMPI(T);
	printf("\n");
	printf("Do you want to get a shortest multiplier using Fincke_Pohst? (Y/N)\n");
	answer = GetYN();
	p = m - 1;
	XX = (MPI **)mmalloc((USL)(p * sizeof(MPI *)));
	for (j = 0; j < p; j++)
		XX[j] = ZEROI();
	if (answer)
	{
		GCDVERBOSE = 1;
		Q = MATI2;
		n = Q->R;
		while (1)
		{
			M = SHORTESTT0(Q, &X);
			for (j = 0; j < p; j++)
			{
				Temp = XX[j];
				XX[j] = ADDI(XX[j], X[j]);
				FREEMPI(Temp);
				FREEMPI(X[j]);
			}
			ffree((char *)X, p * sizeof(MPI *));
			if (M == NULL)
				break;
			else
			{
				for (j = 0; j < Q->C; j++)
				{
					FREEMPI(Q->V[n - 1][j]);
					Q->V[n - 1][j] = COPYI(M->V[0][j]);
				}
				FREEMATI(MATI3);
				MATI3 = M;
			}
		}
		printf("found a shortest multiplier vector:\n");
	}
	else
		Q = MATI2;
	strcpy(buff, "egcdmat.out");
	outfile = fopen(buff, "w");
	FPRINTMATI(outfile,0,Q->R-1,0,Q->C-1,Q);
	fclose(outfile);
	strcpy(buff, "egcdmult.out");
	outfile = fopen(buff, "w");
	if (answer)
		fprintf(outfile, "A Shortest multiplier is ");
	else
	{
		fprintf(outfile, "A not necessarily shortest multiplier is ");
		fprintf(outfile, "b[%u]=", m);
	}
	if (answer)
	{
		printf("b[%u]", m);
		fprintf(outfile, "b[%u]", m);
		for (j = 0; j < p; j++)
		{
			e = XX[j]->S;
			if (e == -1)
			{
				printf("+");
				fprintf(outfile, "+");
				Temp = MINUSI(XX[j]);
				if (!EQONEI(Temp))
				{
					PRINTI(Temp);
					FPRINTI(outfile, Temp);
				}
				printf("b[%u]", j + 1);
				fprintf(outfile, "b[%u]", j + 1);
				FREEMPI(Temp);
			}
			if (e == 1)
			{
				printf("-");
				fprintf(outfile, "-");
				if (!EQONEI(XX[j]))
				{
					PRINTI(XX[j]);
					FPRINTI(outfile, XX[j]);
				}
				printf("b[%u]", j + 1);
				fprintf(outfile, "b[%u]", j + 1);
			}
		}
	}
	if (answer){
		printf("=");
		fprintf(outfile, "=");
		for (i = 0; i < MATI3->C; i++)
		{
			PRINTI(MATI3->V[0][i]); 
			printf(" ");
			FPRINTI(outfile, MATI3->V[0][i]); 
			fprintf(outfile," ");
			
		}
		T = LENGTHSQRI(MATI3, 0);
		printf(": ");
		fprintf(outfile, ": ");
		PRINTI(T);
		FPRINTI(outfile, T);
		printf("\n");
		fprintf(outfile,"\n");
		FREEMPI(T);
	}
	else
	{
		for (i = 0; i < MATI3->C; i++)
		{
			FPRINTI(outfile, MATI3->V[0][i]); 
			fprintf(outfile," ");
			
		}
		T = LENGTHSQRI(MATI3, 0);
		fprintf(outfile, ": ");
		FPRINTI(outfile, T);
		fprintf(outfile,"\n");
		FREEMPI(T);
	}
	fclose(outfile);
	if (answer)
	{
		printf("Do you want to get all the shortest multipliers? (Y/N)\n");
		answer = GetYN();
		if (answer)
			SHORTEST(Q, XX, 1, 1);
	}
	for (j = 0; j < p; j++)
		FREEMPI(XX[j]);
	ffree((char *)XX, p * sizeof(MPI *));
	FREEMATI(MATI3);
	FREEMATI(Q);
	HAVASFLAG = 0;
	GCDVERBOSE = 0;
	return;
}

void FINCKE_POHSTX()
{
	USI answer, i, j, m, n, u;
	MPMATI *MATI1;
	MPMATR *M;
	MPR *C;

	printf("Do you wish to use an existing matrix from a file? (Y/N)\n");
	answer = GetYN();
	if (answer)
	{
		MATI1 = FINPUTMATFILEI_I();
		if (MATI1 == NULL)
			exit (0);
	}
	else
	{
		printf("enter the matrix of integers (first row non-zero) :\n");
		MATI1 = INPUTMATI();
	}
	printf("The matrix entered is\n\n");
	PRINTMATI(0,MATI1->R-1,0,MATI1->C-1,MATI1);
	m = MATI1->R;
	n = MATI1->C;
	M = BUILDMATR(m, n);
	for (i = 0; i < m; i++)
	{
		for (j = 0; j < n; j++)
		{
			elt(M, i, j) = BUILDMPR();
			elt(M, i, j)->N = COPYI(MATI1->V[i][j]);
			elt(M, i, j)->D = ONEI();
		}
	}
	FREEMATI(MATI1);
	printf("Enter the upper bound for Q(X):");
	C = INPUTR(&u);
	FINCKE_POHST(M, C, 1);
	FREEMATR(M);
	FREEMPR(C);
	return;
}

void IMPROVEPX()
{
	USI answer, m1, n1, norig, m, i;
	MPMATI *MATI1, *MATI2;
	char buff[20];
	FILE *outfile, *infile;
	MPI ***temp;
	int k;

	strcpy(buff, "hermitep.out");
	infile = fopen(buff, "r");
	MATI1 = FINPUTMATI(infile);
	fscanf(infile, "%u", &m); /* m = no of rows to be improved */
	fclose(infile);
	norig = MATI1->R - m;/* norig = no of zero rows in HNF */
	if (norig == 0)
	{
		printf("HNF has no non-zero rows - P cannot be improved\n");
		FREEMATI(MATI1);
		return;	
	}
	temp = (MPI ***)mmalloc(m * sizeof(MPI **));
	for (i = 0; i < m; i++)
		temp[i] = MATI1->V[i];
	for (i = 0; i < norig; i++)
		MATI1->V[i] = MATI1->V[i + m];
	for (i = norig; i < MATI1->R; i++)
		MATI1->V[i] = temp[i - norig];
	printf("The matrix entered is\n\n");
	PRINTMATI(0,MATI1->R-1,0,MATI1->C-1,MATI1);
	MLLLVERBOSE = 0;
	printf("VERBOSE? (Y/N)\n");
	answer = GetYN();
	if (answer)
		MLLLVERBOSE = 1;
	printf("enter the parameters m1 and n1 (normally 1 and 1) :");
	scanf("%u %u", &m1, &n1);
	Flush();
	GCDFLAG = 1;
	MATI2 = BASIS_REDUCTION00(MATI1, m1, n1, norig);
	GCDFLAG = 0;
	printf("\n\n");
	strcpy(buff, "improvep.out");
	outfile = fopen(buff, "w");
	for (i = norig; i < MATI2->R; i++)
		temp[i - norig] = MATI2->V[i];
	for (k = norig - 1; k >= 0; k--)
		MATI2->V[k + m] = MATI2->V[k];
	for (i = 0; i < m; i++)
		MATI2->V[i] = temp[i];
	FPRINTMATI(outfile,0,MATI2->R-1,0,MATI2->C-1,MATI2);
	fclose(outfile);
	printf("The improved transformation matrix is\n");
	PRINTMATI(0,MATI2->R-1,0,MATI2->C-1,MATI2);
	FREEMATI(MATI1);
	FREEMATI(MATI2);
	ffree((char *)temp, m * sizeof(MPI **));
	return;
}

void QSORTMPIX()
/*
 * Input: an array A[0],...,A[q] of MPI's.
 * sorts nonnegative MPI's  A[m],...,A[n] in order of size.
 */
{
	USI m, n, p, i, u, q;
	MPI **A;

	printf("enter 0<=m<= n<=q (q+1= no of elements, starting from A[0])  :");
	scanf("%u %u %u", &m, &n, &q);
	p = q + 1;
	A = (MPI **)mmalloc(p * sizeof(MPI *));
	printf("enter the array A[0],...,A[%u] of %u MPIs:",q, p);
	for (i = 0; i < p; i++)
		A[i] = INPUTI(&u);
	for (i = 0; i < p; i++)
	{
		printf("A[%u] = ", i);
		PRINTI(A[i]);
		printf("\n");
	}
	QSORTMPI0(A, m, n);
	for (i = 0; i < p; i++)
	{
		printf("A[%u] = ", i);
		PRINTI(A[i]);
		FREEMPI(A[i]);
		printf("\n");
	}
	ffree((char *)A, p * sizeof(MPI *));
	return;
}

void QSORTMATIX()
{
	USI answer, r, s;
	MPMATI *MATI1;

	printf("Do you wish to use an existing matrix from a file? (Y/N)\n");
	answer = GetYN();
	if (answer)
	{
		MATI1 = FINPUTMATFILEI_I();
		if (MATI1 == NULL)
			exit (0);
	}
	else
	{
		printf("enter the matrix of integers (first row non-zero) :\n");
		MATI1 = INPUTMATI();
	}