nfunc.c

calc大数库

	printf("The matrix entered is\n\n");
	PRINTMATI(0,MATI1->R-1,0,MATI1->C-1,MATI1);
	printf("enter the rows to be sorted :");
	scanf("%u %u", &r, &s);
	QSORTMATI(MATI1, r, s);
	printf("The sorted matrix is\n");
	PRINTMATI(0,MATI1->R-1,0,MATI1->C-1,MATI1);
	FREEMATI(MATI1);
	return;
}

void SCHNORRGCD(MPI *N)
{
	USI answer, m1, n1, n, t, i, j;
	MPMATI *MATI1, *MATI2, *MATI0;
	char buff[20];
	FILE *outfile;

	HAVASFLAG = 1;
	printf("Do you wish to enter your sequence of numbers from an existing column matrix? (Y/N)\n");
	answer = GetYN();
	if (answer)
	{
		MATI1 = FINPUTMATFILEI_I();
		if (MATI1 == NULL)
			exit (0);
	}
	else
	{
		printf("WARNING: Make sure the first integer in the sequence is non-zero) :\n");
		MATI1 = INPUTMATI();
	}
	printf("The matrix entered is\n\n");
	PRINTMATI(0,MATI1->R-1,0,MATI1->C-1,MATI1);
	MLLLVERBOSE = 0;
	printf("VERBOSE? (Y/N)\n");
	answer = GetYN();
	printf("answer = %d\n", answer);
	if (answer)
		MLLLVERBOSE = 1;
	printf("MLLLVERBOSE = %u\n", MLLLVERBOSE);
	printf("to enter alpha=m1/n1: select m1 and n1 (normally 1 and 1) :");
	scanf("%u %u", &m1, &n1);
	Flush();
	n = MATI1->R;
	t = n + 1;
	MATI0 = BUILDMATI(n, t);
	for (i = 0; i < n; i++)
	{
		for (j = 0; j < n; j++)
		{
			if ( i == j)
				MATI0->V[i][j] = ONEI();
			else
				MATI0->V[i][j] = ZEROI();
		}
		MATI0->V[i][n] = MULTI(MATI1->V[i][0], N);
	}
/*
	MATI0 = BUILDMATI(t, t + 1);
	for (i = 0; i < n; i++)
	{
		for (j = 0; j < n; j++)
		{
			if ( i == j)
				MATI0->V[i][j] = ONEI();
			else
				MATI0->V[i][j] = ZEROI();
		}
		MATI0->V[i][n] = MULTI(MATI1->V[i][0], N);
		MATI0->V[i][t] = ZEROI();
	}
	for (j = 0; j < n; j++)
		MATI0->V[n][j] = ONEI();
*/
	FREEMATI(MATI1);
	strcpy(buff, "sgcdmat.out");
	outfile = fopen(buff, "w");
	FPRINTMATI(outfile,0,MATI0->R-1,0,MATI0->C-1,MATI0);
	fclose(outfile);
		
GetReturn();
	MATI2 = BASIS_REDUCTION0(MATI0, m1, n1);
	FREEMATI(MATI0);
	printf("The corresponding reduced basis is\n");
	PRINTMATI(0,MATI2->R-1,0,MATI2->C-1,MATI2);
	strcpy(buff, "sgcdbas.out");
	outfile = fopen(buff, "w");
	FPRINTMATI(outfile,0,MATI2->R-1,0,MATI2->C-1,MATI2);
	fclose(outfile);
	FREEMATI(MATI2);
	HAVASFLAG = 0;
	return;
}

void SHORTESTTTX()
{
/* This is the inhomogenous Fincke-Pohst Algorithm .
 * See Pohst-Zassenhaus .
 */
	USI answer, u;
	MPI *BOUND;
	MPMATI *MATI1;

	printf("Do you wish to enter your matrix from an existing matrix? (Y/N)\n");
	answer = GetYN();
	if (answer)
	{
		MATI1 = FINPUTMATFILEI_I();
		if (MATI1 == NULL)
			exit (0);
	}
	else
		MATI1 = INPUTMATI();
	printf("The matrix entered is\n\n");
	PRINTMATI(0,MATI1->R-1,0,MATI1->C-1,MATI1);
	MLLLVERBOSE = 0;
	printf("VERBOSE? (Y/N)\n");
	answer = GetYN();
	printf("answer = %d\n", answer);
	if (answer)
		MLLLVERBOSE = 1;
	printf("MLLLVERBOSE = %u\n", MLLLVERBOSE);
	printf("enter a positive integer (the distance bound):");
	BOUND = INPUTI(&u);
	Flush();
	if (BOUND->S <= 0)
	{
	  printf("BOUND <= 0\n");
	  return;
	}
	SHORTESTTT(MATI1, BOUND);
	FREEMATI(MATI1);
	FREEMPI(BOUND);
	return;
}

void FIB_MIN()
{
	USI M1, M2, N1, N2, i, m, n;
	MPMATI *MATI1, *MATI2, *Q;
	MPI *A;
	char buff[20];
	FILE *outfile;

	printf("enter positive integers N1 <= N2, the test range:");
	scanf("%u%u", &N1, &N2);
	printf("enter integers (>= 2) M1 <= M2, the sequence range:");
	scanf("%u%u", &M1, &M2);
	Flush();
	strcpy(buff, "fibmin.out");
	outfile = fopen(buff, "w");
	outfile = fopen(buff, "a+");
	FPRINTMATIFLAG = 1;
	MLLLVERBOSE = 0;
for (m = M1; m <= M2; m++)
{
	fprintf(outfile, "====================\n");
	for (n = N1; n <= N2; n++)
	{
		printf("m=%u: n=%u:\n ", m,  n);
		fprintf(outfile, "m=%u: n=%u:", m,  n);
		MATI1 = BUILDMATI(m, 1);
		for (i = 0; i < MATI1->R; i++)
		{
			MATI1->V[i][0] = FIBONACCI(n + i);	
/*
			FPRINTI(outfile, MATI1->V[i][0]);
			fprintf(outfile, " ");
*/
		}
		MATI2 = EXTGCD(MATI1, &A, &Q, 3, 4);
		FREEMATI(MATI1);
		FREEMATI(MATI2);
		FREEMPI(A);
		FPRINTMATI(outfile, Q->R - 1, Q->R-1, 0, Q->C-1, Q);
		FREEMATI(Q);
	}
}
	fprintf(outfile, "====================\n");
	fclose(outfile);
	FPRINTMATIFLAG = 0;
	return;
}

void PRINTW1(USI n, USI r)
{
	USI x;
	x = n - r - 2;
	if (n % 2 == 0)
	{
		if (r <= n - 4)
			printf(" L[%u]       ", x);
		if (r == n - 3)
			printf(" 2           ");
		if (r == n - 2)
			printf(" 1           ");
		if (r >= n - 1)
			printf(" 0           ");
	}
	if (n % 2)
	{
		if (r <= n - 3)
			printf(" L[%u]       ", x);
		if (r == n - 2)
			printf(" 1           ");
		if (r == n - 1)
			printf(" 1           ");
		if (r >= n)
			printf(" 0           ");
	}
}

void PRINTW2(USI n, USI r, USI m)
{
	USI x,y,z;
	x = n - r - 2; y = 2 * m - n; z = r - y - 3;
	if (n % 2 == 0)
	{
		if (r < m)
		{
			if (r <= y + 1)
				printf(" L[%u]       ", x);
			if (r == y + 2)
				printf(" L[%u]+1     ", x);
			if (r == y + 3)
				printf(" L[%u]-1     ", x);
			if (r == y + 4)
				printf(" L[%u]+2     ", x);
			if (r == y + 5)
				printf(" L[%u]-2     ", x);
			if (r < m && r >= y + 6 && r % 2 == 0)
				printf(" L[%u]+L[%u] ", x,z);
			if (r < m && r >= y + 6 && r % 2)
				printf(" L[%u]-L[%u] ", x,z);
		}
		if (r ==  m && n == m + 3)
			printf(" 1           ");
		if (r ==  m && n == m + 5)
			printf(" 2           ");
		if (r ==  m && n >= m + 6)
			printf(" L[%u]       ",n - m - 4);
	}
	if (n % 2)
	{
		if (r < m)
		{
			if (r <= y)
				printf(" L[%u]       ", x);
			if (r == y + 1)
				printf(" L[%u]+1     ", x);
			if (r == y + 2)
				printf(" L[%u]-1     ", x);
			if (r == y + 3)
				printf(" L[%u]+1     ", x);
			if (r == y + 4)
				printf(" L[%u]-2     ", x);
			if (r < m && r >= y + 5 && r % 2 == 0)
				printf(" L[%u]+L[%u] ", x,z);
			if (r < m && r >= y + 6 && r % 2)
				printf(" L[%u]-L[%u] ", x,z);
		}
		if (r ==  m && n == m + 4)
			printf(" 1           ");
		if (r ==  m && n >= m + 5)
			printf(" L[%u]       ",n - m - 4);
	}
}

 void PRINTWW()
{
	USI m, t, w, x, y, z, n, r;

	printf("enter m:");
	scanf("%u", &m);
	x = m + 3; y = 2 * m, z = x - 2; w = m - 1;
	for (n = 3; n <= z; n++)
	{
		for (r = 1; r <= w; r++)
			PRINTW1(n, r);
		printf(" 0           ");
		printf("\n");
	}
	for (r = 1; r <= m; r++)
	{
		t = m - r;
		if (r < m - 1)
			printf(" L[%u]       ", t);
		else if (r == m - 1)
			printf(" L[%u]+1     ", t);
		else
			printf(" 0           ");
	}
	printf("\n");
	for (n = x; n <= y; n++)
	{
		for (r = 1; r <= m; r++)
			PRINTW2(n, r, m);
		printf("\n");
	}
	return;
}

void PRINT_DEFECT()
{
	USI M1, M2, N1, N2, i, m, n;
	MPMATI *MATI, *MATI1, *MATI2, *MATI3, *Q1, *Q2, *Q3, *Q;
	MPI *A;
	char buff[20];
	FILE *outfile;

	printf("enter positive integers N1 <= N2, the test range:");
	scanf("%u%u", &N1, &N2);
	printf("enter integers (>= 2) M1 <= M2, the sequence range:");
	scanf("%u%u", &M1, &M2);
	Flush();
	strcpy(buff, "defect.out");
	outfile = fopen(buff, "w");
	outfile = fopen(buff, "a+");
	FPRINTMATIFLAG = 1;
	MLLLVERBOSE = 0;
for (m = M1; m <= M2; m++)
{
	printf("m=%u:\n ", m);
	fprintf(outfile, "====================\n");
	for (n = N1; n <= N2; n++)
	{
		printf("n=%u:\n ", n);
		fprintf(outfile, "n=%u:", n);
		MATI1 = BUILDMATI(m, 1);
		MATI2 = BUILDMATI(m, 1);
		MATI3 = BUILDMATI(m, 1);
		for (i = 0; i < MATI1->R; i++)
		{
			MATI1->V[i][0] = FIBONACCI(n + i);	
			MATI2->V[i][0] = FIBONACCI(n + 1 + i);	
			MATI3->V[i][0] = FIBONACCI(n + 2 + i);	
		}
		MATI = EXTGCD(MATI1, &A, &Q1, 3, 4);
		FREEMATI(MATI);
		FREEMATI(MATI1);
		FREEMPI(A);
		MATI = EXTGCD(MATI2, &A, &Q2, 3, 4);
		FREEMATI(MATI);
		FREEMATI(MATI2);
		FREEMPI(A);
		MATI = EXTGCD(MATI3, &A, &Q3, 3, 4);
		FREEMATI(MATI);
		FREEMATI(MATI3);
		FREEMPI(A);
		MATI = ADDMATI(Q3, Q2);
		Q = SUBMATI(MATI, Q1);
		FPRINTMATI(outfile, Q->R - 1, Q->R-1, 0, Q->C-1, Q);
		FREEMATI(Q1);
		FREEMATI(Q2);
		FREEMATI(Q3);
		FREEMATI(MATI);
	}
}
	fprintf(outfile, "====================\n");
	fclose(outfile);
	FPRINTMATIFLAG = 0;
	return;
}

void LUCAS_MIN()
{
	USI M1, M2, N1, N2, i, m, n;
	MPMATI *MATI1, *MATI2, *Q;
	MPI *A;
	char buff[20];
	FILE *outfile;

	printf("enter positive integers N1 <= N2, the test range:");
	scanf("%u%u", &N1, &N2);
	printf("enter integers (>= 2) M1 <= M2, the sequence range:");
	scanf("%u%u", &M1, &M2);
	Flush();
	strcpy(buff, "lucasmin.out");
	outfile = fopen(buff, "w");
	outfile = fopen(buff, "a+");
	FPRINTMATIFLAG = 1;
	MLLLVERBOSE = 0;
for (m = M1; m <= M2; m++)
{
	fprintf(outfile, "====================\n");
	for (n = N1; n <= N2; n++)
	{
		printf("m=%u: n=%u:\n ", m,  n);
		fprintf(outfile, "m=%u: n=%u:", m,  n);
		MATI1 = BUILDMATI(m, 1);
		for (i = 0; i < MATI1->R; i++)
		{
			MATI1->V[i][0] = LUCASS(n + i);	
/*
			FPRINTI(outfile, MATI1->V[i][0]);
			fprintf(outfile, " ");
*/
		}
		MATI2 = EXTGCD(MATI1, &A, &Q, 3, 4);
		FREEMATI(MATI1);
		FREEMATI(MATI2);
		FREEMPI(A);
		FPRINTMATI(outfile, Q->R - 1, Q->R-1, 0, Q->C-1, Q);
		FREEMATI(Q);
	}
}
	fprintf(outfile, "====================\n");
	fclose(outfile);
	FPRINTMATIFLAG = 0;
	return;
}

void LLLGCDX()
/*
 * This executes the LLLGCD algorithm of Havas, Majewski and Matthews.
 * 30/1/97.
 */
{
	USI answer, m1, n1, m, n, i, j, p;
	int e;
	MPMATI *MATI1, *MATI2, *Q, *M, *MATI3;
	char buff[20];
	FILE *outfile;
	MPI *A, *T, **XX, **X, *Temp;

	HAVASFLAG = 1;
	printf("Do you wish to enter your sequence of numbers from an existing column matrix? (Y/N)\n");
	answer = GetYN();
	if (answer)
	{
		MATI1 = FINPUTMATFILEI_I();
		if (MATI1 == NULL)
			exit (0);
	}
	else
	{
		printf("WARNING: Make sure the first integer in the sequence is non-zero) :\n");
		MATI1 = INPUTMATI();
	}
	printf("The matrix entered is\n\n");
	PRINTMATI(0,MATI1->R-1,0,MATI1->C-1,MATI1);
	GCDVERBOSE = 0;
	printf("VERBOSE? (Y/N)\n");
	answer = GetYN();
	if (answer)
		GCDVERBOSE = 1;
	printf("to enter alpha=m1/n1: select m1 and n1 (normally 1 and 1) :");
	scanf("%u %u", &m1, &n1);
	Flush();
	m = MATI1->R;
	MATI2 = LLLGCD(MATI1, &A, m, m1, n1);
	FREEMATI(MATI1);
	printf("The multiplier matrix found by LLL is\n");
	PRINTMATI(0,MATI2->R-1,0,MATI2->C-1,MATI2);
	printf("gcd = "); PRINTI(A); printf("\n");
	FREEMPI(A);
	printf("The multiplier vector found by LLL is b[%u]=", m);
	for (i = 0;i < MATI2->R;i++)
	{
		PRINTI(MATI2->V[MATI2->R - 1][i]); 
		printf(" ");
	}
	printf("\n");
	MATI3 = BUILDMATI(1, m);
	for (i = 0;i < m;i++)
		MATI3->V[0][i] = COPYI(MATI2->V[m - 1][i]); 
	T = LENGTHSQRI(MATI3, 0);
	printf("The Euclidean norm squared = ");
	PRINTI(T);
	FREEMPI(T);
	printf("\n");
	printf("Do you want to get a shortest multiplier using Fincke_Pohst? (Y/N)\n");
	answer = GetYN();
	p = m - 1;
	XX = (MPI **)mmalloc((USL)(p * sizeof(MPI *)));
	for (j = 0; j < p; j++)
		XX[j] = ZEROI();
	if (answer)
	{
		GCDVERBOSE = 1;
		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;
			}
		}
		printf("found a shortest multiplier vector:\n");
	}
	else
		Q = MATI2;
	strcpy(buff, "lllgcdmat.out");
	outfile = fopen(buff, "w");
	FPRINTMATI(outfile,0,Q->R-1,0,Q->C-1,Q);
	fclose(outfile);
	strcpy(buff, "lllgcdmult.out");
	outfile = fopen(buff, "w");
	if (answer)
		fprintf(outfile, "A Shortest multiplier is ");
	else
	{
		fprintf(outfile, "A not necessarily shortest multiplier is ");
		fprintf(outfile, "b[%u]=", m);
	}
	if (answer)
	{
		printf("b[%u]", m);
		fprintf(outfile, "b[%u]", m);
		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)
			{