i6r.c

calc大数库

		FLAG = 0;
		for (j = 0; j < n; j++)
		{
			e = (X[j]->N)->S;
			if (e == -1)
			{
				if (FLAG == 0)
					FLAG = 1;
				printf("-");
				fprintf(outfile, "-");
				temp = MINUSR(X[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)
			{
				if (FLAG)
				{
					printf("+");
					fprintf(outfile, "+");
				}
				if (FLAG == 0)
					FLAG = 1;
				if (!EQONER(X[j]))
				{
					PRINTR(X[j]);
					FPRINTR(outfile, X[j]);
				}
				printf("b[%u]", j + 1);
				fprintf(outfile, "b[%u]", j + 1);
			}
		}
		printf("=");
		fprintf(outfile, "=");
		MATI= BUILDMATI(1, A->C);
		for (j = 0; j < A->C; j++)
		{
			MATI->V[0][j] = COPYI((elt(MATR, 0, j))->N);
			PRINTI(MATI->V[0][j]);
			FPRINTI(outfile, MATI->V[0][j]);
			printf(" ");
			fprintf(outfile, " ");
			FREEMPI(MATI->V[0][j]);
			MATI->V[0][j] = SUBI(AA->V[n][j], ((elt(MATR, 0, j))->N));
		}
		printf(": ||b[%u]-LV[%u]||^2=", AA->R, count + 1);
		fprintf(outfile, ": ||b[%u]-LV[%u]||^2=", AA->R,  count + 1);
		lengthi = LENGTHSQRI(MATI, 0);
		PRINTI(lengthi);
		FPRINTI(outfile, lengthi);
		printf("\n");
		fprintf(outfile, "\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));
	fclose(outfile);
	FREEMATR(Q);
	FREEMATR(A);
	FREEMPI(lengthj);
	if (count > 1)
		printf("There are exactly %u lattice vectors\n", count);
	if (count == 1)
		printf("There is exactly one lattice vector\n");
	return;
}

MPMATI *SHORTESTT0(MPMATI *AA, MPI **XX[])
/*
 * Used in EXTGCD(), LLLGCD() and LLLGCD0() in nfunc.c.
 * 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.
 * Problem: If P is the last row, find a shorter multiplier
 * P-x[0]Q[0]-...-x[m-2]Q[m-2]
 */
{
	unsigned int i, j, k, m, n, *l, *u, *t, *x, w;
	unsigned int FOUND = 0, M;
	MPR **X, **Y, **L, **T, **U, **V, *Temp0, *Temp1, *Temp2, *Temp3; 
        MPR *C, *Z, *ONE, *SUM, **N;
	MPI *lengthi, *lengthj, *TempI1, *TempI2;
	MPMATR *A, *B, *BB, *Q, *QQ, *MATR, *VECTOR;
	MPMATI *MATI;
	int tt;

	MATI = NULL;
/* Compute ||AA[m - 1][*]||^2. */
	m = AA->R;
	M = AA->C;
	n = m - 1;
	*XX = (MPI **)mmalloc((USL)(n * sizeof(MPI *)));
	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);
	}
	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, m);
		for (j = 0; j < n; j++)
			elt(VECTOR, 0, j) = MINUSR(X[j]);

		elt(VECTOR, 0, m - 1) = ONER();
		MATR = MULTMATR(VECTOR, A);
		FREEMATR(VECTOR);
		MATI = BUILDMATI(1, M);
		for (j = 0; j < M; j++)
			MATI->V[0][j] = COPYI((elt(MATR, 0, j))->N);
		lengthi = LENGTHSQRI(MATI, 0);
		tt = RSV(lengthj, lengthi);
		FREEMPI(lengthi);
		FREEMATR(MATR);
		if (tt == 1)
		{
			FREEMPI(lengthj);
			FOUND = 1;
			for (j = 0; j < n; j++)
				(*XX)[j] = COPYI(X[j]->N);
			goto exiting;
		}
		else
			FREEMATI(MATI);
		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 (FOUND)
	{
		if (GCDVERBOSE)
		{
			printf("found a shorter multiplier:\n");
			PRINTMATI(0, MATI->R - 1, 0, MATI->C - 1, MATI);
			lengthi = LENGTHSQRI(MATI, 0);
			printf("LENGTHSQUARED = ");
			PRINTI(lengthi);
			FREEMPI(lengthi);
			printf("\n");
			printf("search continuing:\n\n");
		}
		return(MATI);
	}
	else
	{
		FREEMPI(lengthj);
		if (GCDVERBOSE)
			printf("search ended.\n");
		for (j = 0; j < n; j++)
			(*XX)[j] = ZEROI();
		return(NULL);
	}
}

MPI ***SHORTESTX(MPMATI *AA, MPI *I[], USI *counter)
/*
 * Input: An m x m matrix of integers AA = Q arising from LLLGCD. 
 * The last row of AA is assumed to be a multiplier of shortest length.
 * Output: the XX[] defining the multipliers of shortest length.
 * Only used in GCD_CONJ in nfunc.c.
 * Warning: Array XX below only caters for up to GCD_MAX=100000 shortest multipliers.
 * If necessary, increase it and recompile.
 */
{
	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;
	MPI ***XX, *tempI;
	MPMATR *A, *B, *BB, *Q, *QQ, *MATR;
	
	XX = (MPI ***)mmalloc((USL)(GCD_MAX * sizeof(MPI **)));
/* 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));
	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:
		XX[count] = (MPI **)mmalloc((USL)(n * sizeof(MPI *)));
		for (j = 0; j < n; j++)
		{
			tempI = ADDI(X[j]->N, I[j]);
			XX[count][j] = MINUSI(tempI);
			FREEMPI(tempI);
		}
		count++;
		if (count > GCD_MAX)
		{
			fprintf(stderr, "count > %u\n", GCD_MAX);
			exit (1);
		}
		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 (count > 1)
			printf("There are exactly %u shortest multiplier vectors\n", count);
		if (count == 1)
			printf("There is exactly one shortest multiplier vector\n");
*/
		*counter = count;
		return (XX);
}