nfunc.c

calc大数库

				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, 2, 1);
	}
	for (j = 0; j < p; j++)
		FREEMPI(XX[j]);
	ffree((char *)XX, p * sizeof(MPI *));
	FREEMATI(MATI3);
	FREEMATI(Q);
	GCDVERBOSE = 0;
	HAVASFLAG = 0;
	return;
}

void JACOBIGCDX()
/*
 * This executes the LLLGCD algorithm of Havas, Majewski and Matthews.
 * 30/1/97.
 */
{
	USI answer, m, i;
	MPMATI *MATI1, *MATI2, *MATI3;
	MPI *A, *T;

	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);
	m = MATI1->R;
	MLLLVERBOSE = 0;
	printf("VERBOSE? (Y/N)\n");
	answer = GetYN();
	printf("answer = %d\n", answer);
	if (answer)
		MLLLVERBOSE = 1;
	printf("MLLLVERBOSE = %u\n", MLLLVERBOSE);
	MATI2 = JACOBIGCD(MATI1, &A, m);
	MLLLVERBOSE = 0;
	FREEMATI(MATI1);
	HAVASFLAG = 0;
	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 JACOBI are\n");
	for (i = 0;i < MATI2->R;i++)
	{
		PRINTI(MATI2->V[m - 1][i]); 
		printf(" ");
	}
	printf("\n");
	MATI3 = BUILDMATI(1, m);
	for (i = 0;i < m;i++)
		MATI3->V[0][i] = COPYI(MATI2->V[m - 1][i]); 
	T = LENGTHSQRI(MATI3, 0);
	printf("The Euclidean norm squared = ");
	PRINTI(T);
	FREEMPI(T);
	FREEMATI(MATI2);
	FREEMATI(MATI3);
	printf("\n");
	
}

void SCHNORRHERMITE(MPI *N)
{
	USI answer, m1, n1, m, n, t, i, j;
	MPMATI *MATI1, *MATI2, *MATI0;
	MPI *Temp;
	char buff[20];
	FILE *outfile;

	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;
	HERMITEVERBOSE = 0;		
	printf("VERBOSE? (Y/N)\n");
	answer = GetYN();
	printf("answer = %d\n", answer);
	if (answer)
	{
		MLLLVERBOSE = 1;
		HERMITEVERBOSE = 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();
	n = MATI1->R;
	m = MATI1->C;
	t = n + m;
	MATI0 = BUILDMATI(n, t);
	for (i = 0; i < n; i++)
	{
		for (j = 0; j < n; j++)
		{
			if ( i == j)
				MATI0->V[i][j] = ONEI();
			else
				MATI0->V[i][j] = ZEROI();
		}
		for (j = n; j < t; j++)
		{
			Temp = POWERI(N, t - j);
			MATI0->V[i][j] = MULTI(MATI1->V[i][j - n], Temp);
			FREEMPI(Temp);
		}
	}
	FREEMATI(MATI1);
	strcpy(buff, "shermitemat.out");
	outfile = fopen(buff, "w");
	FPRINTMATI(outfile,0,MATI0->R-1,0,MATI0->C-1,MATI0);
	fclose(outfile);
	MATI2 = BASIS_REDUCTION000(MATI0, m1, n1, N);
	HERMITEVERBOSE = 0;		
	FREEMATI(MATI0);
	for (i = 0; i < n; i++)
	{
		for (j = n; j < t; j++)
		{
			Temp = POWERI(N, t - j);
			MATI2->V[i][j] = INTI(MATI2->V[i][j], Temp);
			FREEMPI(Temp);
		}
	}
	printf("The corresponding reduced basis is\n");
	PRINTMATI(0,MATI2->R-1,0,MATI2->C-1,MATI2);
	strcpy(buff, "shermitebas.out");
	outfile = fopen(buff, "w");
	FPRINTMATI(outfile,0,MATI2->R-1,0,MATI2->C-1,MATI2);
	fclose(outfile);
	FREEMATI(MATI2);
	return;
}

void LLLHERMITE1X()
{
	USI answer, rank, m1, n1;
	MPMATI *MATI1, *MATI2, *MATI3;
	char buff[20];
	FILE *outfile;

	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);
	HERMITE1VERBOSE = 0;
	printf("VERBOSE? (Y/N)\n");
	answer = GetYN();
	if (answer)
		HERMITE1VERBOSE = 1;
	printf("enter the parameters m1 and n1 (normally 1 and 1) :");
	scanf("%u %u", &m1, &n1);
	Flush();
	MATI2 = LLLHERMITE1(MATI1, &MATI3, &rank, m1, n1);
	HERMITE1VERBOSE = 0;
	printf("rank = %u\n\n", rank);
	printf("The transformation matrix is\n");
	PRINTMATI(0,MATI2->R-1,0,MATI2->C-1,MATI2);
	strcpy(buff, "lllhermitetrans.out");
	outfile = fopen(buff, "w");
	FPRINTMATI(outfile,0,MATI2->R-1,0,MATI2->C-1,MATI2);
	fclose(outfile);
	printf("The row echelon Form is\n");
	PRINTMATI(0,MATI3->R-1,0,MATI3->C-1,MATI3);
	strcpy(buff, "lllhermitebas.out");
	outfile = fopen(buff, "w");
	FPRINTMATI(outfile,0,MATI3->R-1,0,MATI3->C-1,MATI3);
	fclose(outfile);
	FREEMATI(MATI1);
	FREEMATI(MATI2);
	FREEMATI(MATI3);
	return;
}

void CLEMENS(MPI *N)
/* constructs Clemens' matrix and then applies LLLHERMITE with alpha=3/4. */
{
	MPI *I;
	USI rank;
	USL i, j, n;
	MPMATI *MATI1, *MATI2, *MATI3;

	n = CONVERTI(N);
	MATI1 = BUILDMATI(n, n);
	for (j = 0; j < n; j++)
		MATI1->V[0][j] = ONEI();
	for (i = 1; i < n; i++)
		MATI1->V[i][0] = ZEROI();
	for (i = 1; i < n; i++)
	{
		I = CHANGEI(i);
		for (j = 1; j < n; j++)
			MATI1->V[i][j] = MPOWER_(j, I, N);
		FREEMPI(I);
	}
	PRINTMATI(0, MATI1->R -1, 0, MATI1->C -1, MATI1);
	GetReturn();
	MATI2 = LLLHERMITE1(MATI1, &MATI3, &rank, 3, 4);
	printf("The Hermite normal form = \n");
	PRINTMATI(0, MATI2->R -1, 0, MATI2->C -1, MATI2);
	printf("The unimodular transformation matrix = \n");
	PRINTMATI(0, MATI3->R -1, 0, MATI3->C -1, MATI3);
	FREEMATI(MATI1);
	FREEMATI(MATI2);
	FREEMATI(MATI3);
}

void CLEMENS1(MPI *N)
/* constructs Clemens' matrix and then applies Kannan-Bachem. */
{
	MPI *I;
	USI *alpha, nz;
	USL i, j, n;
	MPMATI *MATI1, *MATI2;

	n = CONVERTI(N);
	MATI1 = BUILDMATI(n, n);
	for (j = 0; j < n; j++)
		MATI1->V[0][j] = ONEI();
	for (i = 1; i < n; i++)
		MATI1->V[i][0] = ZEROI();
	for (i = 1; i < n; i++)
	{
		I = CHANGEI(i);
		for (j = 1; j < n; j++)
			MATI1->V[i][j] = MPOWER_(j, I, N);
		FREEMPI(I);
	}
	alpha = KB_ROW(MATI1, &nz);
	printf("rank = %u\n", nz);
	MATI2 = HERMITE1(MATI1, alpha, nz);
	printf("The Hermite normal form = \n");
	PRINTMATI(0, MATI2->R -1, 0, MATI2->C -1, MATI2);
	FREEMATI(MATI2);
	ffree((char *)alpha, (MATI1->C) * sizeof(USI));
	FREEMATI(MATI1);
}

void CLEMENS2(MPI *N)
/* constructs Clemens' matrix and then applies Kannan-Bachem. */
/* Also produces the transformation matrix. */
{
	MPI *I;
	USI *alpha, nz;
	USL i, j, n;
	MPMATI *MATI1, *MATI2, *MATI3, *MATI4;

	n = CONVERTI(N);
	MATI1 = BUILDMATI(n, n);
	for (j = 0; j < n; j++)
		MATI1->V[0][j] = ONEI();
	for (i = 1; i < n; i++)
		MATI1->V[i][0] = ZEROI();
	for (i = 1; i < n; i++)
	{
		I = CHANGEI(i);
		for (j = 1; j < n; j++)
			MATI1->V[i][j] = MPOWER_(j, I, N);
		FREEMPI(I);
	}
	alpha = KB_ROWP(MATI1, &MATI3, &nz);
	printf("rank = %u\n", nz);
	MATI2 = HERMITE1P(MATI1, MATI3, &MATI4, alpha, nz);
	FREEMATI(MATI3);
	printf("The Hermite normal form = \n");
	PRINTMATI(0, MATI2->R -1, 0, MATI2->C -1, MATI2);
	FREEMATI(MATI2);
	printf("The unimodular transformation matrix = \n");
	PRINTMATI(0, MATI4->R -1, 0, MATI4->C -1, MATI4);
	FREEMATI(MATI4);
	ffree((char *)alpha, (MATI1->C) * sizeof(USI));
	FREEMATI(MATI1);
}

void AXB()
{
	USI answer, m1, n1, rankA, nullA, i, j, k, r, s, flag, p, m, n;
	MPMATI *MATI1, *MATI2, *MATI3, *MATI4, *MATI5, *MATI6, *MATI7;
	MPMATI *Tmp, *Q, *M;
	char buff[20];
	FILE *outfile;
	MPI **XX, **X, *Temp;
	int e;

	printf("Do you wish to use an existing augmented matrix from a file? (Y/N)\n");
	answer = GetYN();
	if (answer)
	{
		Tmp = FINPUTMATFILEI_I();
		if (Tmp == NULL)
			exit (0);
	}
	else
	{
		printf("enter the augmented matrix A|B of integers (first column non-zero) :\n");
		Tmp = INPUTMATI();
	}
	MATI1 = TRANSPOSEI(Tmp);
	FREEMATI(Tmp);
	printf("The matrix entered in transposed form 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;
	MATI3 = IDENTITYI(MATI1->R);
	printf("enter the parameters m1 and n1 (normally 1 and 1) :");
	scanf("%u %u", &m1, &n1);
	Flush();
	MATI2 = BASIS_REDUCTION(MATI1, &MATI3, 0, m1, n1);
	flag = 0;
	r =  MATI3->R - 1;
	s =  MATI3->C - 1;
	k = 0;
	while (k < r){
		if (!EQZEROI(MATI3->V[k][s])){
			flag = 1;
			break;
		}
		k++;
	}
	if (flag){
		printf("No solution of AX=B\n");
		FREEMATI(MATI1);
		FREEMATI(MATI2);
		FREEMATI(MATI3);
		MLLLVERBOSE = 0;
		return;
	}
	printf("\n\n");
	if (MLLLVERBOSE){
		printf("The interim transformation matrix is\n");
		PRINTMATI(0,MATI3->R-1,0,MATI3->C-1,MATI3);
		GetReturn();
	}
/*
	printf("The corresponding reduced basis for Rowspace(A^t) is\n");
	PRINTMATI(0,MATI2->R-1,0,MATI2->C-1,MATI2);
	GetReturn();
	strcpy(buff, "axbimbas.out");
	outfile = fopen(buff, "w");
	FPRINTMATI(outfile,0,MATI2->R-1,0,MATI2->C-1,MATI2);
	fclose(outfile);
*/
	rankA = MATI2->R;
	nullA = (MATI1->R) - (rankA + 1);
	if (!nullA){
		printf("Unique solution\n");
		MATI6 = BUILDMATI(1, MATI1->R - 1);
		for (j = 0; j < MATI1->R - 1; j++)
			MATI6->V[0][j] = MINUSI(MATI3->V[MATI3->R - 1][j]);
		FREEMATI(MATI1);
		FREEMATI(MATI2);
		FREEMATI(MATI3);
		printf("The solution X of XA^t=B^t is\n");
		PRINTMATI(0,MATI6->R-1,0,MATI6->C-1,MATI6);
		strcpy(buff, "axbsol.out");
		outfile = fopen(buff, "w");
		FPRINTMATI(outfile,0,MATI6->R-1,0,MATI6->C-1,MATI6);
		fclose(outfile);
		FREEMATI(MATI6);
		MLLLVERBOSE = 0;
		return;
	}
	MATI4 = BUILDMATI(1 + nullA, MATI1->R - 1);
	for (i = 0; i <= nullA; i++){
		for (j = 0; j < MATI1->R - 1; j++){
			MATI4->V[i][j] = COPYI(MATI3->V[rankA + i][j]);
		}
        }
        GCDFLAG = 1;
	MATI5 = BASIS_REDUCTION0(MATI4, m1, n1);
	MLLLVERBOSE = 0;
	MATI6 = BUILDMATI(1, MATI1->R - 1);
        for (j = 0; j < MATI1->R - 1; j++)
	{
		(MATI5->V[MATI5->R - 1][j])->S = -((MATI5->V[MATI5->R - 1][j])->S);
                MATI6->V[0][j] = COPYI(MATI5->V[MATI5->R - 1][j]);
	}
        GCDFLAG = 0;
	MATI7 = DELETE_ROWI(MATI5->R, MATI5);
	printf("The short basis for nullspace(A^t) is\n");
	PRINTMATI(0,MATI7->R-1,0,MATI7->C-1,MATI7);
	GetReturn();
	strcpy(buff, "axbbas.out");
	outfile = fopen(buff, "w");
	FPRINTMATI(outfile,0,MATI7->R-1,0,MATI7->C-1,MATI7);
	fclose(outfile);
	p = nullA;
	m = p + 1;

	printf("A short solution X of XA^t=B^t is \n");
	printf("b[%u] = ",m); PRINTMATI(0,MATI6->R-1,0,MATI6->C-1,MATI6);
/*
	strcpy(buff, "axbsol.out");
	outfile = fopen(buff, "w");
	FPRINTMATI(outfile,0,MATI6->R-1,0,MATI6->C-1,MATI6);
	fclose(outfile);
*/
	printf("Do you want to get the shortest multipliers using Fincke_Pohst? (Y/N)\n");
	answer = GetYN();
	XX = (MPI **)mmalloc((USL)(p * sizeof(MPI *)));
	for (j = 0; j < p; j++)
		XX[j] = ZEROI();
	if (answer)
	{
		GCDVERBOSE = 1;
		Q = MATI5;
		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(MATI6);
				MATI6 = M;
			}
		}
		printf("found a shortest solution vector:\n");
	}
	else
		Q = MATI5;
	strcpy(buff, "axb.out");
	outfile = fopen(buff, "w");
	if (answer)
		fprintf(outfile, "A shortest solution vector is ");
	else
		fprintf(outfile, "A short multiplier is ");
	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 (!EQONE