i6r.c

calc大数库

MPMATR *TRANSPOSER(MPMATR *Mptr)
/* 
 * returns the transpose of *Mptr.
 */
{
	MPMATR *Nptr;
	unsigned int i, j;

	Nptr = BUILDMATR(Mptr->C, Mptr->R);
	for (j = 0; j <= Nptr->R - 1; j++)
		for (i = 0; i <= Nptr->C - 1; i++)
			elt(Nptr, j, i) = COPYR(elt(Mptr, i, j));
	return (Nptr);
}

MPMATR *MULTMATR(MPMATR *Mptr, MPMATR *Nptr)
/*
 * Here Mptr->C = Nptr->R. 
 * returns (*Mptr) * (*Nptr).
 */
{
	MPR *X, *Y, *Temp;
	unsigned int i, j, k;
	MPMATR *Lptr;

	Lptr = BUILDMATR(Mptr->R, Nptr->C);
	for (i = 0; i <= Mptr->R - 1; i++)
	{
		for (j = 0; j <= Nptr->C - 1; j++)
		{
			X = ZEROR();
			for (k = 0; k <= Mptr->C - 1; k++)
			{
				Y = MULTR(elt(Mptr, i, k), elt(Nptr, k, j));
				Temp = X;
				X = ADDR(X, Y);
				FREEMPR(Temp);
				FREEMPR(Y);
			}
			elt(Lptr, i, j) = X;
		}
	}
	return (Lptr);
}

void SHORTEST(MPMATI *AA, MPI *I[], USI filename, USI number)
/*
 * Input: An m x m matrix of integers AA = Q arising from the EXTGCD
 * LLLGCD or LLLGCD0. The last row of AA is assumed to be a multiplier
 * of shortest length.
 * Output:  All multipliers of shortest length.
 * filename:1->egcdmult.out,2->lllgcdmult.out,3->lllgcd0mult.out.
 * 4->gcd3.out (used in GCD3()).
 * If number = 0, no printing.
 * If number > 0 prints out multipliers on screen and in .out file.
 * Only used in nfunc.c.
 */
{
	unsigned int i, j, k, m, n, *l, *u, *t, *x, w, count = 0;
	MPR **X, **Y, **L, **T, **U, **V, *Temp0, *Temp1, *Temp2, *Temp3; 
        MPR *C, *Z, *ONE, *SUM, **N, *temp, *tempR, **XX;
	MPI *lengthi;
	char buff[25]; 
	FILE *outfile;
	MPMATR *A, *B, *BB, *Q, *QQ, *MATR, *VECTOR;
	MPMATI *MATI;
	int e;
	
/* First get the vector [V[0],...,V[n-2], where V[i]=<A[n-1][*],A[i][*]>.*/
	m = AA->C;
	n = (AA->R) - 1;
	A = BUILDMATR(n + 1, m);
	for (i = 0; i < n + 1; i++)
		for (j = 0; j < m; j++)
		{
			elt(A, i, j) = BUILDMPR();
			elt(A, i, j)->N = COPYI(AA->V[i][j]);
			elt(A, i, j)->D = ONEI();
		}

	V = (MPR **)mmalloc(n * sizeof(MPR *));
	for (i = 0; i < n; i++)
	{
		V[i] = ZEROR();
		for (j = 0; j < m; j++)
		{
			Temp1 = V[i];
			Temp2 = MULTR(elt(A, n, j), elt(A, i, j));
			V[i] = ADDR(V[i], Temp2);
			FREEMPR(Temp1);
			FREEMPR(Temp2);
		}
	}
	B = TRANSPOSER(A);
	BB = MULTMATR(A, B);
	FREEMATR(B);
	Q = CHOLESKYR(BB);
	FREEMATR(BB);
	QQ = COPYMATR(Q);
	for (i = 0; i < n; i++)
	{
		FREEMPR(elt(QQ, i, i));
		elt(QQ, i, i) = ONER();
	}
/* Calculate (U^{-1})^tV = Y */
	MATR = QQ;
	QQ = TRANSPOSER(QQ);
	FREEMATR(MATR);
	Y = (MPR **)mmalloc(n * sizeof(MPR *));
	Y[0] = COPYR(V[0]);
	for (w = 1; w < n; w++)
	{
		SUM = ZEROR();
		for ( k = 0; k < w; k++)
		{
			Temp1 = SUM;
			Temp2 = MULTR(elt(QQ, w, k), Y[k]);
			SUM = ADDR(SUM, Temp2);
			FREEMPR(Temp1);
			FREEMPR(Temp2);
		}
		Y[w] = SUBR(V[w], SUM);
		FREEMPR(SUM);
	}
	FREEMATR(QQ);
	N = (MPR **)mmalloc(n * sizeof(MPR *));
	for (j = 0; j < n; j++)
		N[j] = RATIOR(Y[j], elt(Q, j, j));
	if (filename == 1)
		strcpy(buff, "egcdmult.out");
	else if (filename == 2)
		strcpy(buff, "lllgcdmult.out");
	else if (filename == 3)
		strcpy(buff, "lllgcd0mult.out");
	else if (filename == 4)
		strcpy(buff, "gcd3.out");
	else if (filename == 5)
		strcpy(buff, "axb.out");
	outfile = fopen(buff, "a");
	ONE = ONER();
	X = (MPR **)mmalloc(n * sizeof(MPR *));
	L = (MPR **)mmalloc(n * sizeof(MPR *));
	T = (MPR **)mmalloc(n * sizeof(MPR *));
	U = (MPR **)mmalloc(n * sizeof(MPR *));
	l = (USI *)mmalloc(n * sizeof(USI));
	t = (USI *)mmalloc(n * sizeof(USI));
	u = (USI *)mmalloc(n * sizeof(USI));
	x = (USI *)mmalloc(n * sizeof(USI));
	for (j = 0; j < n; j++)
	{
		l[j] = 0;
		t[j] = 0;
		u[j] = 0;
		x[j] = 0;
	}
	C = ZEROR();
	for (i = 0; i < n; i++)
	{
		Temp1 = C;
		Temp2 = MULTR(Y[i], Y[i]);
		Temp3 = RATIOR(Temp2, elt(Q, i, i));
		C = ADDR(C, Temp3);
		FREEMPR(Temp1);
		FREEMPR(Temp2);
		FREEMPR(Temp3);
	}
	i = n - 1;
	T[i] = COPYR(C); t[i] = 1;
	FREEMPR(C);
	U[i] = ZEROR(); u[i] = 1;
bounds:
	Z = RATIOR(T[i], elt(Q, i, i));
	Temp2 = SUBR(N[i], U[i]);
	if (l[i])
		FREEMPR(L[i]); 
	else
		l[i] = 1;
	L[i] = INTROOT(Z, Temp2);
	FREEMPR(Temp2);
	Temp0 = SUBR(U[i], N[i]);
	Temp2 = INTROOT(Z, Temp0); 
	FREEMPR(Temp0);
	Temp3 = MINUSR(Temp2);
	if (x[i])
		FREEMPR(X[i]);
	else
		x[i] = 1;
	X[i] = SUBR(Temp3, ONE);
	FREEMPR(Z);
	FREEMPR(Temp2);
	FREEMPR(Temp3);
loop:
	Temp1 = X[i];
	X[i] = ADDR(X[i], ONE);
	FREEMPR(Temp1);
	if (COMPAREI(X[i]->N, L[i]->N) == 1)
	{
		i++;
		if (i == n)
			goto exiting;
		else
			goto loop;
	}
	else
	{
		if (i)
		{
			Temp1 = ADDR(X[i], U[i]);
			Temp0 = SUBR(Temp1, N[i]);
			FREEMPR(Temp1);
			Temp2 = MULTR(Temp0, Temp0);
			FREEMPR(Temp0);
			Temp1 = MULTR(elt(Q, i, i), Temp2);
			FREEMPR(Temp2);
			if (t[i - 1])
				FREEMPR(T[i - 1]);
			else
				t[i - 1] = 1;
			T[i - 1] = SUBR(T[i], Temp1);
			FREEMPR(Temp1);
			i--;
			SUM = ZEROR();
			for (j = i + 1; j < n; j++)
			{
				Temp1 = SUM;
				Temp2 = MULTR(elt(Q, i, j), X[j]);
				SUM = ADDR(SUM, Temp2);
				FREEMPR(Temp1);
				FREEMPR(Temp2);
			}
			if (u[i])
				FREEMPR(U[i]);
			else
				u[i] = 1;
			U[i] = COPYR(SUM);
			FREEMPR(SUM);
			goto bounds;
		}
		else
			goto found;
	}
	
	found:
		VECTOR = BUILDMATR(1, AA->R);
		XX = (MPR **)mmalloc((USL)(n * sizeof(MPR *)));
		for (j = 0; j < n; j++)
		{
			tempR = BUILDMPR();
			tempR->N = COPYI(I[j]);
			tempR->D = ONEI();
			XX[j] = ADDR(X[j], tempR);
			FREEMPR(tempR);
		}
if (number)
{
		printf("\tMULTIPLIER[%u] = ", count + 1);
		fprintf(outfile, "\tMULTIPLIER[%u]:", count + 1);
		printf("b[%u]", n + 1);
		fprintf(outfile, "b[%u]", n + 1);
		for (j = 0; j < n; j++)
		{
			e = (XX[j]->N)->S;
			if (e == -1)
			{
				printf("+");
				fprintf(outfile, "+");
				temp = MINUSR(XX[j]);
				if (!EQONER(temp))
				{
					PRINTR(temp);
					FPRINTR(outfile, temp);
				}
				printf("b[%u]", j + 1);
				fprintf(outfile, "b[%u]", j + 1);
				FREEMPR(temp);
			}
			if (e == 1)
			{
				printf("-");
				fprintf(outfile, "-");
				if (!EQONER(XX[j]))
				{
					PRINTR(XX[j]);
					FPRINTR(outfile, XX[j]);
				}
				printf("b[%u]", j + 1);
				fprintf(outfile, "b[%u]", j + 1);
			}
		}
}
		for (j = 0; j < n; j++)
		{
			elt(VECTOR, 0, j) = MINUSR(X[j]);
			FREEMPR(XX[j]);
		}
		ffree((char *)XX, n * sizeof(MPR *));
		elt(VECTOR, 0, n) = ONER();
		MATR = MULTMATR(VECTOR, A);
		FREEMATR(VECTOR);
if (number)
{
		printf("=");
		fprintf(outfile, "=");
}
		MATI = BUILDMATI(1, m);
		for (j = 0; j < m; j++)
			MATI->V[0][j] = COPYI((elt(MATR, 0, j))->N);
if (number)
{
		for (j = 0; j < m; j++)
		{
			PRINTI(MATI->V[0][j]);
			FPRINTI(outfile, MATI->V[0][j]);
			printf(" ");
			fprintf(outfile, " ");
		}
		fprintf(outfile, ": ");
		printf(": ");
}
		lengthi = LENGTHSQRI(MATI, 0);
if (number)
{

		FPRINTI(outfile, lengthi);
		fprintf(outfile, "\n");
		PRINTI(lengthi);
		printf("\n");
}
		FREEMPI(lengthi);
		FREEMATR(MATR);
		FREEMATI(MATI);
		count++;
		goto loop;
exiting:
		FREEMPR(ONE);
		for (j = 0; j < n; j++)
		{
			FREEMPR(N[j]);
			FREEMPR(X[j]);
			FREEMPR(L[j]);
			FREEMPR(T[j]);
			FREEMPR(U[j]);
			FREEMPR(Y[j]);
			FREEMPR(V[j]);
		}
		ffree((char *)N, n * sizeof(MPR *));
		ffree((char *)V, n * sizeof(MPR *));
		ffree((char *)Y, n * sizeof(MPR *));
		ffree((char *)X, n * sizeof(MPR *));
		ffree((char *)L, n * sizeof(MPR *));
		ffree((char *)T, n * sizeof(MPR *));
		ffree((char *)U, n * sizeof(MPR *));
		ffree((char *)l, n * sizeof(USI));
		ffree((char *)t, n * sizeof(USI));
		ffree((char *)u, n * sizeof(USI));
		ffree((char *)x, n * sizeof(USI));
		FREEMATR(Q);
		FREEMATR(A);
if (number)
{
		if (count > 1)
		{
			fprintf(outfile, "\tThere are exactly %u shortest multiplier vectors\n", count);
			printf("\tThere are exactly %u shortest multiplier vectors\n", count);
		}
		if (count == 1)
		{
			printf("\tThere is exactly one shortest multiplier vector\n");
			fprintf(outfile, "\tThere is exactly one shortest multiplier vector\n");
		}
}
		fclose(outfile);
		return;
}

void SHORTESTTT(MPMATI *AA, MPI *BOUND)
/*
 * This is used in INHOMFP().
 * Input: An m x M matrix of integers AA whose first m - 1 rows Q[0],... * ,Q[m-2] are LI and span a lattice L.
 * Output: If P is the last row, finds all x[0],...,x[m-2]
 * satisfying ||P -x[0]Q[0]-...-x[m-2]Q[m-2]||^2<=BOUND.
 */
{
	unsigned int i, j, k, m, n, *l, *u, *t, *x, w, count = 0;
	unsigned int M, FLAG;
	MPR **X, **Y, **L, **T, **U, **V, *Temp0, *Temp1, *Temp2, *Temp3; 
        MPR *C, *Z, *ONE, *SUM, **N, *temp;
	MPI *lengthi, *lengthj, *TempI1, *TempI2;
	MPMATR *A, *B, *BB, *Q, *QQ, *MATR, *VECTOR;
	MPMATI *MATI;
	int e;
	char buff[25]; 
	FILE *outfile;
	
	strcpy(buff, "inhomfp.out");
	outfile = fopen(buff, "w");
/* Compute ||AA[m - 1][*]||^2. */
	m = AA->R;
	M = AA->C;
	n = m - 1;
	lengthj = ZEROI();
	for (j = 0; j < M; j++)
	{
		TempI1 = lengthj;
		TempI2 = MULTI(AA->V[n][j], AA->V[n][j]);
		lengthj = ADDI(lengthj, TempI2);
		FREEMPI(TempI1);
		FREEMPI(TempI2);
	}
/* First get the vector [V[0],...,V[n-2], where V[i]=<A[n-1][*],A[i][*]>.*/
	A = BUILDMATR(m, M);
	for (i = 0; i < m; i++)
		for (j = 0; j < M; j++)
		{
			elt(A, i, j) = BUILDMPR();
			elt(A, i, j)->N = COPYI(AA->V[i][j]);
			elt(A, i, j)->D = ONEI();
		}

	V = (MPR **)mmalloc(n * sizeof(MPR *));
	for (i = 0; i < n; i++)
	{
		V[i] = ZEROR();
		for (j = 0; j < M; j++)
		{
			Temp1 = V[i];
			Temp2 = MULTR(elt(A, n, j), elt(A, i, j));
			V[i] = ADDR(V[i], Temp2);
			FREEMPR(Temp1);
			FREEMPR(Temp2);
		}
	}
	B = TRANSPOSER(A);
	BB = MULTMATR(A, B);
	FREEMATR(B);
	Q = CHOLESKYR(BB);
	FREEMATR(BB);
	QQ = COPYMATR(Q);
	for (i = 0; i < n; i++)
	{
		FREEMPR(elt(QQ, i, i));
		elt(QQ, i, i) = ONER();
	}
/* Calculate (U^{-1})^tV = Y */
	MATR = QQ;
	QQ = TRANSPOSER(QQ);
	FREEMATR(MATR);
	Y = (MPR **)mmalloc(n * sizeof(MPR *));
	Y[0] = COPYR(V[0]);
	for (w = 1; w < n; w++)
	{
		SUM = ZEROR();
		for ( k = 0; k < w; k++)
		{
			Temp1 = SUM;
			Temp2 = MULTR(elt(QQ, w, k), Y[k]);
			SUM = ADDR(SUM, Temp2);
			FREEMPR(Temp1);
			FREEMPR(Temp2);
		}
		Y[w] = SUBR(V[w], SUM);
		FREEMPR(SUM);
	}
	FREEMATR(QQ);
	N = (MPR **)mmalloc(n * sizeof(MPR *));
	for (j = 0; j < n; j++)
		N[j] = RATIOR(Y[j], elt(Q, j, j));
	ONE = ONER();
	X = (MPR **)mmalloc(n * sizeof(MPR *));
	L = (MPR **)mmalloc(n * sizeof(MPR *));
	T = (MPR **)mmalloc(n * sizeof(MPR *));
	U = (MPR **)mmalloc(n * sizeof(MPR *));
	l = (USI *)mmalloc(n * sizeof(USI));
	t = (USI *)mmalloc(n * sizeof(USI));
	u = (USI *)mmalloc(n * sizeof(USI));
	x = (USI *)mmalloc(n * sizeof(USI));
	for (j = 0; j < n; j++)
	{
		l[j] = 0;
		t[j] = 0;
		u[j] = 0;
		x[j] = 0;
	}
	C = ZEROR();
	for (i = 0; i < n; i++)
	{
		Temp1 = C;
		Temp2 = MULTR(Y[i], Y[i]);
		Temp3 = RATIOR(Temp2, elt(Q, i, i));
		C = ADDR(C, Temp3);
		FREEMPR(Temp1);
		FREEMPR(Temp2);
		FREEMPR(Temp3);
	}
	
	Temp0 = BUILDMPR();
	Temp0->N = SUBI(BOUND, lengthj);
	Temp0->D = ONEI();
	Temp1 = C;
	C = ADDR(C, Temp0);
	FREEMPR(Temp0);
	FREEMPR(Temp1);
		
	i = n - 1;
	T[i] = COPYR(C); t[i] = 1;
	FREEMPR(C);
	U[i] = ZEROR(); u[i] = 1;
bounds:
	Z = RATIOR(T[i], elt(Q, i, i));
	Temp2 = SUBR(N[i], U[i]);
	if (l[i])
		FREEMPR(L[i]); 
	else
		l[i] = 1;
	L[i] = INTROOT(Z, Temp2);
	/*printf("at L[%u]\n", i);*/
	FREEMPR(Temp2);
	Temp0 = SUBR(U[i], N[i]);
	Temp2 = INTROOT(Z, Temp0); 
	FREEMPR(Temp0);
	Temp3 = MINUSR(Temp2);
	if (x[i])
		FREEMPR(X[i]);
	else
		x[i] = 1;
	X[i] = SUBR(Temp3, ONE);
	FREEMPR(Z);
	FREEMPR(Temp2);
	FREEMPR(Temp3);
loop:
	Temp1 = X[i];
	X[i] = ADDR(X[i], ONE);
	FREEMPR(Temp1);
	if (COMPAREI(X[i]->N, L[i]->N) == 1)
	{
		i++;
		if (i == n)
			goto exiting;
		else
			goto loop;
	}
	else
	{
		if (i)
		{
			Temp1 = ADDR(X[i], U[i]);
			Temp0 = SUBR(Temp1, N[i]);
			FREEMPR(Temp1);
			Temp2 = MULTR(Temp0, Temp0);
			FREEMPR(Temp0);
			Temp1 = MULTR(elt(Q, i, i), Temp2);
			FREEMPR(Temp2);
			if (t[i - 1])
				FREEMPR(T[i - 1]);
			else
				t[i - 1] = 1;
			T[i - 1] = SUBR(T[i], Temp1);
			FREEMPR(Temp1);
			i--;
			SUM = ZEROR();
			for (j = i + 1; j < n; j++)
			{
				Temp1 = SUM;
				Temp2 = MULTR(elt(Q, i, j), X[j]);
				SUM = ADDR(SUM, Temp2);
				FREEMPR(Temp1);
				FREEMPR(Temp2);
			}
			if (u[i])
				FREEMPR(U[i]);
			else
				u[i] = 1;
			U[i] = COPYR(SUM);
			FREEMPR(SUM);
			goto bounds;
		}
		else
			goto found;
	}
	
	found:
		VECTOR = BUILDMATR(1, AA->R);
		for (j = 0; j < n; j++)
		{
			/*
			printf("X[%u][%u] = ", count, j);
			PRINTR(X[j]);
			printf("\n");
			fprintf(outfile, "X[%u][%u] = ", count, j);
			FPRINTR(outfile, X[j]);
			fprintf(outfile,"\n");
*/
			elt(VECTOR, 0, j) = COPYR(X[j]);
		}
		elt(VECTOR, 0, n) = ZEROR();
		MATR = MULTMATR(VECTOR, A);
		FREEMATR(VECTOR);
		printf("LV[%u] = ", count + 1);
		fprintf(outfile, "LV[%u] = ", count + 1);