lllgcd.c

calc大数库

	{
		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");
*/
	
	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--)
	{
		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);
				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);
			p = m - 1;
			T = LENGTHSQRI(MATI2, 0);
			if (K == m - 2)
			{
				FREEMPI(SHORTER);
				SHORTER = COPYI(T);
				FREEMATI(MATI);
				MATI = COPYMATI(MATI2);
			}
			else
			{
				c = RSV(T, SHORTER);
				if (c < 0)
				{
					FREEMPI(SHORTER);
					SHORTER = COPYI(T);
					FREEMATI(MATI);
					MATI = COPYMATI(MATI2);
					COUNT = K;
				}
			}
			FREEMPI(T);
			FREEMATI(MATI1);
			FREEMATI(MATI2);
	}
	T = LENGTHSQRI(B, m - 1);
	c = RSV(T, SHORTER);
	if (c < 0)
	{
		FREEMPI(SHORTER);
		SHORTER = COPYI(T);
		for (i = 0; i < m; i++)
		{
			FREEMPI(MATI->V[0][i]);
			MATI->V[0][i] = COPYI(B->V[m-1][i]);
		}
		COUNT = K;
	}
	FREEMPI(T);
	count = COUNT + 1;
	for (i = 0; i < m; i++)
	{
		FREEMPI(B->V[m-1][i]);
		B->V[m-1][i] = COPYI(MATI->V[0][i]);
        }

	FREEMATI(MATI);
	FREEMPI(SHORTER);
	FREEMATI(L);
	for (i = 0; i <= m; i++)
		FREEMPI(D[i]);
	ffree((char *)D, (1 + m) * sizeof(MPI *));
	for (i = 0; i < m; i++)
		FREEMPI(A[i]);
	ffree((char *)A, m * sizeof(MPI *));
	for (i = 0; i < m; i++)
		FREEMPI(XX[i]);
	ffree((char *)XX, m * sizeof(MPI *));
	return (B);
}

MPMATI *LLLGCDL(MPMATI *DD, MPI **Aptr, USI m, USI m1, USI n1, MPMATR **LMATRIX)
/* 
 * 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.
 * Also output the L=[mu[i][j]] matrix.
 * matrix B of the LLL extended gcd algorithm of 
 * Havas-Majewski-Matthews is returned.
 */
{
	unsigned int i, j, k, l;
	MPI **A, **D, *X, *Y, *Z, *Tmp, *R, *M1, *N1;
	MPMATI *L, *B;

	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);
	}
	*LMATRIX = ZEROMNR(m , m);
	for (i = 0; i < m; i++)
		for (j = 0; j < m; j++)
		{
			if (i > j)
			{
				FREEMPR(elt(*LMATRIX, i, j));
				elt(*LMATRIX, i, j) = RATIOI(L->V[i][j], D[j + 1]);
			}
		
		}
/*
	printf("*LMATRIX = \n");
	PRINTMATR(0,(*LMATRIX)->R-1,0,(*LMATRIX)->C-1,*LMATRIX);
	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");
*/
	FREEMATI(L);
	for (i = 0; i <= m; i++)
		FREEMPI(D[i]);
	for (i = 0; i < m; i++)
		FREEMPI(A[i]);
	ffree((char *)D, (1 + m) * sizeof(MPI *));
	ffree((char *)A, m * sizeof(MPI *));
	return (B);
}

void GCDCONJECTURE5()
/* We test our LLLGCD conjecture for a[1],a[2],a[3],a[4],a[5] 
 * in the range p1<=a[1]<=q1 etc. and rel prime.
 */
{
	USI p1, q1, p2, q2, p3, q3, p4, q4, p5, q5;
	USI i1, i2, i3, i4, i5, p, record=1;
	USL x, y, z, u;
	USI m1, n1, m, n, i, j, count;
	int r, e, K;
	MPMATI *MATI1, *MATI2, *Q, *M, *MATI3;
	MPI *A, **XX, **X, *Temp, ***COEFF, *tempI;
	MPMATR *LMATRIX;
	MPR **SIGMA, *SUMR, *tR1, *tR2, *tR3, *tempR, *tempR1;
	char buff[20];
	FILE *outfile;

	printf("Enter alpha=m/n,  1/4 < m/n <= 1:  (ie enter  m  and n)");
	scanf("%u%u", &m1, &n1);
	strcpy(buff, "conjecture5.out");
	if(access(buff, R_OK) == 0)
	  unlink(buff);
	printf("Enter  the ranges pi,qi; 2 <= pi < qi < 2^32, i = 1,...,5: ");
	scanf("%u%u%u%u%u%u%u%u%u%u", &p1, &q1, &p2, &q2, &p3, &q3, &p4, &q4, &p5, &q5);
	Flush();
	m = 5;
	p = m - 1;
	for (i1 = p1; i1 <= q1; i1++)
	{
	for (i2 = p2; i2 <= q2; i2++)
	{
	for (i3 = p3; i3 <= q3; i3++)
	{
	printf("(i1,i2,i3)=(%u,%u,%u)\n", i1,i2,i3);
	for (i4 = p4; i4 <= q4; i4++)
	{
	for (i5 = p5; i5 <= q5; i5++)
	{
/*printf("(i1,i2,i3,i4,i5)=(%u,%u,%u,%u,%u)\n", i1,i2,i3,i4,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)
		{
			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);

			m = MATI1->R;
			MATI2 = LLLGCDL(MATI1, &A, m, m1, n1, &LMATRIX);
			FREEMATI(MATI1);
			FREEMPI(A);
			MATI3 = BUILDMATI(1, m);
			for (i = 0;i < m;i++)
				MATI3->V[0][i] = COPYI(MATI2->V[m - 1][i]); 
			p = m - 1;
			XX = (MPI **)mmalloc((USL)(p * sizeof(MPI *)));
			for (j = 0; j < p; j++)
				XX[j] = ZEROI();

			Q = MATI2;
			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(MATI3);
					MATI3 = M;
				}
			}

			/*SHORTEST(Q, XX, 2, 0);*/
			COEFF = SHORTESTX(Q, XX, &count);
if (count>record)
{
			record=count;
			outfile = fopen(buff, "a");
			printf("(i1,i2,i3,i4,i5)=(%u,%u,%u,%u,%u): record=%u\n", i1,i2,i3,i4,i5, record);
			fprintf(outfile, "(i1,i2,i3,i4,i5)=(%u,%u,%u,%u,%u): record=%u\n", i1,i2,i3,i4,i5, record);
			fclose(outfile);
}
			for (j = 0; j < p; j++)
				FREEMPI(XX[j]);
			ffree((char *)XX, p * sizeof(MPI *));
			SIGMA = (MPR **)mmalloc(p * sizeof(MPR *));
			for (j = 0; j < count; j++)
			{
				/*printf("SIGMA[%u]=", j);*/
				for (K = p - 1; K >= 0; K--)
				{
					SUMR = COPYR(elt(LMATRIX, m - 1, K));
					for (r = p - 1; r > K; r--)
					{
						tR2 = COPYR(elt(LMATRIX, r, K));
						tempR = BUILDMPR();
						tempR->N = COPYI(COEFF[j][r]);
						tempR->D = ONEI();
						tR3 = MULTR(tR2, tempR);
						FREEMPR(tempR);
						FREEMPR(tR2);
						tR1 = SUMR;
						SUMR = ADDR(SUMR, tR3);
						FREEMPR(tR1);
						FREEMPR(tR3);
					}
					SIGMA[K] = SUMR;
					/*PRINTR(SIGMA[K]);
					printf(" ");*/
					tempR = BUILDMPR();
					tempR->N = COPYI(COEFF[j][K]);
					tempR->D = ONEI();
					tempR1 = ADDR(tempR, SIGMA[K]);
					tempI = ABSI(tempR1->N);
					e = RSV(tempI, tempR1->D);
					FREEMPR(tempR);
					FREEMPR(tempR1);
					FREEMPI(tempI);
					if (e >= 0)
					{
						fprintf(stderr, "conjecture false:j = %u, K = %d\n", j,  K + 1);
						outfile = fopen(buff, "a"), count;
		fprintf(outfile, "(i1,i2,i3,i4,i5)=(%u,%u,%u,%u,%u): ", i1,i2,i3,i4,i5);
						fprintf(outfile, "conjecture false: j = %u, K = %d\n", j, K + 1);
						fclose(outfile);
						exit (1);
					}
				}
				/*printf("\n");*/
				for (K = 0; K < p; K++)
					FREEMPR(SIGMA[K]);
			}
			ffree((char *)SIGMA, p * sizeof(MPR *));

			ffree((char *)COEFF, GCD_MAX * sizeof(MPI **));
			for (j = 0; j < count; j++)
			{
				ffree((char *)COEFF[j], p * sizeof(MPI *));
				for (i = 0; i < p; i++)
					FREEMPI(COEFF[j][i]);
			}
			FREEMATI(MATI3);
			FREEMATI(Q);
			FREEMATR(LMATRIX);
		}
	}
	}
	}
	}
	}
	return;
}

void GCDCONJECTURE4()
/* We test our LLLGCD conjecture for a[1],a[2],a[3],a[4].
 * in the range p1<=a[1]<=q1 etc. and rel prime.
 */
{
	USI p1, q1, p2, q2, p3, q3, p4, q4;
	USI i1, i2, i3, i4, p, record=1;
	USL x, y, z;
	USI m1, n1, m, n, i, j, count;
	int r, e, K;
	MPMATI *MATI1, *MATI2, *Q, *M, *MATI3;
	MPI *A, **XX, **X, *Temp, ***COEFF, *tempI;
	MPMATR *LMATRIX;
	MPR **SIGMA, *SUMR, *tR1, *tR2, *tR3, *tempR, *tempR1;
	char buff[20];
	FILE *outfile;

	printf("Enter alpha=m/n,  1/4 < m/n <= 1:  (ie enter  m  and n)");
	scanf("%u%u", &m1, &n1);
	strcpy(buff, "conjecture4.out");
	if(access(buff, R_OK) == 0)
	  unlink(buff);
	printf("Enter  the ranges pi,qi; 2 <= pi < qi < 2^32, i = 1,..4: ");
	scanf("%u%u%u%u%u%u%u%u", &p1, &q1, &p2, &q2, &p3, &q3, &p4, &q4);
	Flush();
	m = 4;
	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++)
	{
		printf("(i1,i2,i3,i4)=(%u,%u,%u,%u)\n", i1,i2,i3,i4);
		x = GCDm((USL)i1, (USL)i2);
		y = GCDm(x, (USL)i3);
		z = GCDm(y, (USL)i4);
		if (z == 1)
		{
			printf("(i1,i2,i3,i4)=(%u,%u,%u,%u)\n", i1,i2,i3,i4);
			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);

			m = MATI1->R;
			MATI2 = LLLGCDL(MATI1, &A, m, m1, n1, &LMATRIX);
			FREEMATI(MATI1);
			FREEMPI(A);
			MATI3 = BUILDMATI(1, m);
			for (i = 0;i < m;i++)
				MATI3->V[0][i] = COPYI(MATI2->V[m - 1][i]); 
			p = m - 1;
			XX = (MPI **)mmalloc((USL)(p * sizeof(MPI *)));
			for (j = 0; j < p; j++)
				XX[j] = ZEROI();

			Q = MATI2;
			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]);
					}