i9.c

calc大数库

/* i9.c */
#include <stdio.h>
#include <stdlib.h>
#include "integer.h"
#include "fun.h"

extern long int nettbytes;
/*
MPI *MAXI, *PMAXI, *QMAXI;
*/
extern int answer;

MPI **SMITHI(MPMATI *Mptr, MPMATI **Pptr, MPMATI **Qptr, USI *ptr, MPI *Eptr, USI m1, USI n1)
/*
 * returns the invariant factors of *Mptr.
 * **Pptr and **Qptr are invertible matrices such that
 * **Pptr x *Mptr x **Qptr is the Smith canonical form *Nptr below.
 * *ptr is the number of invariant factors.
 * See Rings, Modules and Linear Algebra by B. Hartley and T. O. Hawkes,
 * Chapman and Hall, 1970.
 * *Eptr is the cutoff above which we bring in MLLL. In any case MLLL is
 * done before each block is processed, as a chance of getting small entries.
 * This modification of the LLL algorithm was outlined by M. Pohst in J. 
 * Symbolic Computation (1987) 4, 123-127. m1/n1 is the LLL parameter alpha, 
 * 1/4 < alpha |= 1.
 * I initially coded this I think in 1992, along the lines of suggestions by
 * George Havas.
 */
{
	unsigned int i, k, l, m, n, p, q, j, t, z, size, r, s;
	unsigned int flag = 1, FLAG = 1;
	MPI *Q, **N, *Tmp;
	MPMATI *Nptr, *Temp;
	int isok;
	
	m = Mptr->R;
	n = Mptr->C;
	Nptr = COPYMATI(Mptr);
	*Pptr = IDENTITYI(m);
	*Qptr = IDENTITYI(n);
/*
	MAXI = MAXELTI(Mptr);
	PMAXI = MAXELTI(*Pptr);
	QMAXI = MAXELTI(*Qptr);
*/
	z = MIN((int) m - 1, (int) n - 1);
	N = (MPI **)mmalloc((z + 1) * sizeof(MPI *));
	for (i = 0; i <= z; i++)
	{
		ZEROTESTI(0, Nptr, &r, &s);
		if (r == Nptr->R && s == Nptr->C)
			break;
		flag = HAVAS_PIVOTI(Nptr, &p, &q, r, s);
		if (p > 0)/* must have a non-zero row 1 for MLLL */
		{
			Nptr = SWAP_ROWSI1(0, p, Nptr);
			*Pptr = SWAP_ROWSI1(i, p + i, *Pptr);
		}
		if (q > 0)
		{
			Nptr = SWAP_COLSI1(0, q, Nptr);
			*Qptr = SWAP_COLSI1(i, q + i, *Qptr);
		}
		if (flag)
			FLAG = 1;
		if (!flag && FLAG)
		{
			size = Nptr->R;
			Temp = Nptr;
			Nptr = BASIS_REDUCTION(Nptr, Pptr, i, m1, n1);
			/*Nptr = BASIS_REDUCTION0(Nptr, 3, 4);*/
			FREEMATI(Temp);
			ZEROTESTI(0, Nptr, &r, &s);
			if (r == Nptr->R && s == Nptr->C)
				break;
			flag = HAVAS_PIVOTI(Nptr, &p, &q, r, s);
			if (p > 0)
			{
				Nptr = SWAP_ROWSI1(0, p, Nptr);
				*Pptr = SWAP_ROWSI1(i, p + i, *Pptr);
			}
			if (q > 0)
			{
				Nptr = SWAP_COLSI1(0, q, Nptr);
				*Qptr = SWAP_COLSI1(i, q + i, *Qptr);
			}
			FLAG = 0;
			size = size - Nptr->R;
			m = m - size;
			z = MIN((int) m - 1, (int) n - 1);
		}
		/* Now elt(Nptr, 0, 0) != 0. */
		while (1)
		{
 		    while(1)
		    { /* PA bit */
			isok = 1;
			for(j = 0; j < Nptr->C; j++)
			{
				if ((Nptr->V[0][j])->S == -1)
				{
					for(t = 0; t < Nptr->R; t++)
						(Nptr->V[t][j])->S = - (Nptr->V[t][j])->S;
					for(t = 0; t < n; t++)
						((*Qptr)->V[t][j + i])->S = - ((*Qptr)->V[t][j + i])->S;
				}
			}
			/* Now all non-zero elements in row 0 are positive. */
			while (1)
			{
				l = ROWSEEKI(0, Nptr);
				if (l == Nptr->C)
					 break;
				else
				{
					/* Case 1, p. 112. */
					Q = INT(Nptr->V[0][l], Nptr->V[0][0]);
					Nptr = COLSUBI(0, l, Q, Nptr); 
/*
					MAXI = UPDATEMAXI(MAXI, Nptr);
*/
					*Qptr = COLSUBI(i, l + i, Q, *Qptr); 
/*
					QMAXI = UPDATEMAXI(QMAXI, *Qptr);
*/
					FREEMPI(Q);
					Nptr = SWAP_COLSI1(0, l, Nptr); 
					*Qptr = SWAP_COLSI1(i, l + i, *Qptr); 
					/* This reduces the size of elt(Nptr, 0, 0). */
				}
			}
			for(j = 1; j < Nptr->R; j++)
			{
				if ((Nptr->V[j][0])->S == -1)
				{
					for(t = 0; t < Nptr->C; t++)
						(Nptr->V[j][t])->S = - (Nptr->V[j][t])->S;
					for(t = 0; t < (*Pptr)->R; t++)
						((*Pptr)->V[j + i][t])->S = - ((*Pptr)->V[j + i][t])->S;
				}
			}
/*
Now all non-zero elements in column 0 are positive.
*/
			while (1)
			{
				k = COLSEEKI(0, Nptr);
				if (k == Nptr->R)
					break;
				else
				{
					isok = 0;
					/* Case 2, p. 112. */
					Q = INT0(Nptr->V[k][0], Nptr->V[0][0]);
					Nptr = ROWSUBI(0, k, Q, Nptr); 
/*
					MAXI = UPDATEMAXI(MAXI, Nptr);
*/
					*Pptr = ROWSUBI(i, k + i, Q, *Pptr);
/*
					PMAXI = UPDATEMAXI(PMAXI, *Pptr);
*/
					FREEMPI(Q);
					Nptr = SWAP_ROWSI1(0, k, Nptr); 
					*Pptr = SWAP_ROWSI1(i, k + i, *Pptr);
					/* This reduces the size of elt(Nptr, 0, 0). */
				}
			}
			/* now l = Nptr->C and k = Nptr->R */
			if (isok)
				break;
		   } /*PA bit*/
/*
Now (Case 3, p.112) elt(Nptr, 0, 0) divides each element to its right as well as below it.
*/
		for (l = 1; l < Nptr->C; l++)
		{
			if ((Nptr->V[0][l])->S != 0)
			{
	/* NOTE: Swapping above may have made elt(Nptr, 0, l) <0 */
				Q = INT(Nptr->V[0][l], Nptr->V[0][0]);
				Nptr = COLSUBI(0, l, Q, Nptr); 
/*
				MAXI = UPDATEMAXI(MAXI, Nptr);
*/
				*Qptr = COLSUBI(i, l + i, Q, *Qptr); 
/*
				QMAXI = UPDATEMAXI(QMAXI, *Qptr);
*/
				FREEMPI(Q);
			}
		}
		for (k = 1; k < Nptr->R; k++)
		{
			if ((Nptr->V[k][0])->S != 0)
			{
				Q = INT0(Nptr->V[k][0], Nptr->V[0][0]);
				FREEMPI(Nptr->V[k][0]);
				Nptr->V[k][0] = ZEROI();
				*Pptr = ROWSUBI(i, k + i, Q, *Pptr);
/*
				PMAXI = UPDATEMAXI(PMAXI, *Pptr);
*/
				FREEMPI(Q);
			}
		}
/*
Now each element to the right of and below elt(Nptr, 0, 0) is zero.
*/
		k = MATSEEKI(0, Nptr);
		if (k != Nptr->R)
		{
			for (t = 0; t < Nptr->C; t++)
			{
				if ((Nptr->V[k][t])->S != 0)
				{
					Tmp = Nptr->V[0][t];
					Nptr->V[0][t] = ADDI(Nptr->V[0][t], Nptr->V[k][t]);
/*
					MAXI = UPDATEMAXI(MAXI, Nptr);
*/
					FREEMPI(Tmp);
				}
			}
			for (t = 0; t < Mptr->R; t++)
			{
				if (((*Pptr)->V[k + i][t])->S != 0)
				{
					Tmp = (*Pptr)->V[i][t];
					(*Pptr)->V[i][t] = ADDI((*Pptr)->V[i][t], (*Pptr)->V[k + i][t]);
/*
					PMAXI = UPDATEMAXI(PMAXI, *Pptr);
*/
					FREEMPI(Tmp);
				}
			}
		}
		else
			break;
	}
/*
Now elt(Nptr, 0, 0) divides all elements elt(Nptr, k, l), k,l > 0.
*/
		N[i] = COPYI(Nptr->V[0][0]);

		if (i < z)
		{
			Temp = Nptr;
			Nptr = DELETE_ROW1I(Nptr);	
			FREEMATI(Temp);
			Temp = Nptr;
			Nptr = DELETE_COL1I(Nptr);	
			FREEMATI(Temp);
		}
		size = Nptr->R;
		Tmp = MAXELTI(Nptr);
		printf("found invariant factor d[%u] = ", i + 1); PRINTI(N[i]); printf("\n");
/*
		printf("MAXI = "); PRINTI(MAXI); printf("\n");
		printf("PMAXI = "); PRINTI(PMAXI); printf("\n");
		printf("QMAXI = "); PRINTI(QMAXI); printf("\n");
*/
		printf("MAX_ENTRY(Nptr) = "); PRINTI(Tmp); printf("\n");
		l = RSV(Tmp, Eptr);
		FREEMPI(Tmp);
		if (l == 1 && !flag)/* we don't want to use MLLL if there are
1's or -1's still present */
		{
			ZEROTESTI(0, Nptr, &r, &s);
			if (r > 0)/* must have a non-zero row 1 for MLLL */
			{
				printf("r = %u > 0\n", r);
				Nptr = SWAP_ROWSI1(0, r, Nptr);
				*Pptr = SWAP_ROWSI1(i + 1, r + i + 1, *Pptr);
			}
			Temp = Nptr;
			Nptr = BASIS_REDUCTION(Nptr, Pptr, i+1, m1, n1);
			/*Nptr = BASIS_REDUCTION0(Nptr, 3, 4);*/
			FREEMATI(Temp);
		}
		size = size - Nptr->R;
		m = m - size;
		z = MIN((int) m - 1, (int) n - 1);
	}
	FREEMATI(Nptr);
	*ptr = i++;
/*
	printf("MAXI = "); PRINTI(MAXI); printf("\n");
	printf("PMAXI = "); PRINTI(PMAXI); printf("\n");
	printf("QMAXI = "); PRINTI(QMAXI); printf("\n");
	FREEMPI(MAXI);
	FREEMPI(PMAXI);
	FREEMPI(QMAXI);
*/
	return (N);
}

MPI *ROWSUMI(MPMATI *Mptr, USI i)
{
	unsigned int j;
	MPI *S, *Temp1I, *Temp2I;

	S = ZEROI();
	for (j = 0; j < Mptr->C; j++)
	{
		Temp1I = ABSI(Mptr->V[i][j]);
		Temp2I = S;
		S = ADD0I(Temp1I, S);
		FREEMPI(Temp1I);
		FREEMPI(Temp2I);
	}
	return (S);
}

MPI *COLSUMI(MPMATI *Mptr, USI j)
{
	unsigned int i;
	MPI *S, *Temp1I, *Temp2I;

	S = ZEROI();
	for (i = 0; i < Mptr->R; i++)
	{
		Temp1I = ABSI(Mptr->V[i][j]);
		Temp2I = S;
		S = ADD0I(Temp1I, S);
		FREEMPI(Temp1I);
		FREEMPI(Temp2I);
	}
	return (S);
}

MPI *MINELTI(MPMATI *Mptr)
/* returns the minimum element of *Mptr */
{
	unsigned int i, j, I, J, flag = 0;

	I = J = 0;
	for (i = 0; i < Mptr->R; i++)
	{
		for (j = 0; j < Mptr->C; j++)
		{
			if ((Mptr->V[i][j])->S)
			{
				I = i;
				J = j;
				flag = 1;
				break;
			}
		}
		if (flag)
			break;
	}
	for (i = I; i < Mptr->R; i++)
	{
		for (j = 0; j < Mptr->C; j++)
		{
			if (RSV(Mptr->V[I][J], Mptr->V[i][j]) == 1 && (Mptr->V[i][j])->S)
			{
				I = i;
				J = j;
			}
		}
	}
	return (ABSI(Mptr->V[I][J]));
}

void MINELTI1(MPMATI *Mptr, USI *iptr, USI *jptr)
/* returns the position of the minimum element of *Mptr */
{
	unsigned int i, j, I, J, flag = 0;

	I = J = 0;
	for (i = 0; i < Mptr->R; i++)
	{
		for (j = 0; j < Mptr->C; j++)
		{
			if ((Mptr->V[i][j])->S)
			{
				I = i;
				J = j;
				flag = 1;
				break;
			}
		}
		if (flag)
			break;
	}
	for (i = 0; i < Mptr->R; i++)
	{
		for (j = 0; j < Mptr->C; j++)
		{
			if (RSV(Mptr->V[I][J], Mptr->V[i][j]) == 1 && (Mptr->V[i][j])->S)
			{
				I = i;
				J = j;
			}
		}
	}
	*iptr = I;
	*jptr = J;
	return;
}

MPI *MAXELTI(MPMATI *Mptr)
/* returns the maximum element of *Mptr */
{
	unsigned int i, j, I = 0, J = 0;

	for (i = 0; i < Mptr->R; i++)
	{
		for (j = 0; j < Mptr->C; j++)
		{
			if (RSV(Mptr->V[i][j], Mptr->V[I][J]) == 1)
			{
				I = i;
				J = j;
			}
		}
	}
	return (ABSI(Mptr->V[I][J]));
}

void MAXELTI1(MPMATI *Mptr, USI *iptr, USI *jptr)
/* returns the position of the maximum element of *Mptr */
{
	unsigned int i, j, I = 0, J = 0;

	for (i = 0; i < Mptr->R; i++)
	{
		for (j = 0; j < Mptr->C; j++)
		{
			if (RSV(Mptr->V[i][j], Mptr->V[I][J]) == 1)
			{
				I = i;
				J = j;
			}
		}
	}
	*iptr = I;
	*jptr = J;
	return;
}

MPI *UPDATEMAXI(MPI *S, MPMATI *Mptr)
{
	MPI *TempI;

	TempI = MAXELTI(Mptr);
	if (RSV(TempI, S) == 1)
	{
		FREEMPI(S);
		return(TempI);
	}	
	else
	{
		FREEMPI(TempI);
		return (S);
	}
}

unsigned int HAVAS_PIVOTI(MPMATI *Mptr, USI *rptr, USI *cptr, USI r, USI s)
{
	MPI **RS, **CS, *S, *Temp1I, *Temp2I, *TempI;
	unsigned int i, j, I, J, m, n, flag;

	m = Mptr->R;
	n = Mptr->C;
	RS = (MPI **)mmalloc(m * sizeof(MPI *));
	CS = (MPI **)mmalloc(n * sizeof(MPI *));
	for (i = 0; i < m; i++)
	{
		TempI = ROWSUMI(Mptr, i);
		if (TempI->S == 0)
		{
			FREEMPI(TempI);
			RS[i] = MINUS_ONEI();
		}
		else
		{
			Temp1I = ONEI();
			RS[i] = SUB0I(TempI, Temp1I);
			FREEMPI(TempI);
			FREEMPI(Temp1I);
		}
	}
	for (j = 0; j < Mptr->C; j++)
	{
		TempI = COLSUMI(Mptr, j);
		if (TempI->S == 0)
		{
			FREEMPI(TempI);
			CS[j] = MINUS_ONEI();
		}
		else
		{
			Temp1I = ONEI();
			CS[j] = SUB0I(TempI, Temp1I);
			FREEMPI(TempI);