lllgcd.c

calc大数库

/* lllgcd.c */
#ifdef _WIN32
#include "unistd_DOS.h"
#else
#include <unistd.h>
#endif
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "integer.h"
#include "fun.h"
extern MPI *MAXI, *PMAXI;

extern unsigned int MLLLVERBOSE;
extern unsigned int GCDVERBOSE;
USI GCD3FLAG;
extern unsigned int GCD_MAX;

MPMATI *LLLGCD(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.
 * matrix B of the LLL extended gcd algorithm of 
 * Havas-Majewski-Matthews is returned.
 * 30/1/97.
 */
{
	unsigned int i, 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);
	}
/*
	if (GCDVERBOSE)
	{
		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");
		temp1  = MULTI(D[1], D[1]);
		printf("a = "); PRINTI(temp1); printf(", ");
		FREEMPI(temp1);
		temp1 = MULTI(D[1], L->V[1][0]);
		temp2 = MULT_I(temp1, 2);
		printf("2h = "); PRINTI(temp2); printf(", ");
		FREEMPI(temp1);
		FREEMPI(temp2);
		temp1 = MULTI(L->V[1][0], L->V[1][0]);
		temp2 = ADD0I(D[2], temp1);
		printf("b = "); PRINTI(temp2); printf(", ");
		FREEMPI(temp1);
		FREEMPI(temp2);
		temp1 = MULTI(D[1], L->V[2][0]);
		temp2 = MULT_I(temp1, 2);
		printf("2g = "); PRINTI(temp2); printf(", ");
		FREEMPI(temp1);
		FREEMPI(temp2);
		temp1 = MULTI(L->V[1][0], L->V[2][0]);
		FREEMPI(temp1);
		temp2 = ADDI(temp1, L->V[2][1]);
		temp1 = MULT_I(temp2, 2);
		printf("2f = "); PRINTI(temp1); printf(", ");
		FREEMPI(temp1);
		FREEMPI(temp2);
	}
*/
/*
	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 SWAP1(USI i, USI k, MPMATI **Lptr, MPI *D[])
{
	MPI *X1, *X2, *X3, *Y1, *Y2, *Tmp;

	X1 = MULTI((*Lptr)->V[i - 1][k - 2], (*Lptr)->V[k - 1][k - 2]);
	Y1 = MULTI((*Lptr)->V[i - 1][k - 1], D[k - 2]);
	Tmp = Y1;
	Y1 = ADDI(Y1, X1);
	FREEMPI(Tmp);
	FREEMPI(X1);
	X2 = MULTI((*Lptr)->V[i - 1][k - 2], D[k]);
	X3 = MINUSI((*Lptr)->V[k - 1][k - 2]);
	Y2 = MULTI((*Lptr)->V[i - 1][k - 1], X3);
	FREEMPI(X3);
	Tmp = Y2;
	Y2 = ADDI(Y2, X2);
	FREEMPI(Tmp);
	FREEMPI(X2);
	FREEMPI((*Lptr)->V[i - 1][k - 2]);
	(*Lptr)->V[i - 1][k - 2] = INT(Y1, D[k - 1]);
	FREEMPI((*Lptr)->V[i - 1][k - 1]);
	(*Lptr)->V[i - 1][k - 1] = INT(Y2, D[k - 1]);
	FREEMPI(Y1);
	FREEMPI(Y2);
	return;
}


void REDUCE(USI k, USI l, MPMATI **Lptr, MPMATI **Bptr, MPI *D[], MPI *A[])
/*
 * updates *Lptr, *Bptr.
 */
{
	unsigned int i, j, m, n;
	MPI *X, *Y, *Q, *Tmp, *Z;
	MPMATI *TmpMATI, *TempMAT;

	m = (*Bptr)->C;
	n = (*Bptr)->R;
	if (A[l - 1]->S)
	{
		Q = NEAREST_INTI(A[k - 1], A[l - 1]);
		if (Q->S)
		{
			Tmp = A[k - 1];
			Z = MULTI(A[l - 1], Q);
			A[k - 1] = SUBI(A[k - 1],Z);
			FREEMPI(Tmp);
			FREEMPI(Z);
		}
	}
	else
	{
		Y = MULT_I((*Lptr)->V[k - 1][l - 1], 2);
		if (RSV(Y, D[l]) == 1)
			Q = NEAREST_INTI((*Lptr)->V[k - 1][l - 1], D[l]);
		else
			Q = ZEROI();
		FREEMPI(Y);
	}
	if (Q->S == 0)
	{
		FREEMPI(Q);
		return;
	}
	X = MINUSI(Q);
	TmpMATI = *Bptr;
	*Bptr = ADD_MULT_ROWI(l - 1, k - 1, X, *Bptr);
	FREEMATI(TmpMATI);
	if (MLLLVERBOSE)
	{
		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);
	}
	if(GCDVERBOSE)
	{
		printf("Row %u -> Row %u + ", k,k);PRINTI(X);printf(" x Row %u\n", l);
		TempMAT = BUILDMATI(n, n + 1);
		for (i = 0; i < n; i++)
		{
			for (j = 0; j <= n; j++)
			{		
				if (j == n)
					TempMAT->V[i][j] = COPYI(A[i]);
				else
					TempMAT->V[i][j] = COPYI((*Bptr)->V[i][j]);
			}
		}
		PRINTMATI(0,TempMAT->R-1,0,TempMAT->C-1,TempMAT);
		printf("\n");
		FREEMATI(TempMAT);
	}

	FREEMPI(X);
	for (j = 1; j < l; j++)
	{
		X = MULTI((*Lptr)->V[l - 1][j - 1], Q);
		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], Q);
	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(Q);
}


void SWAP2(USI k, MPMATI **B1ptr, MPMATI **Lptr, MPI *A[])
{
	MPI *T, *S;
	unsigned int i, j, n;
	MPMATI *TempMAT;

	n = (*B1ptr)->R;
	S = A[k - 2];
	A[k - 2] = A[k - 1];
	A[k - 1] = S;
	
	*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();
	if(GCDVERBOSE)
	{
		printf("Swapping Rows %u and %u\n", k-1,k);
		TempMAT = BUILDMATI(n, n + 1);
		for (i = 0; i < n; i++)
		{
			for (j = 0; j <= n; j++)
			{		
				if (j == n)
					TempMAT->V[i][j] = COPYI(A[i]);
				else
					TempMAT->V[i][j] = COPYI((*B1ptr)->V[i][j]);
			}
		}
		PRINTMATI(0,TempMAT->R-1,0,TempMAT->C-1,TempMAT);
		printf("\n");
		GetReturn();
		FREEMATI(TempMAT);
	}
	if (MLLLVERBOSE)
	{
		printf("Swapping Rows %u and %u\n", k-1,k);
		PRINTMATI(0,(*B1ptr)->R-1,0,(*B1ptr)->C-1,*B1ptr);
	}
	return;
}

MPMATI *JACOBIGCD(MPMATI *DD, MPI **Aptr, USI m)
/* 
 * 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 a version of a method of Jacobi
 * A unimodular transforming matrix B is returned.
 * 31/1/97.
 */
{
	unsigned int i, z;
	int k;
	MPI **A, *Q, *X, *Tmp, **temp;
	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]);
	z = 0;
	while (1)
	{
 		/* A[i] = 0 for i = 0,...,z-1, but A[z] != 0. */
		for (i = z + 1; i < m; i++)
		{
			Q = INT0(A[i], A[z]);
			X = MINUSI(Q);
			FREEMPI(Q);
			L = B;
			B = ADD_MULT_ROWI(z, i, X, B);
			FREEMPI(X);
			FREEMATI(L);
			Tmp = A[i];
			A[i] = MOD0(A[i], A[z]);
			FREEMPI(Tmp);
		}
	if (MLLLVERBOSE)
	{
			printf("After the reduction: B = \n");
			PRINTMATI(0,B->R-1,0,B->C-1,B);
			printf("\n");
			for (i = 0; i < m; i++)
			{
				printf("A[%u] = ", i); PRINTI(A[i]); printf("\n");
			}
				GetReturn();
	}
		Tmp = A[z];
		temp = B->V[z];
		for (k = z; k <= m - 2; k++)
		{
			A[k] = A[k + 1];
			B->V[k] = B->V[k + 1];
		}
		A[m - 1] = Tmp;
		B->V[m - 1] = temp;
	if (MLLLVERBOSE)
	{
		printf("After the rearrangement: B = \n");
		PRINTMATI(0,B->R-1,0,B->C-1,B);
		printf("\n");
			for (i = 0; i < m; i++)
			{
				printf("A[%i] = ", i); PRINTI(A[i]); printf("\n");
			}
			GetReturn();
	}
		while (A[z]->S  == 0)
			z++;
		if (z == m - 1)
			break;
	}
	
	*Aptr = COPYI(A[m - 1]);
	if ((*Aptr)->S == -1)
	{
		(*Aptr)->S = 1;	
		for (i = 0; i < m; i++)
			 (B->V[m - 1][i])->S = -((B->V[m - 1][i])->S);	
	}
	for (i = 0; i < m; i++)
		FREEMPI(A[i]);
	ffree((char *)A, m * sizeof(MPI *));
	return (B);
}

void GCD4()
/* We compute an optimal multipliers for a[1],a[2],a[3],a[4]
 * in the range p1<=a[1]<=q1 etc. and rel prime.
 */
{
	USI l, m, p1, q1, p2, q2, p3, q3, p4, q4, m1, n1, t;
	USI i1, i2, i3, i4, p;
	int e;
	USL x, y, z;
	MPI *GCD, *T1, *T2, **XX, **X, *Temp, *DIFF;
	MPMATI *MATI1, *BB, *MATI3, *Q, *M;
	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, "gcd4.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", &p1, &q1, &p2, &q2, &p3, &q3, &p4, &q4);
	Flush();
	m = 4;
	p = m - 1;
	for (i1 = p1; i1 <= q1; i1++)
	{
	/*for (i2 = (i1+1 > p2 ? i1+1: p2); i2 <= q2; i2++)*/
	for (i2 = p2; i2 <= q2; i2++)
	{
	/*for (i3 = (i2+1 > p3 ? i2+1: p3); i3 <= q3; i3++)*/
	for (i3 = p3; i3 <= q3; i3++)
	{
	/*for (i4 = (i3+1 > p4 ? i3+1: p4); i4 <= q4; i4++)*/
	for (i4 = p4; i4 <= q4; 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);
		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);
		if (RSV(T2, T1) == -1)
		{
			outfile = fopen(buff, "a");
			fprintf(outfile, "(%u,%u,%u,%u): ", i1,i2,i3,i4);
			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 GCD5()
/* We compute an multipliers for a[1],a[2],a[3],a[4],a[5] 
 * in the range p1<=a[1]<=q1 etc. and rel prime.
 */
{
	USI l, m, p1, q1, p2, q2, p3, q3, p4, q4, p5, q5, m1, n1;
	USI i1, i2, i3, i4, i5, p, t;
	int e;
	USL x, y, z, u;
	MPI *GCD, *T1, *T2, **XX, **X, *Temp, *DIFF;
	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++)*/

	printf("Enter alpha=m/n,  1/4 < m/n <= 1:  (ie enter  m  and n)");
	scanf("%u%u", &m1, &n1);
	strcpy(buff, "gcd5.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++)