i9.c

calc大数库

			FREEMPI(Temp1I);
		}
	}
	S = MULTI(RS[r], CS[s]); /* elt(r,s) is the 1st non-zero elt of Nptr */
	I = r; J = s;
	flag = 0;
	/* first minimize RS[i]*CS[j] over ones or minus-ones */
	for (i = r; i < m; i++)
	{
		for (j = 0; j < n; j++)
		{
			if ((Mptr->V[i][j])->D == 0 && (Mptr->V[i][j])->V[0] == 1)
			{
				flag = 1;
				Temp1I = S;
				Temp2I = MULTI(RS[i], CS[j]);
				if(RSV(Temp1I, Temp2I) == 1)
				{
					S = Temp2I;
					FREEMPI(Temp1I);
					I = i; J = j;
				}
				else
					FREEMPI(Temp2I);	
			}
		}
	}
	if (!flag) /* minimize RS[i]*CS[j] over a matrix with no ones or minus-ones */
	{
		for (i = r; i < m; i++)
		{
			for (j = 0; j < n; j++)
			{
				if ((Mptr->V[i][j])->S)
				{
					Temp1I = S;
					Temp2I = MULTI(RS[i], CS[j]);
					if(RSV(Temp1I, Temp2I) == 1)
					{
						S = Temp2I;
						FREEMPI(Temp1I);
						I = i; J = j;
					}
					else
						FREEMPI(Temp2I);	
					
				}
			}
		}
	}
	*rptr = I;
	*cptr = J;
	for (i = 0; i < m; i++)
		FREEMPI(RS[i]);
	for (j = 0; j < n; j++)
		FREEMPI(CS[j]);
	ffree((char *)RS, m * sizeof(MPI *));
	ffree((char *)CS, n * sizeof(MPI *));
	FREEMPI(S);
	return (flag);
}

void ZEROTESTI(USI i, MPMATI *Mptr, USI *p, USI *q)
/*
 * returns the least indices *p >= i, *q >= i such that elt(Mptr, *p, *q) != 0.
 * otherwise *p = Mptr->R and *q = Mptr->C.
 */
{
	unsigned int k, l;
	
	for (k = i; k <= Mptr->R - 1; k++)
	{
		for (l = i; l <= Mptr->C - 1; l++)
		{
			if ((Mptr->V[k][l])->S != 0)
			{
				*p = k;
				*q = l;
				return;
			}
		}
	}
	*p = Mptr->R;
	*q = Mptr->C;
	return;
}

unsigned int MATSEEKI(USI i, MPMATI *Mptr)
/*
 * returns the least row index k, i < k <= Mptr->R - 1, such that 
 * elt(Mptr, k, l) is not divisible by elt(Mptr, i, i), for some l,
 * i < l <= Mptr->C - 1; otherwise returns Mptr->R.
 */
{
	int s;
	unsigned int k, l;
	MPI *P, *Q;
	
	Q = ABSI(Mptr->V[i][i]);
	for (k = i + 1; k <= Mptr->R - 1; k++)
	{
		for (l = i + 1; l <= Mptr->C - 1; l++)
		{
			P = MOD(Mptr->V[k][l], Q);
			s = P->S;
			FREEMPI(P);
			if (s != 0)
			{
				FREEMPI(Q);
				return (k);
			}
		}
	}
	FREEMPI(Q);
	return (Mptr->R);
}

unsigned int ROWSEEKI(USI i, MPMATI *Mptr)
/*
 * returns the least integer l > i such that elt(Mptr, i, l) is not divisible
 * by elt(Mptr, i, i), otherwise returns Mptr->C.
 */
{
	int s;
	unsigned int l;
	MPI *P;

	for (l = i + 1; l <= Mptr->C - 1; l++)
	{
		if ((Mptr->V[i][l])->S != 0)
		{
			P = MOD0(Mptr->V[i][l], Mptr->V[i][i]);
			s = P->S;
			FREEMPI(P);
			if (s != 0)
				return (l);
		}
	}
	return (Mptr->C);
}

unsigned int COLSEEKI(USI i, MPMATI *Mptr)
/*
 * returns the least integer k > i such that Mptr->V[k][i] is not divisible
 * by Mptr->V[i][i], otherwise returns Mptr->R.
 */
{
	int s;
	unsigned int k;
	MPI *P;

	for (k = i + 1; k <= Mptr->R - 1; k++)
	{
		if ((Mptr->V[k][i])->S != 0)
		{
			P = MOD0(Mptr->V[k][i], Mptr->V[i][i]);
			s = P->S;
			FREEMPI(P);
			if (s != 0)
				return (k);
		}
	}
	return (Mptr->R);
}

MPMATI *ROWSUBI0(USI p, USI q, MPI *Aptr, MPMATI *Mptr)
/*
 * subtracts *Aptr times the p-th row of *Mptr from the q-th.
 * 0 <= p, q <= Mprt->R - 1, where elt(Mptr, i, j) = 0 if j < p.
 */
{
	MPI *Y, *Temp;
	unsigned int l;
	if (Aptr->S != 0)
	{
		for (l = p; 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 *COLSUBI0(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, where elt(Mptr, i, j) = 0 if i < p.
 */
{
	MPI *Y, *Temp;
	unsigned int k;

	if (Aptr->S != 0)
	{
		for (k = p; 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);
}

MPI **SMITHI1(MPMATI *Mptr, USI *ptr, MPI *Eptr, USI m1, USI n1)
/*
 * returns the invariant factors of *Mptr.
 * *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.
 */
{
	unsigned int i, k, l, m, n, p, q, j, t, z, size, r, s;
	unsigned int flag = 1;
	MPI *Q, **N, *Tmp;
	MPMATI *Nptr, *Temp;
	int isok;
	
	m = Mptr->R;
	n = Mptr->C;
	Nptr = COPYMATI(Mptr);
/*
	MAXI = MAXELTI(Mptr);
*/
	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);
		if (q > 0)
			Nptr = SWAP_COLSI1(0, q, Nptr);
/*
		if (flag)
			FLAG = 1;
		if (!flag && FLAG)
		{
			size = Nptr->R;
			Temp = Nptr;
			printf("FLAG = 1 and flag = 0: doing MLLL:\n");
			Nptr = BASIS_REDUCTION0(Nptr, m1, n1);
			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);
			if (q > 0)
				Nptr = SWAP_COLSI1(0, q, Nptr);
			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;
				}
			}
			/* 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);
*/
					FREEMPI(Q);
					Nptr = SWAP_COLSI1(0, l, Nptr); 
					/* 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;
				}
			}
/*
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);
*/
					FREEMPI(Q);
					Nptr = SWAP_ROWSI1(0, k, Nptr); 
					/* 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);
*/
				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();
				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);
				}
			}
		}
		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("MAX_ENTRY(Nptr) = "); PRINTI(Tmp); printf("\n");
		printf("MAXI has %u digits\n", LENGTHI(MAXI));
*/
		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);
			}
			Temp = Nptr;
			Nptr = BASIS_REDUCTION0(Nptr, m1, n1);
			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");
	FREEMPI(MAXI);
*/
	return (N);
}

MPMATI *DELETE_ROW1I(MPMATI *Mptr)
/*
 * deletes the first row of *Mptr.
 */
{
	unsigned int i, j, m, n;
	MPMATI *Nptr;

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

MPMATI *DELETE_COL1I(MPMATI *Mptr)
/*
 * deletes  the first column of *Mptr.
 */
{
	unsigned int i, j, m, n;
	MPMATI *Nptr;

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

unsigned int PIVOTI(MPMATI *Mptr, USI *rptr, USI *cptr, USI r, USI s)
{
	unsigned int i, j, m, n;

	m = Mptr->R;
	n = Mptr->C;
	*rptr = r; *cptr = s;
	for (i = r; i < m; i++)
	{
		for (j = 0; j < n; j++)
		{
			if ((Mptr->V[i][j])->D == 0 && (Mptr->V[i][j])->V[0] == 1)
			{
				*rptr = i; *cptr = j;
				return (1);
			}
		}
	}
	return (0);
}