i6r.c

calc大数库

/* i6R.c */
/* matrix operations using MPR's */

#include <stdio.h>
/*#include <unistd.h>*/
#include <string.h>
#include <stdlib.h>
#include "integer.h"
#include "fun.h"
MPMATR *TRANSPOSER(MPMATR *Mptr);
MPMATR *MULTMATR(MPMATR *Mptr, MPMATR *Nptr);
extern USI GCDVERBOSE;
extern USI GCD3FLAG;
extern USI GCD_MAX;

MPMATR *BUILDMATR(USI m, USI n)
/*
 * Allocates space for an m x n matrix of MPR's.
 */
{
	MPMATR *N;

	N=(MPMATR *)mmalloc(sizeof(MPMATR));
	N->R = m;
	N->C = n;
	N->V = (MPR **)mmalloc(m * n * sizeof(MPR *));
	return (N);
}

MPMATR *COPYMATR(MPMATR *Mptr)
/*
 * a replacement for the assignment *Nptr = *Mptr.
 */
{
	unsigned int i, t;
	MPMATR *Nptr;

	t = Mptr->R * Mptr->C;
	Nptr = BUILDMATR(Mptr->R, Mptr->C);
	for (i = 0; i < t; i++)
			Nptr->V[i] = COPYR(Mptr->V[i]);
	return (Nptr);
}

void FREEMATR(MPMATR *Mptr)
/*
 * frees the storage allocated to the two-dimensional array Mptr->V.
 */
{
	unsigned int i, t;

	if (Mptr == NULL)
		return;
	t = Mptr->R * Mptr->C;
	for (i = 0; i < t; i++)
		FREEMPR(Mptr->V[i]);
	ffree((char *)(Mptr->V), t * sizeof(MPR *));
	ffree((char *)Mptr, sizeof(MPMATR));
	return;
}

MPMATR *CHOLESKYR(MPMATR *A)
/*
 * Input: A positive definite matrix A.
 * Output: The Cholesky decomposition of A.
 * See U.Finke and M. Pohst, "Improved methods for calculating vectors of
 * short length in a lattice, including a complexity analysis." Math. Comp.
 * 44, 463-471, 1985
 */
{
	MPMATR *Q;
	unsigned int i, j, k, l, m;
	MPR *temp1, *temp2;

	m = A->R;
	Q = ZEROMNR(m, m);
	for (i = 0; i < m; i++)
	{
		for (j = i; j < m; j++)
		{
			FREEMPR(elt(Q, i, j));
			elt(Q, i, j) = COPYR(elt(A, i, j));
		}
	}
	for (i = 0; i < m - 1; i++)
	{
		for (j = i + 1; j < m; j++)
		{
			FREEMPR(elt(Q, j, i));
			elt(Q, j, i) = COPYR(elt(Q, i, j));
			temp1 = elt(Q, i, j);
			elt(Q, i, j) = RATIOR(elt(Q, i, j), elt(Q, i, i));
			FREEMPR(temp1);
		}
		for (k = i + 1; k < m; k++)
		{
			for (l = k; l < m; l++)
			{
				temp1 = elt(Q, k, l);
				temp2 = MULTR(elt(Q, k, i), elt(Q, i, l));
				elt(Q, k, l) = SUBR(temp1, temp2);
				FREEMPR(temp1);
				FREEMPR(temp2);
			}
		}
	}
	for (j = 0; j < m - 1; j++)
	{
		for (i = j + 1; i < m; i++)
		{
			FREEMPR(elt(Q, i, j));
			elt(Q, i, j) = ZEROR();
		}
	}
	return (Q);
}

void FINCKE_POHST(MPMATR *A, MPR *C, USI filename)
/*
 * Input: A matrix of integers A with LI rows spanning a lattice L.
 * Output:  The integer vectors X with ||X||^2 <= C, highest nonzero coord <=0.
 * See "Improved methods for calculating vectors of short length in a lattice,
 * including a complexity analysis, U. Fincke and M. Pohst,
 * Mathematics of Computation, 44, 1985, 463-471.
 */
{
	unsigned int i, j, n, flag, *l, *u, *t, *x, count = 0;
	unsigned int FLAG;
	MPR **X, **L, **T, **U,  *Temp1, *Temp2, *Temp3, *Z, *ONE, *SUM;
	MPR *temp;
	char buff[20]; 
	FILE *outfile;
	MPMATR *B, *BB, *Q, *VECTOR, *MATR;
	MPMATI *MATI;
	int e;
	
	B = TRANSPOSER(A);
	BB = MULTMATR(A, B);
	FREEMATR(B);
	Q = CHOLESKYR(BB);
	FREEMATR(BB);
	if (filename == 1)
		strcpy(buff, "fp.out");
	else
		strcpy(buff, "slv.out");
	outfile = fopen(buff, "w");
	n = Q->R;
	i = n - 1;
	ONE = ONER();
	X = (MPR **)mmalloc(n * sizeof(MPR *));
	L = (MPR **)mmalloc(n * sizeof(MPR *));
	T = (MPR **)mmalloc(n * sizeof(MPR *));
	U = (MPR **)mmalloc(n * sizeof(MPR *));
	l = (USI *)mmalloc(n * sizeof(USI));
	t = (USI *)mmalloc(n * sizeof(USI));
	u = (USI *)mmalloc(n * sizeof(USI));
	x = (USI *)mmalloc(n * sizeof(USI));
	for (j = 0; j < n; j++)
	{
		l[j] = 0;
		t[j] = 0;
		u[j] = 0;
		x[j] = 0;
	}
	T[i] = COPYR(C); t[i] = 1;
	U[i] = ZEROR(); u[i] = 1;
bounds:
	Z = RATIOR(T[i], elt(Q, i, i));
	Temp2 = MINUSR(U[i]);
	if (l[i])
		FREEMPR(L[i]); 
	else
		l[i] = 1;
	L[i] = INTROOT(Z, Temp2);
	FREEMPR(Temp2);
	Temp2 = INTROOT(Z, U[i]); 
	Temp3 = MINUSR(Temp2);
	if (x[i])
		FREEMPR(X[i]);
	else
		x[i] = 1;
	X[i] = SUBR(Temp3, ONE);
	FREEMPR(Z);
	FREEMPR(Temp2);
	FREEMPR(Temp3);
loop:
	Temp1 = X[i];
	X[i] = ADDR(X[i], ONE);
	FREEMPR(Temp1);
	if (COMPAREI(X[i]->N, L[i]->N) == 1)
	{
		i++;
		goto loop;
	}
	else
	{
		if (i)
		{
			Temp1 = ADDR(X[i], U[i]);
			Temp2 = MULTR(Temp1, Temp1);
			FREEMPR(Temp1);
			Temp1 = MULTR(elt(Q, i, i), Temp2);
			FREEMPR(Temp2);
			if (t[i - 1])
				FREEMPR(T[i - 1]);
			else
				t[i - 1] = 1;
			T[i - 1] = SUBR(T[i], Temp1);
			FREEMPR(Temp1);
			i--;
			SUM = ZEROR();
			for (j = i + 1; j < n; j++)
			{
				Temp1 = SUM;
				Temp2 = MULTR(elt(Q, i, j), X[j]);
				SUM = ADDR(SUM, Temp2);
				FREEMPR(Temp1);
				FREEMPR(Temp2);
			}
			if (u[i])
				FREEMPR(U[i]);
			else
				u[i] = 1;
			U[i] = COPYR(SUM);
			FREEMPR(SUM);
			goto bounds;
		}
		else
			goto found;
	}
	
	found:
		flag = 0;
		for (j = 0; j < n; j++)
		{
			if (!EQZEROR(X[j]))
			{
				flag = 1;
				break;
			}
		}
		if (flag == 0)
		{
			FREEMPR(ONE);
			for (j = 0; j < n; j++)
			{
				FREEMPR(X[j]);
				FREEMPR(L[j]);
				FREEMPR(T[j]);
				FREEMPR(U[j]);
			}
			ffree((char *)X, n * sizeof(MPR *));
			ffree((char *)L, n * sizeof(MPR *));
			ffree((char *)T, n * sizeof(MPR *));
			ffree((char *)U, n * sizeof(MPR *));
			ffree((char *)l, n * sizeof(USI));
			ffree((char *)t, n * sizeof(USI));
			ffree((char *)u, n * sizeof(USI));
			ffree((char *)x, n * sizeof(USI));
			fclose(outfile);
			FREEMATR(Q);
			return;
		}
		else
		{
			VECTOR = BUILDMATR(1, A->R);
		printf("LATTICE_VECTOR[%u] = ", count + 1);
		fprintf(outfile, "LATTICE_VECTOR[%u]:", count + 1);
		FLAG = 0;
		for (j = 0; j < A->R; j++)
		{
			e = (X[j]->N)->S;
			if (e == -1)
			{
				if (FLAG)
				{
					printf("+");
					fprintf(outfile, "+");
				}
				if (FLAG == 0)
					FLAG = 1;
				temp = MINUSR(X[j]);
				if (!EQONER(temp))
				{
					PRINTR(temp);
					FPRINTR(outfile, temp);
				}
				printf("b[%u]", j + 1);
				fprintf(outfile, "b[%u]", j + 1);
				FREEMPR(temp);
			}
			if (e == 1)
			{
				if (FLAG == 0)
					FLAG = 1;
				printf("-");
				fprintf(outfile, "-");
				if (!EQONER(X[j]))
				{
					PRINTR(X[j]);
					FPRINTR(outfile, X[j]);
				}
				printf("b[%u]", j + 1);
				fprintf(outfile, "b[%u]", j + 1);
			}
			elt(VECTOR, 0, j) = MINUSR(X[j]);
		}
		MATR = MULTMATR(VECTOR, A);
		FREEMATR(VECTOR);
		printf("=");
		fprintf(outfile, "=");
		MATI = BUILDMATI(1, A->C);
		for (j = 0; j < A->C; j++)
		{
			MATI->V[0][j] = COPYI((elt(MATR, 0, j))->N);
			PRINTI(MATI->V[0][j]);
			FPRINTI(outfile, MATI->V[0][j]);
			printf(" ");
			fprintf(outfile, " ");
		}
		printf(": ");
		fprintf(outfile, ": ");
		FREEMATR(MATR);
		FREEMATI(MATI);
		Temp1 = ADDR(X[0], U[0]);
			Temp2 = MULTR(Temp1, Temp1);
			FREEMPR(Temp1);
			Temp1 = MULTR(elt(Q, 0, 0), Temp2);
			FREEMPR(Temp2);
			Temp2 = SUBR(C, T[0]);
			Temp3 = ADDR(Temp2, Temp1);
			FREEMPR(Temp1);
			FREEMPR(Temp2);
			PRINTR(Temp3); 
			FPRINTR(outfile, Temp3); 
			printf("\n");
			fprintf(outfile, "\n");
			FREEMPR(Temp3);
			count++;
			goto loop;
		}
}

MPR *SLVECTOR(MPMATR *A, MPR *C, MPR **VALUE)
/*
 * Input: A matrix of integers A with LI LLL reduced rows spanning a lattice L.
 * C is the length-squared of a small lattive vector.
 * Output: if 0 is returned, this means that VALUE is the shortest length.
 * Then in nfunc.c, FINCKE-POHST is applied to get all the shortest vectors in
 * L with highest nonzero coord > 0 are printed.
 * Otherwise a shorter length is returned and VALUE will be NULL.
 * See "Improved methods for calculating vectors of short length in a lattice,
 * including a complexity analysis, U. Fincke and M. Pohst,
 * Mathematics of Computation, 44, 1985, 463-471.
 */
{
	unsigned int i, j, n, flag, *l, *u, *t, *x, count = 0;
	MPR **X, **L, **T, **U,  *Temp1, *Temp2, *Temp3, *Z, *ONE, *SUM;
	MPMATR *B, *BB, *Q;
	int s1, s2;
	
	*VALUE = (MPR *)NULL;
	B = TRANSPOSER(A);
	BB = MULTMATR(A, B);
	FREEMATR(B);
	Q = CHOLESKYR(BB);
	FREEMATR(BB);
	n = Q->R;
	i = n - 1;
	ONE = ONER();
	X = (MPR **)mmalloc(n * sizeof(MPR *));
	L = (MPR **)mmalloc(n * sizeof(MPR *));
	T = (MPR **)mmalloc(n * sizeof(MPR *));
	U = (MPR **)mmalloc(n * sizeof(MPR *));
	l = (USI *)mmalloc(n * sizeof(USI));
	t = (USI *)mmalloc(n * sizeof(USI));
	u = (USI *)mmalloc(n * sizeof(USI));
	x = (USI *)mmalloc(n * sizeof(USI));
	for (j = 0; j < n; j++)
	{
		l[j] = 0;
		t[j] = 0;
		u[j] = 0;
		x[j] = 0;
	}
	T[i] = COPYR(C); t[i] = 1;
	U[i] = ZEROR(); u[i] = 1;
bounds:
	Z = RATIOR(T[i], elt(Q, i, i));
	Temp2 = MINUSR(U[i]);
	if (l[i])
		FREEMPR(L[i]); 
	else
		l[i] = 1;
	L[i] = INTROOT(Z, Temp2);
	FREEMPR(Temp2);
	Temp2 = INTROOT(Z, U[i]); 
	Temp3 = MINUSR(Temp2);
	if (x[i])
		FREEMPR(X[i]);
	else
		x[i] = 1;
	X[i] = SUBR(Temp3, ONE);
	FREEMPR(Z);
	FREEMPR(Temp2);
	FREEMPR(Temp3);
loop:
	Temp1 = X[i];
	X[i] = ADDR(X[i], ONE);
	FREEMPR(Temp1);
	if (COMPAREI(X[i]->N, L[i]->N) == 1)
	{
		i++;
		goto loop;
	}
	else
	{
		if (i)
		{
			Temp1 = ADDR(X[i], U[i]);
			Temp2 = MULTR(Temp1, Temp1);
			FREEMPR(Temp1);
			Temp1 = MULTR(elt(Q, i, i), Temp2);
			FREEMPR(Temp2);
			if (t[i - 1])
				FREEMPR(T[i - 1]);
			else
				t[i - 1] = 1;
			T[i - 1] = SUBR(T[i], Temp1);
			FREEMPR(Temp1);
			i--;
			SUM = ZEROR();
			for (j = i + 1; j < n; j++)
			{
				Temp1 = SUM;
				Temp2 = MULTR(elt(Q, i, j), X[j]);
				SUM = ADDR(SUM, Temp2);
				FREEMPR(Temp1);
				FREEMPR(Temp2);
			}
			if (u[i])
				FREEMPR(U[i]);
			else
				u[i] = 1;
			U[i] = COPYR(SUM);
			FREEMPR(SUM);
			goto bounds;
		}
		else
			goto found;
	}
	
	found:
		flag = 0;
		for (j = 0; j < n; j++)
		{
			if (!EQZEROR(X[j]))
			{
				flag = 1;
				break;
			}
		}
		if (flag == 0)
		{
			FREEMPR(ONE);
			for (j = 0; j < n; j++)
			{
				FREEMPR(X[j]);
				FREEMPR(L[j]);
				FREEMPR(T[j]);
				FREEMPR(U[j]);
			}
			ffree((char *)X, n * sizeof(MPR *));
			ffree((char *)L, n * sizeof(MPR *));
			ffree((char *)T, n * sizeof(MPR *));
			ffree((char *)U, n * sizeof(MPR *));
			ffree((char *)l, n * sizeof(USI));
			ffree((char *)t, n * sizeof(USI));
			ffree((char *)u, n * sizeof(USI));
			ffree((char *)x, n * sizeof(USI));
			FREEMATR(Q);
			return (ZEROR());
		}
		else
		{
			Temp1 = ADDR(X[0], U[0]);
			Temp2 = MULTR(Temp1, Temp1);
			FREEMPR(Temp1);
			Temp1 = MULTR(elt(Q, 0, 0), Temp2);
			FREEMPR(Temp2);
			Temp2 = SUBR(C, T[0]);
			Temp3 = ADDR(Temp2, Temp1);
			FREEMPR(Temp1);
			FREEMPR(Temp2);
			s1 = (Temp3->N)->S;
			if (s1)
			{
				if (*VALUE != NULL)
					FREEMPR(*VALUE);
				*VALUE = COPYR(Temp3);
			}
			s2 = RSV(Temp3->N, C->N);
			if (s1 && s2 == -1)
			{
				FREEMPR(ONE);
				for (j = 0; j < n; j++)
				{
					FREEMPR(X[j]);
					FREEMPR(L[j]);
					FREEMPR(T[j]);
					FREEMPR(U[j]);
				}
				ffree((char *)X, n * sizeof(MPR *));
				ffree((char *)L, n * sizeof(MPR *));
				ffree((char *)T, n * sizeof(MPR *));
				ffree((char *)U, n * sizeof(MPR *));
				ffree((char *)l, n * sizeof(USI));
				ffree((char *)t, n * sizeof(USI));
				ffree((char *)u, n * sizeof(USI));
				ffree((char *)x, n * sizeof(USI));
				FREEMATR(Q);
				return(Temp3);
			}
			FREEMPR(Temp3);
			count++;
			goto loop;
		}
}

MPR *INTROOT(MPR *Z, MPR *U)
/*
 * Returns [sqrt(Z)+U]. First ANSWER = [sqrt(Z)] + [U]. One then
 * tests if Z < ([sqrt(Z)] + 1 -{U})^2. If this does not hold, ANSWER++.
 * For use in FINCKE_POHST().
 */
{
	int t;
	MPR *X, *Y, *ANSWER, *R, *S, *T1, *T2, *ONE, *Temp;

	if (EQZEROR(Z))
		return (INTR(U));
	ONE = ONER();
	X = MTHROOTR(Z, 2, 0);
	Y = INTR(U);
	ANSWER = ADDR(X, Y);
	R = SUBR(U, Y);
	FREEMPR(Y);
	S = SUBR(ONE, R);
	FREEMPR(R);
	T1 = ADDR(X, S);
	FREEMPR(S);
	FREEMPR(X);
	T2 = MULTR(T1, T1);
	FREEMPR(T1);
	t = COMPARER(T2, Z);
	FREEMPR(T2);
	if (t <= 0)
	{
		Temp = ANSWER;
		ANSWER = ADDR(ANSWER, ONE);
		FREEMPR(Temp);
	}
	FREEMPR(ONE);
	return (ANSWER);
}

int COMPARER(MPR *Aptr, MPR *Bptr)
/*
 *				    1 if *Aptr > *Bptr
 *            Returns               0 if *Aptr = *Bptr
 *       			   -1 if *Aptr < *Bptr
 */
{
	MPI *T1, *T2;
	int t;

	T1 = MULTI(Aptr->N, Bptr->D);
	T2 = MULTI(Aptr->D, Bptr->N);
	t = COMPAREI(T1, T2);
	FREEMPI(T1);
	FREEMPI(T2);
	return (t);
}

MPMATR *ZEROMNR(USI m, USI n)
/*
 * returns the zero  m x n matrix.
 */
{
	unsigned int i, j;
	MPMATR *Mptr;

	Mptr = BUILDMATR(m, n);
	for (i = 0; i <= m - 1; i++)
	{
		for (j = 0; j <= n - 1; j++)
			elt(Mptr, i, j) = ZEROR();
	}
	return (Mptr);
}