i6i.c

calc大数库

/*
 * sorts the array T[0],...,T[nz-1] into ascending order, thereby returning
 * a permutation of 0,...,nz-1 used in PERMUTE_MATI below.
 */
{
	unsigned int i, j, temp, *P;

	P = (unsigned int *) mmalloc(nz * sizeof(unsigned int));
	for (i = 0; i < nz ; i++)
		P[i] = i;
	for (i = 0; i < nz - 1; i++)
	{
		for (j = i + 1; j < nz; j++)
			if (T[i] > T[j])
			{
				temp = T[i];
				T[i] = T[j];
				T[j] = temp;
				temp = P[i];
				P[i] = P[j];
				P[j] = temp;
			}
	}
	for (i = 0; i < nz; i++)
		T[i] = P[i];
	ffree ((char *) P, nz * sizeof(unsigned int));
	return;
}

MPMATI *HERMITE1(MPMATI *Aptr, USI *T, USI nz)
/*
 * Takes the output of KB_ROW() and permutes the rows to get the Hermite normal
 * form from Aptr.
 */
{
	MPMATI *Tmp;
	unsigned int i, j;
	
	SORTI(T, nz);
	Tmp = COPYMATI(Aptr);
	for (i = 0; i < nz; i++)
	{
		for (j = 0; j < Aptr->C; j++)
		{
			FREEMPI(Tmp->V[i][j]);
			Tmp->V[i][j] = COPYI(Aptr->V[T[i]][j]);
		}
	}
	return (Tmp);	
}

MPMATI *ROWSUBI(USI p, USI q, MPI *Aptr, MPMATI *Mptr)
/*
 * Updates Mptr: subtracts *Aptr times the p-th row of *Mptr from the q-th.
 * 0 <= p, q <= Mprt->R - 1.
 */
{
	MPI *Y, *Temp;
	unsigned int l;
	if (Aptr->S != 0)
	{
		for (l = 0; l <= Mptr->C - 1; l++) 
		{
			
			if ((Mptr->V[p][l])->S != 0)
			{
				Y = MULTI(Aptr, Mptr->V[p][l]);
				Temp = Mptr->V[q][l];
				Mptr->V[q][l] = SUBI(Mptr->V[q][l], Y);
				FREEMPI(Y);
				FREEMPI(Temp);
			}
		}
	}
	return (Mptr);
}


MPMATI *COLSUBI(USI p, USI q, MPI *Aptr, MPMATI *Mptr)
/*
 * subtracts *Aptr times the p-th column of *Mptr from the q-th.
 * 0 <= p, q <= Mprt->C - 1.
 */
{
	MPI *Y, *Temp;
	unsigned int k;

	if (Aptr->S != 0)
	{
		for (k = 0; k <= Mptr->R - 1; k++) 
		{
			if ((Mptr->V[k][p])->S != 0)
			{
				Y = MULTI(Aptr, Mptr->V[k][p]);
				Temp = Mptr->V[k][q];
				Mptr->V[k][q] = SUBI(Mptr->V[k][q], Y);
				FREEMPI(Y);
				FREEMPI(Temp);
			}
		}
	}
	return (Mptr);
}

MPMATI *IDENTITYI(USI n)
/*
 * returns the identity matrix of size n.
 */
{
	unsigned int i, j;
	MPMATI *Mptr;

	Mptr = BUILDMATI(n, n);
	for (i = 0; i <= n - 1; i++)
	{
		for (j = 0; j <= n - 1; j++)
		{
			if ( i == j)
				Mptr->V[i][j] = ONEI();
			else
				Mptr->V[i][j] = ZEROI();
		}
	}
	return (Mptr);
}

MPMATI *MULTMATI(MPMATI *Mptr, MPMATI *Nptr)
/*
 * Here Mptr->C = Nptr->R. 
 * returns (*Mptr) * (*Nptr).
 */
{
	MPI *X, *Y, *Temp;
	unsigned int i, j, k;
	MPMATI *Lptr;

	Lptr = BUILDMATI(Mptr->R, Nptr->C);
	for (i = 0; i <= Mptr->R - 1; i++)
	{
		for (j = 0; j <= Nptr->C - 1; j++)
		{
			X = ZEROI();
			for (k = 0; k <= Mptr->C - 1; k++)
			{
				Y = MULTI(Mptr->V[i][k], Nptr->V[k][j]);
				Temp = X;
				X = ADDI(X, Y);
				FREEMPI(Temp);
				FREEMPI(Y);
			}
			Lptr->V[i][j] = X;
		}
	}
	return (Lptr);
}

MPMATI *ADD_MULT_ROWI(USI p, USI q, MPI *Aptr, MPMATI *Mptr)
/*
 * adding *Aptr times row p to row q of *Mptr.
 */
{
	unsigned int j;
	MPI *X, *Temp;
	MPMATI *Nptr;

	Nptr = COPYMATI(Mptr);
	for (j = 0; j <= Nptr->C - 1; j++)
	{
		X = MULTI(Nptr->V[p][j], Aptr);
		Temp = Nptr->V[q][j];
		Nptr->V[q][j] = ADDI(Temp, X);
		FREEMPI(X);
		FREEMPI(Temp);
	}
	return (Nptr);
}

void *ADD_MULT_ROWI0(USI p, USI q, MPI *Aptr, MPMATI *Mptr)
/*
 * adding *Aptr times row p to row q of *Mptr and overwrites *Mptr.
 */
{
	unsigned int j;
	MPI *X, *Temp;

	for (j = 0; j <= Mptr->C - 1; j++)
	{
		X = MULTI(Mptr->V[p][j], Aptr);
		Temp = Mptr->V[q][j];
		Mptr->V[q][j] = ADDI(Temp, X);
		FREEMPI(X);
		FREEMPI(Temp);
	}
	return (Mptr);
}

MPMATI *TRANSPOSEI(MPMATI *Mptr)
/* 
 * returns the transpose of *Mptr.
 */
{
	MPMATI *Nptr;
	unsigned int i, j;

	Nptr = BUILDMATI(Mptr->C, Mptr->R);
	for (j = 0; j < Nptr->R; j++)
		for (i = 0; i < Nptr->C; i++)
			Nptr->V[j][i]= COPYI(Mptr->V[i][j]);
	return (Nptr);
}

MPMATI *DELETE_ROWI(USI r, MPMATI *Mptr)
/*
 * deletes row r of *Mptr.
 */
{
	unsigned int i, j;
	MPMATI *Nptr;

	r--;
	Nptr = BUILDMATI(Mptr->R - 1, Mptr->C);
	for (i = 0; i <= Nptr->R - 1; i++)
	{
		if (i < r)
		{
			for (j = 0; j <= Nptr->C - 1; j++)
				Nptr->V[i][j] = COPYI(Mptr->V[i][j]);
		}
		else
		{
			for (j = 0; j <= Nptr->C - 1; j++)
				Nptr->V[i][j] = COPYI(Mptr->V[i + 1][j]);
		}
	}
	return (Nptr);
}

unsigned int *KB_ROWP(MPMATI *Aptr, MPMATI **Pptr, USI *nz)
/*
 * The Kannan-Bachem reduction to column permuted Hermite normal form, for
 * use in their Smith Normal form algorithm. The transforming unimodular matrix
 * **Pptr is returned.
 */
{
	unsigned int *T, i, j, k, s, m, n, tmp, is_zero;
	MPI *D, *P, *Q, *R, *S, *Tmp0, *Tmp1, *Tmp2 , *Tmp3;
        MPI *Tmp0P, *Tmp1P, *Tmp2P, *Tmp3P;

	m = Aptr->R;
	n = Aptr->C;
	*Pptr = IDENTITYI(m);
	T = (unsigned int *)mmalloc(n * sizeof(unsigned int));
	for (i = 0; i < n; i++)
		T[i] = i;
	s = m - 1;
	/* s = index of last row of A not known to be zero. */
	i = 0;
	while (i <= s)
	{
if(MLLLVERBOSE)
		printf("i = %u\n", i);
		for (k = 0; k < i; k++)
		{
			D = EUCLIDI(Aptr->V[k][T[k]], Aptr->V[i][T[k]], &P, &Q);
			Tmp0 = INT(Aptr->V[i][T[k]], D);
			R = MINUSI(Tmp0);
			FREEMPI(Tmp0);
			S = INT(Aptr->V[k][T[k]], D);
			FREEMPI(D);
			for (j = 0; j < n; j++)
			{
				Tmp0 = Aptr->V[k][j];
				Tmp1 = MULTI(P, Aptr->V[k][j]);
				Tmp2 = MULTI(Q, Aptr->V[i][j]);
				Aptr->V[k][j] = ADDI(Tmp1, Tmp2);
				FREEMPI(Tmp1);
				FREEMPI(Tmp2);
				Tmp3 = Aptr->V[i][j];
				Tmp1 = MULTI(R, Tmp0);
				FREEMPI(Tmp0);
				Tmp2 = MULTI(S, Aptr->V[i][j]);
				Aptr->V[i][j] = ADDI(Tmp1, Tmp2);
				FREEMPI(Tmp3);
				FREEMPI(Tmp1);
				FREEMPI(Tmp2);
			}
			for (j = 0; j < m; j++)
			{
				Tmp0P = (*Pptr)->V[k][j];
				Tmp1P = MULTI(P, (*Pptr)->V[k][j]);
				Tmp2P = MULTI(Q, (*Pptr)->V[i][j]);
				(*Pptr)->V[k][j] = ADDI(Tmp1P, Tmp2P);
				FREEMPI(Tmp1P);
				FREEMPI(Tmp2P);
				Tmp3P = (*Pptr)->V[i][j];
				Tmp1P = MULTI(R, Tmp0P);
				FREEMPI(Tmp0P);
				Tmp2P = MULTI(S, (*Pptr)->V[i][j]);
				(*Pptr)->V[i][j] = ADDI(Tmp1P, Tmp2P);
				FREEMPI(Tmp3P);
				FREEMPI(Tmp1P);
				FREEMPI(Tmp2P);
			}
			FREEMPI(P);
			FREEMPI(Q);
			FREEMPI(R);
			FREEMPI(S);
			RODIP(Aptr, *Pptr, T, k);
		}
		is_zero = 1;
		for (j = 0; j < n; j++)
		{
			if ((Aptr->V[i][j])->S)
			{
				is_zero = 0;
				break;
			}
		}
		if (is_zero)
		{
			if (i < s)
			{
				SWAP_ROWSI1(i, s, Aptr);
				SWAP_ROWSI1(i, s, *Pptr);
			}
			s--;
		}
		else
		{
			for (j = 0; j < n; j++)
			{
				if ((Aptr->V[i][T[j]])->S)
					break;
			}
			if (j != i)
			{
				tmp = T[j];
				T[j] = T[i];
				T[i] = tmp;
			}
			RODIP(Aptr, *Pptr, T, i);
			i++;
		}
	}
	*nz = s+1;
	return (T);
}
 
void RODIP(MPMATI *Aptr, MPMATI *Pptr, USI *T, USI j)
/* Kannan_Bachem subroutine: input an integer j, 0 <= j < Aptr->R such that
 * elt(Aptr, j, T[j]) != 0 and elt(Aptr, j, T[k]) = 0 for 0 <= k < j.
 * entries elt(Aptr, i, T[j]), 0 <= i < j, are reduced mod elt(Aptr, j, T[j]).
 * *Pptr is updated.
 */
{
	unsigned int i, k, m, n;
	MPI *Q;

	m = Aptr->R;
	n = Aptr->C;
	if ((Aptr->V[j][T[j]])->S == -1)
	{
		for(k = 0; k < n; k++)
		{
			if ((Aptr->V[j][k])->S != 0)
				(Aptr->V[j][k])->S = -((Aptr->V[j][k])->S);
		}
		for(k = 0; k < m; k++)
		{
			if ((Pptr->V[j][k])->S != 0)
				(Pptr->V[j][k])->S = -((Pptr->V[j][k])->S);
		}
	}	
	for (i = 0; i < j; i++)
	{
		Q = INTI(Aptr->V[i][T[j]], Aptr->V[j][T[j]]);
		Aptr = ROWSUBI(j, i, Q, Aptr); 
		Pptr = ROWSUBI(j, i, Q, Pptr); 
		FREEMPI(Q);
	}
}

MPMATI *HERMITE1P(MPMATI *Aptr, MPMATI *Pptr, MPMATI **Qptr, USI *T, USI nz)
/*
 * Takes the output of KB_ROWP() and permutes the rows to get the Hermite normal
 * form from Aptr. *Pptr is updated to **Qptr.
 */
{
	MPMATI *Tmp;
	unsigned int i, j;
	
	SORTI(T, nz);
	Tmp = COPYMATI(Aptr);
	*Qptr = COPYMATI(Pptr);
	for (i = 0; i < nz; i++)
	{
		for (j = 0; j < Aptr->C; j++)
		{
			FREEMPI(Tmp->V[i][j]);
			Tmp->V[i][j] = COPYI(Aptr->V[T[i]][j]);
		}
		for (j = 0; j < Aptr->R; j++)
		{
			FREEMPI((*Qptr)->V[i][j]);
			(*Qptr)->V[i][j] = COPYI(Pptr->V[T[i]][j]);
		}
	}
	return (Tmp);	
}

MPMATI *ADDMATI(MPMATI *Mptr, MPMATI *Nptr)
/*
 * returns *Mptr + *Nptr.
 */
{
	unsigned int i, j;
	MPMATI *Lptr;

	Lptr = BUILDMATI(Mptr->R, Mptr->C);
	for (i = 0; i <= Lptr->R - 1; i++)
		for (j = 0; j <= Lptr->C - 1; j++)
			Lptr->V[i][j] = ADDI(Mptr->V[i][j], Nptr->V[i][j]);
	return (Lptr);
}

MPMATI *SUBMATI(MPMATI *Mptr, MPMATI *Nptr)
/*
 * returns *Mptr - *Nptr.
 */
{
	unsigned int i, j;
	MPMATI *Lptr;

	Lptr = BUILDMATI(Mptr->R, Mptr->C);
	for (i = 0; i <= Lptr->R - 1; i++)
		for (j = 0; j <= Lptr->C - 1; j++)
			Lptr->V[i][j] = SUBI(Mptr->V[i][j], Nptr->V[i][j]);
	return (Lptr);
}