lll.c

calc大数库

					printf("gcdflag set\n");
				K = norig + 1;
			}
			else
				K = k;
			for (l = K - 2; l >= 1; l--)
			{
				Flag = 	STEP40(k, l, &L, &B1ptr, D);
				if (Flag)
				{
					FREEMPI(Z);
					sigma++;
					if (MLLLVERBOSE)
						printf("relation vector number %u found\n", sigma);
					tau = k++;
					goto found; /* STEP 9 of POHST. */
				}
			}
			k++;
		}
		if (!gcdflag)
			FREEMPI(Z);
	}
	FREEMPI(M1);
	FREEMPI(N1);
	FREEMATI(C);
	FREEMATI(L);
	for (i = 0; i <= norigg; i++)
		FREEMPI(D[i]);
	ffree((char *)D, (1 + norigg) * sizeof(MPI *));
	/* Now to improve the last n - norig rows of *B1ptr, using Gauss
lattice reduction - pointed out by George Havas in Sims' book - 8/11/94: */


/*
	strcpy(buff, "progress");
	outfile = fopen(buff, "a+");
	FPRINTMATI(outfile,0,B1ptr->R-1,0,B1ptr->C-1,B1ptr);
*/


	REPEAT = 1;
	temp1 = BUILDMATI(1, B1ptr->C);
	for (k = 0; k < B1ptr->C; k++)
		temp1->V[0][k] = ZEROI();
	iterate = 0;
	QSORTMATI(B1ptr, 0, norig - 1);
	while (REPEAT)
	{
		REPEAT = 0;
		for (j = 1; j < norig; j++)
		{
			t4 = DOTRI(B1ptr, j, j);
			for (i = 0; i < j; i++)
			{
				t1 = DOTRI(B1ptr, j, i);
				t2 = DOTRI(B1ptr, i, i);
				t3 = NEAREST_INTI(t1, t2);
				FREEMPI(t1);
				FREEMPI(t2);
				if (t3->S == 0)
				{
					FREEMPI(t3);
					continue;
				}
				t1 = t3;
				t3 = MINUSI(t3);
				FREEMPI(t1);
				for (k = 0; k < B1ptr->C; k++)
				{
					FREEMPI(temp1->V[0][k]);
					X = MULTI(t3, B1ptr->V[i][k]);
					temp1->V[0][k] = ADDI(X, B1ptr->V[j][k]);
					FREEMPI(X);
				}
				t1 = DOTRI(temp1, 0, 0);
				tt = RSV(t1, t4);
				tt = RSV(t1, t4);
				if (tt == -1)
				{
					ADD_MULT_ROWI0(i, j, t3, B1ptr);
					if (MLLLVERBOSE)
					{
						printf("Gauss-improving nullspace row %u: length squared was = ", j + 1);
						PRINTI(t4);
						printf("\n");
						printf("lengthsquared now = ");
						PRINTI(t1);
						printf("\n");
					}
/*
					fprintf(outfile, "Gauss-improving nullspace row %u: length squared was = ", j + 1);
					FPRINTI(outfile, t4);
					fprintf(outfile, "\n");
					fprintf(outfile, "lengthsquared now = ");
					FPRINTI(outfile, t1);
					fprintf(outfile, "\n");
					fflush(outfile);
*/
					FREEMPI(t4);
					t4 = t1;
					REPEAT = 1;
				}
				else
					FREEMPI(t1);
				FREEMPI(t3);
			}
				FREEMPI(t4);
		}
		iterate++;
		QSORTMATI(B1ptr, 0, norig - 1);
	}
	REPEAT = 1;
	iterate = 0;
	while (REPEAT)
	{
		REPEAT = 0;
		for (j = norig; j < B1ptr->R; j++)
		{
			t4 = DOTRI(B1ptr, j, j);
			for (i = 0; i < norig; i++)
			{
				t1 = DOTRI(B1ptr, j, i);
				t2 = DOTRI(B1ptr, i, i);
				t3 = NEAREST_INTI(t1, t2);
				FREEMPI(t1);
				FREEMPI(t2);
				if (t3->S == 0)
				{
					FREEMPI(t3);
					continue;
				}
				t1 = t3;
				t3 = MINUSI(t3);
				FREEMPI(t1);
				for (k = 0; k < B1ptr->C; k++)
				{
					FREEMPI(temp1->V[0][k]);
					X = MULTI(t3, B1ptr->V[i][k]);
					temp1->V[0][k] = ADDI(X, B1ptr->V[j][k]);
					FREEMPI(X);
				}
				t1 = DOTRI(temp1, 0, 0);
				tt = RSV(t1, t4);
				tt = RSV(t1, t4);
				if (tt == -1)
				{
					ADD_MULT_ROWI0(i, j, t3, B1ptr);
					if (MLLLVERBOSE)
					{
						printf("Gauss-improving row %u: length squared was = ", j + 1);
						PRINTI(t4);
						printf("\n");
						printf("lengthsquared now = ");
						PRINTI(t1);
						printf("\n");
					}
/*
					fprintf(outfile, "Gauss-improving row %u: length squared was = ", j + 1);
					FPRINTI(outfile, t4);
					fprintf(outfile, "\n");
					fprintf(outfile, "lengthsquared now = ");
					FPRINTI(outfile, t1);
					fprintf(outfile, "\n");
					fflush(outfile);
*/
					FREEMPI(t4);
					t4 = t1;
					REPEAT = 1;
				}
				else
					FREEMPI(t1);
				FREEMPI(t3);
			}
			FREEMPI(t4);
		}
		iterate++;
	}
	FREEMATI(temp1);
/*
	fclose(outfile);
*/
	return (B1ptr);
}

MPMATI *BASIS_REDUCTION000(MPMATI *Bptr, USI m1, USI n1, MPI *N)
/*
 * Input: *Bptr, a matrix of MPI's, whose first row is not zero.
 * Output: a pointer to an MPMATI whose rows form a reduced basis for 
 * the lattice spanned by the rows of *Bptr. This basis is reduced in the 
 * sense of the paper "Factoring polynomials with rational coefficients" by
 * A. K. Lenstra, H. W. Lenstra and L. Lovasz, Math. Ann. 261, 515-534 (1982)
 * using the modified version in "Solving exponential Diophantine equations
 * using lattice basis reduction algorithms" by B. M. M. De Weger, J. No. Theory
 * 26, 325-367 (1987). No change of basis matrix is returned.
 * De Weger's algorithm has been changed to cater for arbitrary matrices. The
 * the rows are now in general linearly dependent. 
 * We use the fact that the Gram Schmidt process detects the first row 
 * which is a linear combination of the preceding rows. We employ a modification
 * of the LLL algorithm outlined by M. Pohst in J. Symbolic Computation (1987)4,
 * 123-127.  We call this the MLLL algorithm.
 * If we are using this algorithm to find small multipliers for the extended 
 * gcd problem, GCDFLAG is set in EXTGCD() and gcdflag is set below.
 * m1 / n1 is usually taken to be 3 / 4.
 */
{
	unsigned int i, k, l, n, m, t, flag = 0, Flag = 0, gcdflag = 0;
	unsigned int flagg, beta, K1 = 0, tau = 2, sigma = 0, rho;
	MPI **D, *X, *Y, *Z, *H, *Tmp, *R, *M1, *N1;
	MPMATI *C, *L, *B1ptr;
	unsigned int norig;

	Z = NULL;
	m = Bptr->C;
	n = Bptr->R;
	norig = n;
	B1ptr = COPYMATI(Bptr);
	D = (MPI **)mmalloc((1 + n) * sizeof(MPI *));
	D[0] = ONEI();
	for (i = 1; i <= n; i++)
		D[i] = ZEROI();
	C = ZEROMNI(n, m);
	L = ZEROMNI(n, n);

	found:
	n = B1ptr->R;
	i = (K1 == 0) ? 1 : K1;
	/* K1 = no. of consecutive rows of *B1ptr that don't need updating 
for the Gram Schmidt process. */
	while (i <= B1ptr->R)
	{
		BASIS_UPDATE(i, m, &C, &L, B1ptr, D);
		flag = 1;
		for (t = 0; t < m; t++)
		{
			if (!EQZEROI(C->V[i - 1][t]))
			{
				flag = 0;	
				break;
			}
		}	
		if (flag)
			break;
		X = ZEROI();
		for (t = 0; t < m; t++)
		{
			H = MULTI(C->V[i - 1][t], C->V[i - 1][t]);
			Tmp = X;
			X = ADDI(X, H);
			FREEMPI(Tmp);
			FREEMPI(H);
		}
		FREEMPI(D[i]);
		D[i] = INT(X, D[i - 1]);	
		FREEMPI(X);
		i++;
		if (MLLLVERBOSE)
			printf("i = %u\n", i);
	}
	beta =  (flag) ? i : i - 1;
	rho = K1 = i - 1;
	if (MLLLVERBOSE)
		printf("completed updating the basis\n");
/* Here K1 = no. of LI rows in *B1ptr found by Gram Schmidt process.
   flag = 0 means all the rho = beta rows of *B1ptr are LI;
   flag = 1 means that the first rho = beta - 1 rows of *B1ptr are LI, but the
   beta-th row is a LC of the preceding rows. So beta = number of rows of *B1ptr
   currently being examined by the LLL algorithm. */
	k = tau;
	M1 = CHANGE(m1);
	N1 = CHANGE(n1);
	while (k <= beta)
	{
		if (MLLLVERBOSE)
			printf("beta - k = %u\n", beta -k);
		l = k - 1;
		Flag = STEP4000(k, l, &L, &B1ptr, D, N);
		if (k >= norig && GCDFLAG)
		{
			gcdflag = 1;
			goto FOUND;
		}
		if (Flag)/* STEP 9 of POHST. */
		{
			sigma++;
			if (MLLLVERBOSE)
				printf("relation vector number %u found\n", sigma);
			tau = k++;
			goto found;
		}
		X = MULTI(D[k - 2], D[k]);
		Y = MULTI(D[k - 1], D[k - 1]);
		Tmp = Y;
		Y = MULT_I(Y, m1);
		FREEMPI(Tmp);
		R = MULTI(L->V[k - 1][k - 2], L->V[k - 1][k - 2]);
		Z = ADD0I(X, R);
		Tmp = Z;
		Z = MULT_I(Z, n1);
		FREEMPI(Tmp);
		if (RSV(Y, Z) == 1)/*& STEP 5 of POHST. */
		{
			flagg = 0;
			if (EQZEROI(D[k]) && EQZEROI(R))
			{/* CASE B=0 of STEP 7 of POHST. */
				FREEMPI(D[k - 1]);
				D[k - 1] = ZEROI();
				STEP8000(k, &B1ptr, &L, N);
				if (k - 1 < K1)
					K1 = k - 1;
				/* The swap may have changed 2nd last row */
				/* of *B1ptr. */
				for (t = 0; t < m; t++)
				{
					FREEMPI(C->V[k - 2][t]);
					C->V[k - 2][t] = ZEROI();
				}
				beta--;
				flagg = 1;
				if (k > 2)
					k--;
				FREEMPI(X);
				FREEMPI(Y);
				FREEMPI(R);
				FREEMPI(Z);
				continue;
			}
			if (flagg == 0)
			{
				for (i = k + 1; i <= beta; i++)
					STEP7(i, k, &L, D);
			}
			STEP8000(k, &B1ptr, &L, N);
			if (k - 2 < K1)
				K1 = k - 2;
			/* swap will change last two rows of *B1ptr. */
			FREEMPI(R);
			FREEMPI(Y);
			Y = MULTI(L->V[k - 1][k - 2], L->V[k - 1][k - 2]);
			Tmp = Y;
			Y = ADD0I(Y, X);
			FREEMPI(X);
			FREEMPI(Tmp);
			Tmp = D[k - 1];
			D[k - 1] = INT0(Y, D[k - 1]);
			FREEMPI(Tmp);
			FREEMPI(Y);
			if (k > 2)
				k--;
		}
		else
		{ /* STEP 6 of POHST. */
			FREEMPI(R);
			FREEMPI(X);
			FREEMPI(Y);
FOUND:
			for (l = k - 2; l >= 1; l--)
			{
				Flag = 	STEP4000(k, l, &L, &B1ptr, D, N);
				if (Flag)
				{
					FREEMPI(Z);
					sigma++;
					if (MLLLVERBOSE)
						printf("relation vector number %u found\n", sigma);
					tau = k++;
					goto found; /* STEP 9 of POHST. */
				}
			}
			k++;
		}
		if (!gcdflag)
			FREEMPI(Z);
		else
			break;
	}
	FREEMPI(M1);
	FREEMPI(N1);
	FREEMATI(C);
	FREEMATI(L);
	for (i = 0; i <= norig; i++)
		FREEMPI(D[i]);
	ffree((char *)D, (1 + norig) * sizeof(MPI *));
	if (MLLLVERBOSE)
	{
		printf("number of basis vectors found = %u ;\n", rho);
		printf("number of relation vectors found = %u .\n", sigma);
	}
	return (B1ptr);
}

unsigned int STEP4000(USI k, USI l, MPMATI **Lptr, MPMATI **Bptr, MPI *D[], MPI *N)
/*
 * updates *Lptr, *Bptr.
 * returns 1 if row k of *Bptr becomes zero, returns zero otherwise.
 */
{
	unsigned int i,  j, flag = 1, t, m, n;
	MPI *X, *Y, *R, *Tmp, *Temp;
	MPMATI *TmpMATI;

	m = (*Bptr)->C;
	Y = MULT_I((*Lptr)->V[k - 1][l - 1], 2);
	if (RSV(Y, D[l]) == 1)
	{
		R = NEAREST_INTI((*Lptr)->V[k - 1][l - 1], D[l]);
		X = MINUSI(R);
		TmpMATI = *Bptr;
		*Bptr = ADD_MULT_ROWI(l - 1, k - 1, X, *Bptr);
		FREEMATI(TmpMATI);
	n = (*Bptr)->R;
	for (i = 0; i < n; i++)
	{
		for (j = n; j < m; j++)
		{
			Temp = POWERI(N, m - j);
			(*Bptr)->V[i][j] = INTI((*Bptr)->V[i][j], Temp);
			FREEMPI(Temp);
		}
	}
		if (HERMITEVERBOSE)
		{
		printf("Row %u -> Row %u + ", k,k);PRINTI(X);printf(" x Row %u\n", l);
		PRINTMATI(0,(*Bptr)->R-1,0,(*Bptr)->C-1,*Bptr);
		GetReturn();
		}
	for (i = 0; i < n; i++)
	{
		for (j = n; j < m; j++)
		{
			Temp = POWERI(N, m - j);
			(*Bptr)->V[i][j] = MULTI((*Bptr)->V[i][j], Temp);
			FREEMPI(Temp);
		}
	}
/*
		MAXI = UPDATEMAXI(MAXI, *Bptr);
*/
		FREEMPI(X);
		for (j = 1; j < l; j++)
		{
			X = MULTI((*Lptr)->V[l - 1][j - 1], R);
			Tmp = (*Lptr)->V[k - 1][j - 1];
			(*Lptr)->V[k - 1][j - 1] = SUBI((*Lptr)->V[k - 1][j - 1], X);
			FREEMPI(Tmp);
			FREEMPI(X);
		}
		X = MULTI(D[l], R);
		Tmp = (*Lptr)->V[k - 1][l - 1];
		(*Lptr)->V[k - 1][l - 1] = SUBI((*Lptr)->V[k - 1][l - 1], X);

		FREEMPI(Tmp);
		FREEMPI(X);
		FREEMPI(R);
	}
	for (t = 0; t < m; t++)
	{
		if (!EQZEROI((*Bptr)->V[k - 1][t]))
		{
			flag = 0;	
			break;
		}
	}	
	if (flag)
	{
		TmpMATI = *Bptr;
		*Bptr = DELETE_ROWI(k, *Bptr);
		FREEMATI(TmpMATI);
	}
	FREEMPI(Y);
	return (flag);
}
void STEP8000(USI k, MPMATI **B1ptr, MPMATI **Lptr, MPI *N)
{
	MPI *T, *Temp;
	USI i, j, n, m;

	*B1ptr = SWAP_ROWSI1(k - 2, k - 1, *B1ptr);
	T = COPYI((*Lptr)->V[k - 1][k - 2]);
	*Lptr = SWAP_ROWSI1(k - 2, k - 1, *Lptr);
	FREEMPI((*Lptr)->V[k - 1][k - 2]);
	(*Lptr)->V[k - 1][k - 2] = T;
	FREEMPI((*Lptr)->V[k - 2][k - 2]);
	(*Lptr)->V[k - 2][k - 2] = ZEROI();

	n = (*B1ptr)->R;
	m = (*B1ptr)->C;
	for (i = 0; i < n; i++)
	{
		for (j = n; j < m; j++)
		{
			Temp = POWERI(N, m - j);
			(*B1ptr)->V[i][j] = INTI((*B1ptr)->V[i][j], Temp);
			FREEMPI(Temp);
		}
	}
	if (HERMITEVERBOSE)
	{
		printf("Swapping Rows %u and %u\n", k-1,k);
		PRINTMATI(0,(*B1ptr)->R-1,0,(*B1ptr)->C-1,*B1ptr);
		GetReturn();
	}
	for (i = 0; i < n; i++)
	{
		for (j = n; j < m; j++)
		{
			Temp = POWERI(N, m - j);
			(*B1ptr)->V[i][j] = MULTI((*B1ptr)->V[i][j], Temp);
			FREEMPI(Temp);
		}
	}

	return;
}