lllgcd.c

calc大数库

	{
	for (i4 = p4; i4 <= q4; i4++)
	{
	for (i5 = p5; i5 <= q5; i5++)
	{
		x = GCDm((USL)i1, (USL)i2);
		y = GCDm(x, (USL)i3);
		z = GCDm(y, (USL)i4);
		u = GCDm(z, (USL)i5);
		if (u == 1)
		{
		printf("(i1,i2,i3,i4,i5)=(%u,%u,%u,%u,%u)\n", i1,i2,i3,i4,i5);
		MATI1 = BUILDMATI(m, 1);
		MATI1->V[0][0] = CHANGE((USL) i1);
		MATI1->V[1][0] = CHANGE((USL) i2);
		MATI1->V[2][0] = CHANGE((USL) i3);
		MATI1->V[3][0] = CHANGE((USL) i4);
		MATI1->V[4][0] = CHANGE((USL) i5);
		BB = LLLGCD0M(MATI1, &GCD, m, m1, n1);
		FREEMATI(MATI1);
		FREEMPI(GCD);
		MATI3 = BUILDMATI(1, m);
		for (l = 0; l < m; l++)
			MATI3->V[0][l] = COPYI(BB->V[m - 1][l]); 
		T1 = LENGTHSQRI(MATI3, 0);
		Q = BB;
		XX = (MPI **)mmalloc((USL)(p * sizeof(MPI *)));
		for (t = 0; t < p; t++)
			XX[t] = ZEROI();
		while (1)
		{
			M = SHORTESTT0(Q, &X);
			for (t = 0; t < p; t++)
			{
				Temp = XX[t];
				XX[t] = ADDI(XX[t], X[t]);
				FREEMPI(Temp);
				FREEMPI(X[t]);
			}
			ffree((char *)X, p * sizeof(MPI *));
			if (M == NULL)
				break;
			else
			{
				for (l = 0; l < Q->C; l++)
				{
					FREEMPI(Q->V[m - 1][l]);
					Q->V[m - 1][l] = COPYI(M->V[0][l]);
				}
				FREEMATI(MATI3);
				MATI3 = M;
			}
		}
		T2 = LENGTHSQRI(MATI3, 0);
		outfile = fopen(buff, "a");
		if (RSV(T2, T1) == -1)
		{
			fflush(outfile);
			fprintf(outfile, "(%u,%u,%u,%u,%u): ", i1,i2,i3,i4,i5);
			fprintf(outfile, "b[%u]", m);
			for (t = 0; t < p; t++)
			{
				e = XX[t]->S;
				if (e == -1)
				{
					fprintf(outfile, "+");
					Temp = MINUSI(XX[t]);
					if (!EQONEI(Temp))
						FPRINTI(outfile, Temp);
					fprintf(outfile, "b[%u]", t + 1);
					FREEMPI(Temp);
				}
				if (e == 1)
				{
					fprintf(outfile, "-");
					if (!EQONEI(XX[t]))
						FPRINTI(outfile, XX[t]);
					fprintf(outfile, "b[%u]", t + 1);
				}
			}
			fprintf(outfile, ": ");
			DIFF = SUB0I(T1, T2);
			FPRINTI(outfile, DIFF); 
			FREEMPI(DIFF);
			fprintf(outfile, "\n");
		}
		fclose(outfile);
		for (t = 0; t < p; t++)
			FREEMPI(XX[t]);
		ffree((char *)XX, p * sizeof(MPI *));
		FREEMATI(BB);
		FREEMATI(MATI3);
		FREEMPI(T1);
		FREEMPI(T2);
		}
	}
	}
	}
	}
	}
	return;
}

void LLLGCDMINUS(MPMATI **Pptr, MPMATI **Lptr, USI j)
{
	USI r, s, m;
	MPI *Temp;
	MPMATI *TempMATI;

	m = (*Pptr)->R;
	for (r = 0; r < m; r++)
		for (s = 0; s < r; s++)
		{
			if (r == j || s == j)
			{
				Temp = (*Lptr)->V[r][s];
				(*Lptr)->V[r][s] = MINUSI(Temp);
				FREEMPI(Temp);
			}
		}
	
	TempMATI = *Pptr;
	Temp = MINUS_ONEI();
	*Pptr = SCALAR_ROWI(j, Temp, *Pptr);
	FREEMPI(Temp);
	FREEMATI(TempMATI);
	if (GCDVERBOSE)
	{
		printf("Row %u -> - Row %u\n", j,j);
		printf("P = \n");
		PRINTMATI(0,(*Pptr)->R-1,0,(*Pptr)->C-1,*Pptr);
		GetReturn();
	}
	return;
}

void GCD6()
/* We compute an optimal multipliers for a[1],a[2],a[3],a[4],a[5],a[6]
 * in the range p1<=a[1]<=q1 etc. and rel prime.
 */
{
	USI l, m, p1, q1, p2, q2, p3, q3, p4, q4, p5, q5, p6, q6, m1, n1;
	USI i1, i2, i3, i4, i5, i6, p, t;
	int e;
	USL x, y, z, u, v, w;
	MPI *GCD, *T1, *T2, **XX, **X, *Temp;
	MPMATI *MATI1, *BB, *MATI3, *Q, *M;
	char buff[20];
	FILE *outfile;

	/*for (i2 = (i1+1 > p2 ? i1+1: p2); i2 <= q2; i2++)*/
	/*for (i3 = (i2+1 > p3 ? i2+1: p3); i3 <= q3; i3++)*/
	/*for (i4 = (i3+1 > p4 ? i3+1: p4); i4 <= q4; i4++)*/
	/*for (i5 = (i4+1 > p5 ? i4+1: p5); i5 <= q5; i5++)*/
	/*for (i6 = (i5+1 > p6 ? i5+1: p6); i6 <= q6; i6++)*/

	printf("Enter alpha=m/n,  1/4 < m/n <= 1:  (ie enter  m  and n)");
	scanf("%u%u", &m1, &n1);
	strcpy(buff, "gcd6.out");
	if(access(buff, R_OK) == 0)
	  unlink(buff);
	printf("Enter  the ranges pi,qi; 2 <= pi < qi < 2^32, i = 1,...,6: ");
	scanf("%u%u%u%u%u%u%u%u%u%u%u%u", &p1, &q1, &p2, &q2, &p3, &q3, &p4, &q4, &p5, &q5, &p6, &q6);
	Flush();
	m = 6;
	p = m - 1;
	for (i1 = p1; i1 <= q1; i1++)
	{
	for (i2 = p2; i2 <= q2; i2++)
	{
	for (i3 = p3; i3 <= q3; i3++)
	{
	for (i4 = p4; i4 <= q4; i4++)
	{
	for (i5 = p5; i5 <= q5; i5++)
	{
	for (i6 = p6; i6 <= q6; i6++)
	{
		x = GCDm((USL)i1, (USL)i2);
		y = GCDm(x, (USL)i3);
		z = GCDm(y, (USL)i4);
		u = GCDm(z, (USL)i5);
		v = GCDm(u, (USL)i5);
		w = GCDm(v, (USL)i6);
		if (v == 1)
		{
		printf("(i1,i2,i3,i4,i5,i6)=(%u,%u,%u,%u,%u,%u)\n", i1,i2,i3,i4,i5,i6);
		MATI1 = BUILDMATI(m, 1);
		MATI1->V[0][0] = CHANGE((USL) i1);
		MATI1->V[1][0] = CHANGE((USL) i2);
		MATI1->V[2][0] = CHANGE((USL) i3);
		MATI1->V[3][0] = CHANGE((USL) i4);
		MATI1->V[4][0] = CHANGE((USL) i5);
		MATI1->V[5][0] = CHANGE((USL) i6);
		BB = LLLGCD(MATI1, &GCD, m, m1, n1);
		FREEMATI(MATI1);
		FREEMPI(GCD);
		MATI3 = BUILDMATI(1, m);
		for (l = 0; l < m; l++)
			MATI3->V[0][l] = COPYI(BB->V[m - 1][l]); 
		T1 = LENGTHSQRI(MATI3, 0);
		Q = BB;
		XX = (MPI **)mmalloc((USL)(p * sizeof(MPI *)));
		for (t = 0; t < p; t++)
			XX[t] = ZEROI();
		while (1)
		{
			M = SHORTESTT0(Q, &X);
			for (t = 0; t < p; t++)
			{
				Temp = XX[t];
				XX[t] = ADDI(XX[t], X[t]);
				FREEMPI(Temp);
				FREEMPI(X[t]);
			}
			ffree((char *)X, p * sizeof(MPI *));
			if (M == NULL)
				break;
			else
			{
				for (l = 0; l < Q->C; l++)
				{
					FREEMPI(Q->V[m - 1][l]);
					Q->V[m - 1][l] = COPYI(M->V[0][l]);
				}
				FREEMATI(MATI3);
				MATI3 = M;
			}
		}
		T2 = LENGTHSQRI(MATI3, 0);
		if (RSV(T2, T1) == -1)
		{
			outfile = fopen(buff, "a");
			fprintf(outfile, "(%u,%u,%u,%u,%u,%u): ", i1,i2,i3,i4,i5,i6);
			fprintf(outfile, "b[%u]", m);
			for (t = 0; t < p; t++)
			{
				e = XX[t]->S;
				if (e == -1)
				{
					fprintf(outfile, "+");
					Temp = MINUSI(XX[t]);
					if (!EQONEI(Temp))
						FPRINTI(outfile, Temp);
					fprintf(outfile, "b[%u]", t + 1);
					FREEMPI(Temp);
				}
				if (e == 1)
				{
					fprintf(outfile, "-");
					if (!EQONEI(XX[t]))
						FPRINTI(outfile, XX[t]);
					fprintf(outfile, "b[%u]", t + 1);
				}
			}
			fprintf(outfile, "\n");
			fclose(outfile);
		}
		for (t = 0; t < p; t++)
			FREEMPI(XX[t]);
		ffree((char *)XX, p * sizeof(MPI *));
		FREEMATI(BB);
		FREEMATI(MATI3);
		FREEMPI(T1);
		FREEMPI(T2);
		}
	}
	}
	}
	}
	}
	}
	return;
}

MPMATI *LLLGCD0(MPMATI *DD, MPI **Aptr, USI m, USI m1, USI n1)
/* 
 * Input: an m x 1 vector of positive MPI's.
 * Output: *Aptr = gcd of the vector of MPI's. Also we return a small set of
 * multipliers using the LLL method of Havas and Matthews.
 * Normally m1=1,n1=1.
 * Also returns matrix B of the LLL extended gcd algorithm of 
 * Havas-Majewski-Matthews is returned.
 * S is the shortest length squared all the m short vectors returned.
 * 30/1/97.
 */
{
	unsigned int i, j, k, l, flag, p;
	MPI **A, **D, *X, *Y, *Z, *Tmp, *R, *M1, *N1;
	MPI **XX, *temp, *mu, *SUM, *SHORTER;
	MPI *temp1, *T, *tempI2, *Temp;
	MPMATI *L, *B, *MATI1, *MATI2, *MATI;
	MPR *tempR1, *tR1, *tR2, *tR3, *SUMR, **RHO;
	int e, r, kk, K, c, COUNT, count;
	char buff[20];
	FILE *outfile;
	
	B = IDENTITYI(m);
	A = (MPI **)mmalloc(m * sizeof(MPI *));
	for (i = 0; i < m; i++)
		A[i] = COPYI(DD->V[i][0]);
	D = (MPI **)mmalloc((1 + m) * sizeof(MPI *));
	for (i = 0; i <= m; i++)
		D[i] = ONEI();

	L = ZEROMNI(m, m);
	M1 = CHANGE(m1);
	N1 = CHANGE(n1);

	k = 2;
	while (k <= m)
	{
		REDUCE(k, k - 1, &L, &B, D, A);
		X = MULTI(D[k - 2], D[k]);
		Y = MULTI(D[k - 1], D[k - 1]);
		Tmp = Y;
		Y = MULTI(Y, M1);
		FREEMPI(Tmp);
		R = MULTI(L->V[k - 1][k - 2], L->V[k - 1][k - 2]);
		Z = ADD0I(X, R);
		FREEMPI(R);
		Tmp = Z;
		Z = MULTI(Z, N1);
		FREEMPI(Tmp);
		if (!EQZEROI(A[k-2]) || (EQZEROI(A[k-1]) && (RSV(Y, Z) == 1)))
		{
			for (i = k + 1; i <= m; i++)
				SWAP1(i, k, &L, D);
			SWAP2(k, &B, &L, A);
			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
		{ 
			FREEMPI(X);
			FREEMPI(Y);
			for (l = k - 2; l >= 1; l--)
				REDUCE(k, l, &L, &B, D, A);
			k++;
		}
		FREEMPI(Z);		
	}
	FREEMPI(M1);
	FREEMPI(N1);
	*Aptr = COPYI(A[m - 1]);
	if ((*Aptr)->S == -1)
	{	
		LLLGCDMINUS(&B, &L, m - 1);
		(*Aptr)->S = -((*Aptr)->S);
	}
	printf("L = \n");
	PRINTMATI(0,L->R-1,0,L->C-1,L);
	for (i = 0; i <= m; i++)
	{
		printf("D[%u] = ", i);PRINTI(D[i]);printf(", ");
	}
	printf("\n");
	GetReturn();

	RHO = (MPR **)mmalloc((m - 1) * sizeof(MPR *));
	for (K = m - 2; K >= 0; K--)
	{
		SUMR = RATIOI(L->V[m - 1][K], D[K + 1]);
		for (r = m-2; r > K; r--)
		{
			tR1 = SUMR;
			tR2 = RATIOI(L->V[r][K], D[K + 1]);
			tR3 = MULTR(tR2, RHO[r]);
			FREEMPR(tR2);
			SUMR = ADDR(SUMR, tR3);
			FREEMPR(tR1);
			FREEMPR(tR3);
		}
		RHO[K] = MINUSR(SUMR);
		FREEMPR(SUMR);
	}
	strcpy(buff, "lllgcd0mult.out");
	outfile = fopen(buff, "w");
	for (K = 0; K < m - 1; K++)
	{
		fprintf(outfile, "RHO[%d]=", K + 1);
		printf("RHO[%d]=", K + 1);
		PRINTR(RHO[K]);
		FPRINTR(outfile, RHO[K]);
		printf("\n");
		fprintf(outfile, "\n");
	}

	for (K = 0; K < m - 1; K++)
	{
		fprintf(outfile, "RHO[%d]=", K + 1);
		printf("RHO[%d]=", K + 1);
		PRINTDR(2, RHO[K]);
		FPRINTDR(outfile, 2, RHO[K]);
		FREEMPR(RHO[K]);
		fprintf(outfile, "\n");
		printf("\n");
	}
	ffree((char *)RHO, (m - 1)* sizeof(MPR *));
	GetReturn();
	XX = (MPI **)mmalloc(m * sizeof(MPI *));
	for (i = 0; i < m - 1; i++)
		XX[i] = ZEROI();
	XX[m-1] = ONEI();
	COUNT = m - 2;
	SHORTER = ZEROI();
	MATI = BUILDMATI(1, m);
	for (j = 0; j < m; j++)
		MATI->V[0][j] = ZEROI();
	for (K = m - 2; K >= 0; K--)
	{
		fprintf(outfile, "X[%d]= ", K + 1);
		printf("X[%d]=", K + 1);
		flag = 1;
		for (k = 0; k < m - 1; k++)
		{
			FREEMPI(XX[k]);
			XX[k] = ZEROI();
		}
		for (kk = K; kk >= 0; kk--)
		{
			mu = COPYI(L->V[m - 1][kk]);
			if (flag)
			{
				flag = 0;
				if (mu->S == 0)
				{
					FREEMPI(mu);
					continue;
				}
				else
				{
/*
					tempI1 = ADDI(L->V[m - 1][kk], L->V[m - 1][kk]);
					tempI2= ABSI(tempI1);
					t = EQUALI(tempI2, D[kk]);
					FREEMPI(tempI1);
					FREEMPI(tempI2);
					if (t)
					{
						FREEMPI(mu);
						continue;
					}
*/
					temp = XX[kk];
					if (mu->S == 1)
						XX[kk] = MINUS_ONEI();
					else
						XX[kk] = ONEI();
					FREEMPI(temp);
				}
			}
			else
			{
				SUM = COPYI(mu);
				for (r = m-2; r >= kk+1; r--)
				{
					temp = SUM;
					tempI2 = MULTI(XX[r], L->V[r][kk]);
					SUM = ADDI(SUM, tempI2);
					FREEMPI(tempI2);
					FREEMPI(temp);
				}
				tempR1 = RATIOI(SUM, D[kk + 1]);
				FREEMPI(SUM);
				temp1 = ABS_NEAREST_INTR(tempR1);
/*
				if (K == 15 && kk == 14)
				{
					printf("tempR1=");PRINTR(tempR1);printf("\n");
					printf("temp1=");PRINTI(temp1);printf("\n");
					GetReturn();
					temp = temp1;
					temp1 = ZEROI();
					FREEMPI(temp);
				}
*/
				FREEMPR(tempR1);
				temp = XX[kk];
				XX[kk] = MINUSI(temp1);
				FREEMPI(temp1);
				FREEMPI(temp);
			}
			FREEMPI(mu);
		}
			MATI1 = BUILDMATI(1, m);
			for (j = 0; j < m; j++)
				MATI1->V[0][j] = COPYI(XX[j]);
			MATI2 = MULTMATI(MATI1, B);
			for (i = 0; i < m; i++)
			{
				PRINTI(MATI2->V[0][i]); printf(" ");
				FPRINTI(outfile, MATI2->V[0][i]); fprintf(outfile, " ");
				
			}
			printf("=b[%u]", m);
			fprintf(outfile, "=b[%u]", m);
			p = m - 1;
			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);
				}
			}
			printf(": ");
			fprintf(outfile, ": ");
			T = LENGTHSQRI(MATI2, 0);
			PRINTI(T);
			FPRINTI(outfile, T);
			printf("\n");
			fprintf(outfile, "\n");
			if (K == m - 2)
			{
				FREEMPI(SHORTER);
				SHORTER = COPYI(T);
				FREEMATI(MATI);
				MATI = COPYMATI(MATI2);
			}
			else
			{
				c = RSV(T, SHORTER);
				if (c < 0)
				{