www.pudn.com > calc.rar > lllhermi.c

/* lllhermi.c */ 
#include  
#include  
#include "integer.h" 
#include "fun.h" 
extern MPI *MAXI, *PMAXI; 
 
extern unsigned int MLLLVERBOSE; 
extern unsigned int HERMITE1VERBOSE; 
unsigned int keith97; 
 
USI FLAGCOLI(MPMATI *Aptr); 
 
 
MPMATI *LLLHERMITE1(MPMATI *DD, MPMATI **Aptr, USI *rank, USI m1, USI n1) 
/*  
 * Input: an m x n matrix DD of MPI's. 
 * Output: *Aptr = HNF(DD). Also we return a small transformation matrix B 
 * if nullity > 0. 
 * The limiting form of  LLL applied to [I_m|N^mD_1|N^{m-1}D_2|...] is used. 
 * Normally m1/n1 = 3/4. 
 * 23/2/97. 
 */ 
{ 
	unsigned int *COL; 
	unsigned int i, j, k, l, m, n, col1, col2, r; 
	int t; 
	MPI **D, *X, *Y, *Z, *Tmp, *R, ***temp1, ***temp2; 
	MPI *M1, *N1; 
	MPMATI *L, *B; 
 
	keith97 = FLAGCOLI(DD); 
	m = DD->R; 
	n = DD->C; 
	B = IDENTITYI(m); 
	*Aptr = COPYMATI(DD); 
	if (!keith97){ 
		(B->V[m - 1][m - 1])->S = -1; 
		for (j = 0; j < n; j++) 
			((*Aptr)->V[m - 1][j])->S = - (((*Aptr)->V[m - 1][j])->S); 
	} 
	D = (MPI **)mmalloc((1 + m) * sizeof(MPI *)); 
	for (i = 0; i <= m; i++) 
		D[i] = ONEI(); 
 
	L = ZEROMNI(m, m); 
 
	k = 2; 
	M1 = CHANGE(m1); 
	N1 = CHANGE(n1); 
	while (k <= m) 
	{ 
		REDUCE1(k, k - 1, &L, &B, D, Aptr, &col1, &col2); 
		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 ((RSV(Y, Z) == 1 && col1 == n && col2 == n ) || (col1 <= col2 && col1 < n )) 
		{ 
			for (i = k + 1; i <= m; i++) 
				SWAP1(i, k, &L, D); 
			SWAP21(k, &B, &L, Aptr); 
			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--) 
				REDUCE1(k, l, &L, &B, D, Aptr, &col1, &col2); 
			k++; 
		} 
		FREEMPI(Z);		 
	} 
	FREEMPI(M1); 
	FREEMPI(N1); 
	FREEMATI(L); 
	for (i = 0; i <= m; i++) 
		FREEMPI(D[i]); 
 
	COL = (USI *)mmalloc((1 + m) * sizeof(USI)); 
	COL[0] = 0; 
	for (t = m - 1; t >= 0; t--) 
	{ 
		r = COLSEEKI0((*Aptr), t, COL[m-1-t]); 
		if (r == n)  
			break; 
		else 
			COL[m-t] = r; 
	} 
	if (t == -1)  
		*rank = m; 
	else 
		*rank = m - t - 1; 
 
	temp1 = (MPI ***)mmalloc((*rank) * sizeof(MPI **)); 
	temp2 = (MPI ***)mmalloc((*rank) * sizeof(MPI **)); 
	for (i = 0; i < *rank; i++) 
	{ 
		temp1[i] = (*Aptr)->V[m - 1 - i]; 
		temp2[i] = B->V[m - 1 - i]; 
	} 
	for (i = m - 1; i >= *rank; i--) 
	{ 
		(*Aptr)->V[i] = (*Aptr)->V[i - *rank]; 
		B->V[i] = B->V[i - *rank]; 
	} 
	for (i = 0; i < *rank; i++) 
	{ 
		(*Aptr)->V[i] = temp1[i]; 
		B->V[i] = temp2[i]; 
	} 
 
	ffree((char *)COL, (1 + m) * sizeof(USI *)); 
	ffree((char *)temp1, (*rank) * sizeof(MPI **)); 
	ffree((char *)temp2, (*rank) * sizeof(MPI **)); 
	ffree((char *)D, (1 + m) * sizeof(MPI *)); 
	return (B); 
} 
 
 
void REDUCE1(USI k, USI l, MPMATI **Lptr, MPMATI **Bptr, MPI *D[], MPMATI **Aptr, USI *col1, USI *col2) 
/* 
 * updates *Lptr, *Bptr, *Aptr. 
 */ 
{ 
	unsigned int j, n, t; 
	MPI *X, *Y, *Q, *Tmp; 
	MPMATI *TempMATI; 
 
	n = (*Aptr)->C; 
	*col1 =  n; 
	for (j = 0; j < n; j++) 
	{ 
		t = ((*Aptr)->V[l - 1][j])->S; 
		if (t) 
		{ 
			*col1 = j; 
			if (t == -1) 
				LLLMINUS(Aptr, Bptr, Lptr, l - 1); 
			break; 
		} 
	} 
	*col2 =  n; 
	for (j = 0; j < n; j++) 
	{ 
		t = ((*Aptr)->V[k - 1][j])->S; 
		if (t) 
		{ 
			*col2 = j; 
/* 
			if ((k == (*Aptr)->R) && !keith97) 
			{ 
				if (t == -1) 
					LLLMINUS(Aptr, Bptr, Lptr, k - 1); 
			} 
*/ 
			break; 
		} 
	} 
 
	if (*col1 < n) 
/*		Q = NEAREST_INTI((*Aptr)->V[k - 1][*col1], (*Aptr)->V[l - 1][*col1]);*/ 
		Q = INTI((*Aptr)->V[k - 1][*col1], (*Aptr)->V[l - 1][*col1]); 
	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); 
	TempMATI = *Bptr; 
	*Bptr = ADD_MULT_ROWI(l - 1, k - 1, X, *Bptr); 
	FREEMATI(TempMATI); 
	TempMATI = *Aptr; 
	*Aptr = ADD_MULT_ROWI(l - 1, k - 1, X, *Aptr); 
	FREEMATI(TempMATI); 
	if (HERMITE1VERBOSE) 
	{ 
		printf("Row %u -> Row %u + ", k,k);PRINTI(X);printf(" x Row %u\n", l); 
		printf("P = \n"); 
		PRINTMATI(0,(*Bptr)->R-1,0,(*Bptr)->C-1,*Bptr); 
		printf("A = \n"); 
		PRINTMATI(0,(*Aptr)->R-1,0,(*Aptr)->C-1,*Aptr); 
		GetReturn(); 
	} 
	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); 
	return; 
} 
 
 
void SWAP21(USI k, MPMATI **Bptr, MPMATI **Lptr, MPMATI **Aptr) 
{ 
	MPI *T; 
 
	*Aptr = SWAP_ROWSI1(k - 2, k - 1, *Aptr); 
	*Bptr = SWAP_ROWSI1(k - 2, k - 1, *Bptr); 
	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 (HERMITE1VERBOSE) 
	{ 
		printf("Swapping Rows %u and %u\n", k-1,k); 
		printf("P = \n"); 
		PRINTMATI(0,(*Bptr)->R-1,0,(*Bptr)->C-1,*Bptr); 
		printf("A = \n"); 
		PRINTMATI(0,(*Aptr)->R-1,0,(*Aptr)->C-1,*Aptr); 
		GetReturn(); 
		return; 
	} 
} 
 
void LLLMINUS(MPMATI **Aptr, MPMATI **Pptr, MPMATI **Lptr, USI j) 
{ 
	USI r, s, m; 
	MPI *Temp; 
	MPMATI *TempMATI; 
 
	m = (*Aptr)->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 = *Aptr; 
	Temp = MINUS_ONEI(); 
	*Aptr = SCALAR_ROWI(j, Temp, *Aptr); 
	FREEMPI(Temp); 
	FREEMATI(TempMATI); 
	TempMATI = *Pptr; 
	Temp = MINUS_ONEI(); 
	*Pptr = SCALAR_ROWI(j, Temp, *Pptr); 
	FREEMPI(Temp); 
	FREEMATI(TempMATI); 
	if (HERMITE1VERBOSE) 
	{ 
		printf("Row %u -> - Row %u\n", j + 1,j + 1); 
		printf("P = \n"); 
		PRINTMATI(0,(*Pptr)->R-1,0,(*Pptr)->C-1,*Pptr); 
		printf("A = \n"); 
		PRINTMATI(0,(*Aptr)->R-1,0,(*Aptr)->C-1,*Aptr); 
		GetReturn(); 
	} 
	return; 
} 
 
MPMATI *SCALAR_ROWI(USI p, MPI *Aptr, MPMATI *Mptr) 
/* 
 * multiplying the pth row of *Mptr by *Aptr. 
 */ 
{ 
	unsigned int j; 
	MPI *Temp; 
	MPMATI *Nptr; 
 
	Nptr = COPYMATI(Mptr); 
	for (j = 0; j <= Nptr->C - 1; j++) 
	{ 
		Temp = Nptr->V[p][j]; 
		Nptr->V[p][j] = MULTI(Temp, Aptr); 
		FREEMPI(Temp); 
	} 
	return (Nptr); 
} 
 
USI COLSEEKI0(MPMATI *A, USI i, USI j) 
/* Scans cols >= i, rows >= j for first nonzero column if any. Returns  
 * this column number if ther is one, else returns A->C; 
 */ 
{ 
	USI k, l, n; 
 
	n = A->C; 
	for (k = j; k < n; k++) 
	{ 
		for (l = 0; l <= i; l++) 
		{ 
			if (A->V[l][k]->S) 
				return k; 
		} 
	} 
	return k; 
} 
 
USI FLAGCOLI(MPMATI *Aptr) 
/* 
 * Returns 1 if the first nonzero column of Aptr contains more than one 
 * entry, 0 otherwise.  Only applied to a nonzero matrix with at least  
 * two rows. 
 */ 
{ 
	USI i, j, k, flag = 0, m, n; 
 
	m = Aptr->R; 
	n = Aptr->C; 
	for (j = 0; j < n; j++) 
	{ 
		for (i = 0; i < m; i++) 
		{ 
			if ((Aptr->V[i][j])->S) 
			{ 
				flag = 1; 
				break; 
			} 
		} 
		if (flag) 
		{ 
			for (k = i + 1; k < m; k++) 
			{ 
				if ((Aptr->V[k][j])->S) 
					return (1); 
			} 
			return (0); 
		} 
	} 
	return (0); 
}