Another LLL implementation

This is from Calc; see
http://www.maths.uq.edu.au/~krm/krm_calc.html#[49]

/* The Pohst Algorithm updating only from where necessary */
#include <stdio.h>
#include <stdlib.h>
#include "integer.h"
#include "fun.h"
extern MPI *MAXI, *PMAXI;

extern unsigned int MLLLVERBOSE;
extern unsigned int HERMITEVERBOSE;
unsigned int GCDFLAG;

MPMATI *BASIS_REDUCTION(MPMATI *Bptr, MPMATI **Eptr, USI rowstage,
USI m1, USI n1)
/*
* Input: *Bptr, a matrix of MPI's, whose first row is not zero.
* Output: a pointer to an MPMATI whose rows form a reduced basis for
* the lattice spanned by the rows of *Bptr. This basis is reduced in
the
* sense of the paper "Factoring polynomials with rational
coefficients" by
* A. K. Lenstra, H. W. Lenstra and L. Lovasz, Math. Ann. 261, 515-
534 (1982)
* using the modified version in "Solving exponential Diophantine
equations
* using lattice basis reduction algorithms" by B. M. M. De Weger, J.
No. Theory
* 26, 325-367 (1987). A change of basis matrix **Eptr is also
returned.
* De Weger's algorithm has been changed to cater for arbitrary
matrices. The
* the rows are now in general linearly dependent.
* We use the fact that the Gram Schmidt process detects the first
row
* which is a linear combination of the preceding rows. We employ a
modification
* of the LLL algorithm outlined by M. Pohst in J. Symbolic
Computation (1987)4,
* 123-127. We call this the MLLL algorithm.
* The last sigma rows of the matrix **Eptr are relation vectors.
* m1 / n1 is usually taken to be 3 / 4.
*/
{
unsigned int i, k, l, n, m, t, flag = 0, Flag = 0;
unsigned int flagg, beta, K1 = 0, tau = 2, sigma = 0, rho;
MPI **D, *X, *Y, *Z, *H, *Tmp, *R, *M1, *N1;
MPMATI *C, *L, *B1ptr;

m = Bptr->C;
n = Bptr->R;

/* We initial Eptr outside the function whenever we call the
function. */
/* This is because we have to do so in SMITH(). */

/* MAXI = MAXELTI(Bptr);
PMAXI = MAXELTI(*Eptr); */
B1ptr = COPYMATI(Bptr);
D = (MPI **)mmalloc((1 + n) * sizeof(MPI *));
D[0] = ONEI();
for (i = 1; i <= n; i++)
D[i] = ZEROI();
C = ZEROMNI(n, m);
L = ZEROMNI(n, n);

found:
n = B1ptr->R;
i = (K1 == 0) ? 1 : K1;
/* K1 = no. of consecutive rows of *B1ptr that don't need
updating
for the Gram Schmidt process. */
while (i <= B1ptr->R)
{
BASIS_UPDATE(i, m, &C, &L, B1ptr, D);
flag = 1;
for (t = 0; t < m; t++)
{
if (!EQZEROI(C->V[i - 1][t]))
{
flag = 0;
break;
}
}
if (flag)
break;
X = ZEROI();
for (t = 0; t < m; t++)
{
H = MULTI(C->V[i - 1][t], C->V[i - 1][t]);
Tmp = X;
X = ADDI(X, H);
FREEMPI(Tmp);
FREEMPI(H);
}
FREEMPI(D[i]);
D[i] = INT(X, D[i - 1]);
FREEMPI(X);
i++;
if (MLLLVERBOSE)
printf("i = %u\n", i);
}
beta = (flag) ? i : i - 1;
rho = K1 = i - 1;
if (MLLLVERBOSE)
printf("completed updating the basis\n");
/* Here K1 = no. of LI rows in *B1ptr found by Gram Schmidt process.
flag = 0 means all the rho = beta rows of *B1ptr are LI;
flag = 1 means that the first rho = beta - 1 rows of *B1ptr are
LI, but the
beta-th row is a LC of the preceding rows. So beta = number of
rows of *B1ptr
currently being examined by the LLL algorithm. */
k = tau;
if (MLLLVERBOSE)
printf("beta = %u\n", beta);
M1 = CHANGE(m1);
N1 = CHANGE(n1);
while (k <= beta)
{
if (MLLLVERBOSE)
printf("beta - k = %u\n", beta -k);
l = k - 1;
Flag = STEP4(k, l, &L, &B1ptr, Eptr, D, rowstage);
if (Flag)/* STEP 9 of POHST. */
{
sigma++;
if (MLLLVERBOSE)
printf("relation vector number %u
found\n", sigma);
tau = k++;
goto found;
}
X = MULTI(D[k - 2], D[k]);
Y = MULTI(D[k - 1], D[k - 1]);
Tmp = Y;
Y = MULT_I(Y, m1);
FREEMPI(Tmp);
R = MULTI(L->V[k - 1][k - 2], L->V[k - 1][k - 2]);
Z = ADD0I(X, R);
Tmp = Z;
Z = MULT_I(Z, n1);
FREEMPI(Tmp);
if (RSV(Y, Z) == 1)/*& STEP 5 of POHST. */
{
flagg = 0;
if (EQZEROI(D[k]) && EQZEROI(R))
{/* CASE B=0 of STEP 7 of POHST. */
FREEMPI(D[k - 1]);
D[k - 1] = ZEROI();
STEP8(k, &B1ptr, &L, Eptr, rowstage);
if (k - 1 < K1)
K1 = k - 1;
/* The swap may have changed 2nd last
row */
/* of *B1ptr. */
for (t = 0; t < m; t++)
{
FREEMPI(C->V[k - 2][t]);
C->V[k - 2][t] = ZEROI();
}
beta--;
flagg = 1;
if (k > 2)
k--;
FREEMPI(X);
FREEMPI(Y);
FREEMPI(R);
FREEMPI(Z);
continue;
}
if (flagg == 0)
{
for (i = k + 1; i <= beta; i++)
STEP7(i, k, &L, D);
}
STEP8(k, &B1ptr, &L, Eptr, rowstage);
if (k - 2 < K1)
K1 = k - 2;
/* swap will change last two rows of *B1ptr.
*/
FREEMPI(R);
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
{ /* STEP 6 of POHST. */
FREEMPI(R);
FREEMPI(X);
FREEMPI(Y);
for (l = k - 2; l >= 1; l--)
{
Flag = STEP4(k, l, &L, &B1ptr, Eptr,
D, rowstage);
if (Flag)
{
FREEMPI(Z);
sigma++;
if (MLLLVERBOSE)
printf("relation vector
number %u found\n", sigma);
tau = k++;
goto found; /* STEP 9 of
POHST. */
}
}
k++;
}
FREEMPI(Z);
}
FREEMPI(M1);
FREEMPI(N1);
FREEMATI(C);

printf("L = \n");
PRINTMATI(0,L->R-1,0,L->C-1,L);
for (i = 0; i <= Bptr->R; i++)
{
printf("D[%u] = ", i);PRINTI(D[i]);printf(", ");
}
printf("\n");
FREEMATI(L);
for (i = 0; i <= Bptr->R; i++)
FREEMPI(D[i]);
ffree((char *)D, (1 + Bptr->R) * sizeof(MPI *));
if (MLLLVERBOSE)
{
printf("number of basis vectors found = %u ;\n", rho);
printf("number of relation vectors found = %u .\n",
sigma);
}
return (B1ptr);
}

unsigned int STEP4(k, l, Lptr, Bptr, Eptr, D, i)
/*
* updates *Lptr, *Bptr and *Eptr.
* returns 1 if row k of *Bptr becomes zero, returns zero otherwise.
*/
unsigned int k, l, i;
MPI *D[];
MPMATI **Lptr, **Bptr, **Eptr;
{
unsigned int j, flag = 1, t, m, n;
MPI *X, *Y, *R, *Tmp;
MPMATI *TmpMATI;

m = (*Bptr)->C;
n = (*Eptr)->R;
Y = MULT_I((*Lptr)->V[k - 1][l - 1], 2);
if (RSV(Y, D[l]) == 1)
{
R = NEAREST_INTI((*Lptr)->V[k - 1][l - 1], D[l]);
X = MINUSI(R);
TmpMATI = *Bptr;
*Bptr = ADD_MULT_ROWI(l - 1, k - 1, X, *Bptr);
/*
MAXI = UPDATEMAXI(MAXI, *Bptr);
*/
FREEMATI(TmpMATI);
TmpMATI = *Eptr;
*Eptr = ADD_MULT_ROWI(l + i - 1, k + i - 1, X, *Eptr);
/*
PMAXI = UPDATEMAXI(PMAXI, *Eptr);
*/
FREEMATI(TmpMATI);
FREEMPI(X);
for (j = 1; j < l; j++)
{
X = MULTI((*Lptr)->V[l - 1][j - 1], R);
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], R);
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(R);
}
for (t = 0; t < m; t++)
{
if (!EQZEROI((*Bptr)->V[k - 1][t]))
{
flag = 0;
break;
}
}
if (flag)
{
TmpMATI = *Bptr;
*Bptr = DELETE_ROWI(k, *Bptr);
FREEMATI(TmpMATI);
for (j = k - 1; j < n - i - 1; j++)
*Eptr = SWAP_ROWSI1(j + i, j + i + 1, *Eptr);
}
FREEMPI(Y);
return (flag);
}

void STEP8(USI k, MPMATI **B1ptr, MPMATI **Lptr, MPMATI **Eptr, USI i)
{
MPI *T;

*B1ptr = SWAP_ROWSI1(k - 2, k - 1, *B1ptr);
*Eptr = SWAP_ROWSI1(k + i - 2, k + i - 1, *Eptr);
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();
return;
}

void STEP7(USI i, USI k, MPMATI **Lptr, MPI *D[])
{
MPI *X1, *X2, *X3, *Y1, *Y2, *Tmp;

X1 = MULTI((*Lptr)->V[i - 1][k - 2], (*Lptr)->V[k - 1][k -
2]);
Y1 = MULTI((*Lptr)->V[i - 1][k - 1], D[k - 2]);
Tmp = Y1;
Y1 = ADDI(Y1, X1);
FREEMPI(Tmp);
FREEMPI(X1);
X2 = MULTI((*Lptr)->V[i - 1][k - 2], D[k]);
X3 = MINUSI((*Lptr)->V[k - 1][k - 2]);
Y2 = MULTI((*Lptr)->V[i - 1][k - 1], X3);
FREEMPI(X3);
Tmp = Y2;
Y2 = ADDI(Y2, X2);
FREEMPI(Tmp);
FREEMPI(X2);
FREEMPI((*Lptr)->V[i - 1][k - 2]);
(*Lptr)->V[i - 1][k - 2] = INT(Y1, D[k - 1]);
FREEMPI((*Lptr)->V[i - 1][k - 1]);
(*Lptr)->V[i - 1][k - 1] = INT(Y2, D[k - 1]);
FREEMPI(Y1);
FREEMPI(Y2);
return;
}

void BASIS_UPDATE(USI i, USI m, MPMATI **Cptr, MPMATI **Lptr, MPMATI
*B1ptr, MPI *D[])
{
unsigned int j, k, t;
MPI *X, *Tmp, *X1, *X2, *H;

for (k = 1; k <= m; k++)
{
FREEMPI((*Cptr)->V[i - 1][k - 1]);
(*Cptr)->V[i - 1][k - 1] = COPYI(B1ptr->V[i - 1][k -
1]);
}
for (j = 1; j < i; j++)
{
X = ZEROI();
for (t = 0; t < m; t++)
{
if (((*Cptr)->V[j - 1][t])->S != 0 && (B1ptr-
>V[i - 1][t])->S != 0)
{
H = MULTI((*Cptr)->V[j - 1][t], B1ptr-
>V[i - 1][t]);
Tmp = X;
X = ADDI(X, H);
FREEMPI(H);
FREEMPI(Tmp);
}
}
FREEMPI((*Lptr)->V[i - 1][j - 1]);
(*Lptr)->V[i - 1][j - 1] = X;
for (t = 0; t < m; t++)
{
if (((*Cptr)->V[i - 1][t])->S == 0)
X1 = ZEROI();
else
X1 = MULTI((*Cptr)->V[i - 1][t], D
[j]);
if (((*Cptr)->V[j - 1][t])->S != 0 &&
((*Lptr)->V[i - 1][j - 1])->S != 0)
X2 = MULTI((*Cptr)->V[j - 1][t],
(*Lptr)->V[i - 1][j - 1]);
else
X2 = ZEROI();
Tmp = X1;
X1 = SUBI(X1, X2);
FREEMPI(Tmp);
FREEMPI(X2);
FREEMPI((*Cptr)->V[i - 1][t]);
(*Cptr)->V[i - 1][t] = INT(X1, D[j - 1]);
FREEMPI(X1);
}
}
return;
}

void CSWAP_UPDATE(USI k, USI m, MPI *S, MPMATI **Cptr, MPI *D[])
{
unsigned int t;
MPI *Tmp1, *Tmp2, *Tmp3, *Tmp4;

for (t = 0; t < m; t++)
{
Tmp1 = MULTI((*Cptr)->V[k - 1][t], D[k - 2]);
Tmp2 = MULTI((*Cptr)->V[k - 2][t], S);
Tmp3 = ADDI(Tmp1, Tmp2);
FREEMPI(Tmp1);
FREEMPI(Tmp2);

Tmp1 = MULTI((*Cptr)->V[k - 2][t], D[k]);
Tmp2 = MULTI((*Cptr)->V[k - 1][t], S);
Tmp4 = SUBI(Tmp1, Tmp2);
FREEMPI(Tmp1);
FREEMPI(Tmp2);
FREEMPI((*Cptr)->V[k - 2][t]);
FREEMPI((*Cptr)->V[k - 1][t]);
(*Cptr)->V[k - 2][t] = INT(Tmp3, D[k - 1]);
(*Cptr)->V[k - 1][t] = INT(Tmp4, D[k - 1]);
FREEMPI(Tmp3);
FREEMPI(Tmp4);
}
return;
}

MPMATI *BASIS_REDUCTION0(MPMATI *Bptr, USI m1, USI n1)
/*
* Input: *Bptr, a matrix of MPI's, whose first row is not zero.
* Output: a pointer to an MPMATI whose rows form a reduced basis for
* the lattice spanned by the rows of *Bptr. This basis is reduced in
the
* sense of the paper "Factoring polynomials with rational
coefficients" by
* A. K. Lenstra, H. W. Lenstra and L. Lovasz, Math. Ann. 261, 515-
534 (1982)
* using the modified version in "Solving exponential Diophantine
equations
* using lattice basis reduction algorithms" by B. M. M. De Weger, J.
No. Theory
* 26, 325-367 (1987). No change of basis matrix is returned.
* De Weger's algorithm has been changed to cater for arbitrary
matrices. The
* the rows are now in general linearly dependent.
* We use the fact that the Gram Schmidt process detects the first
row
* which is a linear combination of the preceding rows. We employ a
modification
* of the LLL algorithm outlined by M. Pohst in J. Symbolic
Computation (1987)4,
* 123-127. We call this the MLLL algorithm.
* If we are using this algorithm to find small multipliers for the
extended
* gcd problem, GCDFLAG is set in EXTGCD() and gcdflag is set below.
* m1 / n1 is usually taken to be 3 / 4.
*/
{
unsigned int i, k, l, n, m, t, flag = 0, Flag = 0, gcdflag =
0;
unsigned int flagg, beta, K1 = 0, tau = 2, sigma = 0, rho;
MPI **D, *X, *Y, *Z, *H, *Tmp, *R, *M1, *N1;
MPMATI *C, *L, *B1ptr;
unsigned int norig;

Z = NULL;
m = Bptr->C;
n = Bptr->R;
norig = n;
B1ptr = COPYMATI(Bptr);
D = (MPI **)mmalloc((1 + n) * sizeof(MPI *));
D[0] = ONEI();
for (i = 1; i <= n; i++)
D[i] = ZEROI();
C = ZEROMNI(n, m);
L = ZEROMNI(n, n);

found:
n = B1ptr->R;
i = (K1 == 0) ? 1 : K1;
/* K1 = no. of consecutive rows of *B1ptr that don't need
updating
for the Gram Schmidt process. */
while (i <= B1ptr->R)
{
BASIS_UPDATE(i, m, &C, &L, B1ptr, D);
flag = 1;
for (t = 0; t < m; t++)
{
if (!EQZEROI(C->V[i - 1][t]))
{
flag = 0;
break;
}
}
if (flag)
break;
X = ZEROI();
for (t = 0; t < m; t++)
{
H = MULTI(C->V[i - 1][t], C->V[i - 1][t]);
Tmp = X;
X = ADDI(X, H);
FREEMPI(Tmp);
FREEMPI(H);
}
FREEMPI(D[i]);
D[i] = INT(X, D[i - 1]);
FREEMPI(X);
i++;
if (MLLLVERBOSE)
printf("i = %u\n", i);
}
beta = (flag) ? i : i - 1;
rho = K1 = i - 1;
if (MLLLVERBOSE)
printf("BASIS0 completed updating the basis\n");
/* Here K1 = no. of LI rows in *B1ptr found by Gram Schmidt process.
flag = 0 means all the rho = beta rows of *B1ptr are LI;
flag = 1 means that the first rho = beta - 1 rows of *B1ptr are
LI, but the
beta-th row is a LC of the preceding rows. So beta = number of
rows of *B1ptr
currently being examined by the LLL algorithm. */
k = tau;
M1 = CHANGE(m1);
N1 = CHANGE(n1);
while (k <= beta)
{
if (MLLLVERBOSE)
printf("beta - k = %u\n", beta -k);
l = k - 1;
Flag = STEP40(k, l, &L, &B1ptr, D);
if (k >= norig && GCDFLAG)
{
gcdflag = 1;
goto FOUND;
}
if (Flag)/* STEP 9 of POHST. */
{
sigma++;
if (MLLLVERBOSE)
printf("relation vector number %u
found\n", sigma);
tau = k++;
goto found;
}
X = MULTI(D[k - 2], D[k]);
Y = MULTI(D[k - 1], D[k - 1]);
Tmp = Y;
Y = MULT_I(Y, m1);
FREEMPI(Tmp);
R = MULTI(L->V[k - 1][k - 2], L->V[k - 1][k - 2]);
Z = ADD0I(X, R);
Tmp = Z;
Z = MULT_I(Z, n1);
FREEMPI(Tmp);
if (RSV(Y, Z) == 1)/*& STEP 5 of POHST. */
{
flagg = 0;
if (EQZEROI(D[k]) && EQZEROI(R))
{/* CASE B=0 of STEP 7 of POHST. */
FREEMPI(D[k - 1]);
D[k - 1] = ZEROI();
STEP80(k, &B1ptr, &L);
if (k - 1 < K1)
K1 = k - 1;
/* The swap may have changed 2nd last
row */
/* of *B1ptr. */
for (t = 0; t < m; t++)
{
FREEMPI(C->V[k - 2][t]);
C->V[k - 2][t] = ZEROI();
}
beta--;
flagg = 1;
if (k > 2)
k--;
FREEMPI(X);
FREEMPI(Y);
FREEMPI(R);
FREEMPI(Z);
continue;
}
if (flagg == 0)
{
for (i = k + 1; i <= beta; i++)
STEP7(i, k, &L, D);
}
STEP80(k, &B1ptr, &L);
if (k - 2 < K1)
K1 = k - 2;
/* swap will change last two rows of *B1ptr.
*/
FREEMPI(R);
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
{ /* STEP 6 of POHST. */
FREEMPI(R);
FREEMPI(X);
FREEMPI(Y);
FOUND:
for (l = k - 2; l >= 1; l--)
{
Flag = STEP40(k, l, &L, &B1ptr, D);
if (Flag)
{
FREEMPI(Z);
sigma++;
if (MLLLVERBOSE)
printf("relation
vector number %u found\n", sigma);
tau = k++;
goto found; /* STEP 9 of
POHST. */
}
}
k++;
}
if (!gcdflag)
FREEMPI(Z);
else
break;
}
FREEMPI(M1);
FREEMPI(N1);
FREEMATI(C);
/*
printf("L = \n");
PRINTMATI(0,L->R-1,0,L->C-1,L);
for (i = 0; i <= norig; i++)
{
printf("D[%u] = ", i);PRINTI(D[i]);printf(", ");
}
printf("\n");
*/
FREEMATI(L);
for (i = 0; i <= norig; i++)
FREEMPI(D[i]);
ffree((char *)D, (1 + norig) * sizeof(MPI *));
if (MLLLVERBOSE)
{
printf("number of basis vectors found = %u ;\n", rho);
printf("number of relation vectors found = %u .\n",
sigma);
}
return (B1ptr);
}

unsigned int STEP40(k, l, Lptr, Bptr, D)
/*
* updates *Lptr, *Bptr.
* returns 1 if row k of *Bptr becomes zero, returns zero otherwise.
*/
unsigned int k, l;
MPI *D[];
MPMATI **Lptr, **Bptr;
{
unsigned int j, flag = 1, t, m;
MPI *X, *Y, *R, *Tmp;
MPMATI *TmpMATI;

m = (*Bptr)->C;
Y = MULT_I((*Lptr)->V[k - 1][l - 1], 2);
if (RSV(Y, D[l]) == 1)
{
R = NEAREST_INTI((*Lptr)->V[k - 1][l - 1], D[l]);
X = MINUSI(R);
TmpMATI = *Bptr;
*Bptr = ADD_MULT_ROWI(l - 1, k - 1, X, *Bptr);
if (MLLLVERBOSE)
{
printf("Row %u -> Row %u + ", k,k);PRINTI
(X);printf(" x Row %u\n", l);
PRINTMATI(0,(*Bptr)->R-1,0,(*Bptr)->C-
1,*Bptr);
GetReturn();
}
FREEMATI(TmpMATI);
/*
MAXI = UPDATEMAXI(MAXI, *Bptr);
*/
FREEMPI(X);
for (j = 1; j < l; j++)
{
X = MULTI((*Lptr)->V[l - 1][j - 1], R);
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], R);
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(R);
}
for (t = 0; t < m; t++)
{
if (!EQZEROI((*Bptr)->V[k - 1][t]))
{
flag = 0;
break;
}
}
if (flag)
{
TmpMATI = *Bptr;
*Bptr = DELETE_ROWI(k, *Bptr);
FREEMATI(TmpMATI);
}
FREEMPI(Y);
return (flag);
}

unsigned int STEP400(k, l, Lptr, Bptr, D, gcdflag, norig)
/*
* updates *Lptr, *Bptr.
* returns 1 if row k of *Bptr becomes zero, returns zero otherwise.
*/
unsigned int k, l, gcdflag, norig;
MPI *D[];
MPMATI **Lptr, **Bptr;
{
unsigned int j, flag = 1, t, m;
MPI *X, *Y, *R, *Tmp;
MPMATI *TmpMATI;

m = (*Bptr)->C;
Y = MULT_I((*Lptr)->V[k - 1][l - 1], 2);
if (RSV(Y, D[l]) == 1)
{
R = NEAREST_INTI((*Lptr)->V[k - 1][l - 1], D[l]);
X = MINUSI(R);
TmpMATI = *Bptr;
*Bptr = ADD_MULT_ROWI(l - 1, k - 1, X, *Bptr);
FREEMATI(TmpMATI);
/*
MAXI = UPDATEMAXI(MAXI, *Bptr);
*/
FREEMPI(X);
for (j = 1; j < l; j++)
{
X = MULTI((*Lptr)->V[l - 1][j - 1], R);
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], R);
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(R);
}
for (t = 0; t < m; t++)
{
if (!EQZEROI((*Bptr)->V[k - 1][t]))
{
flag = 0;
break;
}
}
if (flag)
{
TmpMATI = *Bptr;
*Bptr = DELETE_ROWI(k, *Bptr);
FREEMATI(TmpMATI);
}
FREEMPI(Y);
return (flag);
}

void STEP80(k, B1ptr, Lptr)
MPMATI **B1ptr, **Lptr;
unsigned int k;
{
MPI *T;

*B1ptr = SWAP_ROWSI1(k - 2, k - 1, *B1ptr);
if (MLLLVERBOSE)
{
printf("Swapping Rows %u and %u\n", k-1,k);
PRINTMATI(0,(*B1ptr)->R-1,0,(*B1ptr)->C-1,*B1ptr);
GetReturn();
}
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();
return;
}

MPMATI *EXTGCD(MPMATI *Dptr, MPI **Aptr, MPMATI **Q, USI m1, USI n1)
/*
* Input: an n x 1 vector of MPI's.
* Output: *Aptr = gcd of the vector of MPI's. Also we return a small
set of
* multipliers using the LLL method of Havas and Matthews.
* parameters m1 and n1 were put in at the request of George Havas on
11/7/94.
* Normally m1/n1 = 3/4.
* matrix Q of the LLL extended gcd paper of Havas-Matthews is
returned.
*/
{
MPMATI *Temp, *P, *SM, *R;
USI i, m, n, *alpha, nz;
MPI **Tptr;

n = Dptr->R;
m = n - 1;
alpha = KB_ROWP(Dptr, &R, &nz);
Temp = HERMITE1P(Dptr, R, &P, alpha, nz);
FREEMATI(R);
ffree((char *)alpha, (Dptr->C) * sizeof(USI));
*Aptr = COPYI(Temp->V[0][0]);
FREEMATI(Temp);
Tptr = P->V[0];
for (i = 0; i < m; i++)
P->V[i] = P->V[i+1];
P->V[n - 1] = Tptr;
GCDFLAG = 1;
*Q = BASIS_REDUCTION0(P, m1, n1);
SM = BUILDMATI(1, n);
for (i = 0; i < n; i++)
SM->V[0][i] = COPYI((*Q)->V[m][i]);
FREEMATI(P);
GCDFLAG = 0;
return (SM);
}

MPI *LENGTHSQRI(MPMATI *Mptr, USI i)
/*
* Returns the square of the length of row i of matrix *Mptr.
*/
{
MPI *SUM, *T, *T1;
USI j;

SUM = ZEROI();
for (j = 0; j < Mptr->C; j++)
{
T = SUM;
T1 = MULTI(Mptr->V[i][j], Mptr->V[i][j]);
SUM = ADD0I(SUM, T1);
FREEMPI(T);
FREEMPI(T1);
}
return (SUM);
}

MPI *LENGTHSQCI(MPMATI *Mptr, USI j)
/*
* Returns the square of the length of column j of matrix *Mptr.
*/
{
MPI *SUM, *T, *T1;
USI i;

SUM = ZEROI();
for (i = 0; i < Mptr->R; i++)
{
T = SUM;
T1 = MULTI(Mptr->V[i][j], Mptr->V[i][j]);
SUM = ADD0I(SUM, T1);
FREEMPI(T);
FREEMPI(T1);
}
return (SUM);
}

MPMATI *BASIS_REDUCTION00(MPMATI *Bptr, USI m1, USI n1, USI norig)
/*
* Input: *Bptr, a matrix of MPI's, whose first row is not zero.
* Output: a pointer to an MPMATI whose rows form a reduced basis for
* the lattice spanned by the rows of *Bptr. This basis is reduced in
the
* sense of the paper "Factoring polynomials with rational
coefficients" by
* A. K. Lenstra, H. W. Lenstra and L. Lovasz, Math. Ann. 261, 515-
534 (1982)
* using the modified version in "Solving exponential Diophantine
equations
* using lattice basis reduction algorithms" by B. M. M. De Weger, J.
No. Theory
* 26, 325-367 (1987). No change of basis matrix is returned.
* De Weger's algorithm has been changed to cater for arbitrary
matrices. The
* the rows are now in general linearly dependent.
* We use the fact that the Gram Schmidt process detects the first
row
* which is a linear combination of the preceding rows. We employ a
modification
* of the LLL algorithm outlined by M. Pohst in J. Symbolic
Computation (1987)4,
* 123-127. We call this the MLLL algorithm.
* If we are using this algorithm to find small multipliers for the
extended
* gcd problem, GCDFLAG is set in IMPROVEP() and gcdflag is set below.
* m1 / n1 is usually taken to be 3 / 4.
* For use in IMPROVEP().
*/
{
unsigned int i, k, n, m, t, flag = 0, Flag = 0, gcdflag = 0;
unsigned int flagg, beta, K1 = 0, tau = 2, sigma = 0, rho,
norigg;
unsigned int K, j, REPEAT, iterate;
int l, tt;
MPI **D, *X, *Y, *Z, *H, *Tmp, *R, *M1, *N1;
MPI *t1, *t2, *t3, *t4;
MPMATI *C, *L, *B1ptr, *temp1;

Z = NULL;
m = Bptr->C;
norigg = n = Bptr->R;
B1ptr = COPYMATI(Bptr);
D = (MPI **)mmalloc((1 + n) * sizeof(MPI *));
D[0] = ONEI();
for (i = 1; i <= n; i++)
D[i] = ZEROI();
C = ZEROMNI(n, m);
L = ZEROMNI(n, n);

found:
n = B1ptr->R;
i = (K1 == 0) ? 1 : K1;
/* K1 = no. of consecutive rows of *B1ptr that don't need
updating
for the Gram Schmidt process. */
while (i <= B1ptr->R)
{
BASIS_UPDATE(i, m, &C, &L, B1ptr, D);
flag = 1;
for (t = 0; t < m; t++)
{
if (!EQZEROI(C->V[i - 1][t]))
{
flag = 0;
break;
}
}
if (flag)
break;
X = ZEROI();
for (t = 0; t < m; t++)
{
H = MULTI(C->V[i - 1][t], C->V[i - 1][t]);
Tmp = X;
X = ADDI(X, H);
FREEMPI(Tmp);
FREEMPI(H);
}
FREEMPI(D[i]);
D[i] = INT(X, D[i - 1]);
FREEMPI(X);
i++;
if (MLLLVERBOSE)
printf("i = %u\n", i);
}
/*
strcpy(buff, "L.out");
outfile = fopen(buff, "w");
FPRINTMATI(outfile,0,L->R-1,0,L->C-1,L);
fclose(outfile);
for (t = 0;t < B1ptr->R; t++){
printf("D[%u] = ", t);PRINTI(D[t]);printf("\n");
}
GetReturn();
*/
beta = (flag) ? i : i - 1;
rho = K1 = i - 1;
if (MLLLVERBOSE)
printf("completed updating the basis\n");
/* Here K1 = no. of LI rows in *B1ptr found by Gram Schmidt process.
flag = 0 means all the rho = beta rows of *B1ptr are LI;
flag = 1 means that the first rho = beta - 1 rows of *B1ptr are
LI, but the
beta-th row is a LC of the preceding rows. So beta = number of
rows of *B1ptr
currently being examined by the LLL algorithm. */
k = tau;
M1 = CHANGE(m1);
N1 = CHANGE(n1);
while (k <= beta)
{
if (MLLLVERBOSE)
printf("beta - k = %u\n", beta -k);
if (gcdflag)
l = norig;
else
l = k - 1;
Flag = STEP40(k, l, &L, &B1ptr, D);
if (k > norig && GCDFLAG)
{
gcdflag = 1;
if (MLLLVERBOSE)
printf("improving row %u\n", k);
goto FOUND;
}
if (Flag)/* STEP 9 of POHST. */
{
sigma++;
if (MLLLVERBOSE)
printf("relation vector number %u
found\n", sigma);
tau = k++;
goto found;
}
X = MULTI(D[k - 2], D[k]);
Y = MULTI(D[k - 1], D[k - 1]);
Tmp = Y;
Y = MULT_I(Y, m1);
FREEMPI(Tmp);
R = MULTI(L->V[k - 1][k - 2], L->V[k - 1][k - 2]);
Z = ADD0I(X, R);
Tmp = Z;
Z = MULT_I(Z, n1);
FREEMPI(Tmp);
if (RSV(Y, Z) == 1)/*& STEP 5 of POHST. */
{
flagg = 0;
if (EQZEROI(D[k]) && EQZEROI(R))
{/* CASE B=0 of STEP 7 of POHST. */
FREEMPI(D[k - 1]);
D[k - 1] = ZEROI();
STEP80(k, &B1ptr, &L);
if (k - 1 < K1)
K1 = k - 1;
/* The swap may have changed 2nd last
row */
/* of *B1ptr. */
for (t = 0; t < m; t++)
{
FREEMPI(C->V[k - 2][t]);
C->V[k - 2][t] = ZEROI();
}
beta--;
flagg = 1;
if (k > 2)
k--;
FREEMPI(X);
FREEMPI(Y);
FREEMPI(R);
FREEMPI(Z);
continue;
}
if (flagg == 0)
{
for (i = k + 1; i <= beta; i++)
STEP7(i, k, &L, D);
}
STEP80(k, &B1ptr, &L);
if (k - 2 < K1)
K1 = k - 2;
/* swap will change last two rows of *B1ptr.
*/
FREEMPI(R);
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
{ /* STEP 6 of POHST. */
FREEMPI(R);
FREEMPI(X);
FREEMPI(Y);
FOUND:
if (gcdflag)
{
if (MLLLVERBOSE)
printf("gcdflag set\n");
K = norig + 1;
}
else
K = k;
for (l = K - 2; l >= 1; l--)
{
Flag = STEP40(k, l, &L, &B1ptr, D);
if (Flag)
{
FREEMPI(Z);
sigma++;
if (MLLLVERBOSE)
printf("relation
vector number %u found\n", sigma);
tau = k++;
goto found; /* STEP 9 of
POHST. */
}
}
k++;
}
if (!gcdflag)
FREEMPI(Z);
}
FREEMPI(M1);
FREEMPI(N1);
FREEMATI(C);
FREEMATI(L);
for (i = 0; i <= norigg; i++)
FREEMPI(D[i]);
ffree((char *)D, (1 + norigg) * sizeof(MPI *));
/* Now to improve the last n - norig rows of *B1ptr, using
Gauss
lattice reduction - pointed out by George Havas in Sims' book -
8/11/94: */

/*
strcpy(buff, "progress");
outfile = fopen(buff, "a+");
FPRINTMATI(outfile,0,B1ptr->R-1,0,B1ptr->C-1,B1ptr);
*/

REPEAT = 1;
temp1 = BUILDMATI(1, B1ptr->C);
for (k = 0; k < B1ptr->C; k++)
temp1->V[0][k] = ZEROI();
iterate = 0;
QSORTMATI(B1ptr, 0, norig - 1);
while (REPEAT)
{
REPEAT = 0;
for (j = 1; j < norig; j++)
{
t4 = DOTRI(B1ptr, j, j);
for (i = 0; i < j; i++)
{
t1 = DOTRI(B1ptr, j, i);
t2 = DOTRI(B1ptr, i, i);
t3 = NEAREST_INTI(t1, t2);
FREEMPI(t1);
FREEMPI(t2);
if (t3->S == 0)
{
FREEMPI(t3);
continue;
}
t1 = t3;
t3 = MINUSI(t3);
FREEMPI(t1);
for (k = 0; k < B1ptr->C; k++)
{
FREEMPI(temp1->V[0][k]);
X = MULTI(t3, B1ptr->V[i][k]);
temp1->V[0][k] = ADDI(X,
B1ptr->V[j][k]);
FREEMPI(X);
}
t1 = DOTRI(temp1, 0, 0);
tt = RSV(t1, t4);
tt = RSV(t1, t4);
if (tt == -1)
{
ADD_MULT_ROWI0(i, j, t3,
B1ptr);
if (MLLLVERBOSE)
{
printf("Gauss-
improving nullspace row %u: length squared was = ", j + 1);
PRINTI(t4);
printf("\n");
printf("lengthsquared
now = ");
PRINTI(t1);
printf("\n");
}
/*
fprintf(outfile, "Gauss-
improving nullspace row %u: length squared was = ", j + 1);
FPRINTI(outfile, t4);
fprintf(outfile, "\n");
fprintf
(outfile, "lengthsquared now = ");
FPRINTI(outfile, t1);
fprintf(outfile, "\n");
fflush(outfile);
*/
FREEMPI(t4);
t4 = t1;
REPEAT = 1;
}
else
FREEMPI(t1);
FREEMPI(t3);
}
FREEMPI(t4);
}
iterate++;
QSORTMATI(B1ptr, 0, norig - 1);
}
REPEAT = 1;
iterate = 0;
while (REPEAT)
{
REPEAT = 0;
for (j = norig; j < B1ptr->R; j++)
{
t4 = DOTRI(B1ptr, j, j);
for (i = 0; i < norig; i++)
{
t1 = DOTRI(B1ptr, j, i);
t2 = DOTRI(B1ptr, i, i);
t3 = NEAREST_INTI(t1, t2);
FREEMPI(t1);
FREEMPI(t2);
if (t3->S == 0)
{
FREEMPI(t3);
continue;
}
t1 = t3;
t3 = MINUSI(t3);
FREEMPI(t1);
for (k = 0; k < B1ptr->C; k++)
{
FREEMPI(temp1->V[0][k]);
X = MULTI(t3, B1ptr->V[i][k]);
temp1->V[0][k] = ADDI(X,
B1ptr->V[j][k]);
FREEMPI(X);
}
t1 = DOTRI(temp1, 0, 0);
tt = RSV(t1, t4);
tt = RSV(t1, t4);
if (tt == -1)
{
ADD_MULT_ROWI0(i, j, t3,
B1ptr);
if (MLLLVERBOSE)
{
printf("Gauss-
improving row %u: length squared was = ", j + 1);
PRINTI(t4);
printf("\n");
printf("lengthsquared
now = ");
PRINTI(t1);
printf("\n");
}
/*
fprintf(outfile, "Gauss-
improving row %u: length squared was = ", j + 1);
FPRINTI(outfile, t4);
fprintf(outfile, "\n");
fprintf
(outfile, "lengthsquared now = ");
FPRINTI(outfile, t1);
fprintf(outfile, "\n");
fflush(outfile);
*/
FREEMPI(t4);
t4 = t1;
REPEAT = 1;
}
else
FREEMPI(t1);
FREEMPI(t3);
}
FREEMPI(t4);
}
iterate++;
}
FREEMATI(temp1);
/*
fclose(outfile);
*/
return (B1ptr);
}

MPMATI *BASIS_REDUCTION000(MPMATI *Bptr, USI m1, USI n1, MPI *N)
/*
* Input: *Bptr, a matrix of MPI's, whose first row is not zero.
* Output: a pointer to an MPMATI whose rows form a reduced basis for
* the lattice spanned by the rows of *Bptr. This basis is reduced in
the
* sense of the paper "Factoring polynomials with rational
coefficients" by
* A. K. Lenstra, H. W. Lenstra and L. Lovasz, Math. Ann. 261, 515-
534 (1982)
* using the modified version in "Solving exponential Diophantine
equations
* using lattice basis reduction algorithms" by B. M. M. De Weger, J.
No. Theory
* 26, 325-367 (1987). No change of basis matrix is returned.
* De Weger's algorithm has been changed to cater for arbitrary
matrices. The
* the rows are now in general linearly dependent.
* We use the fact that the Gram Schmidt process detects the first
row
* which is a linear combination of the preceding rows. We employ a
modification
* of the LLL algorithm outlined by M. Pohst in J. Symbolic
Computation (1987)4,
* 123-127. We call this the MLLL algorithm.
* If we are using this algorithm to find small multipliers for the
extended
* gcd problem, GCDFLAG is set in EXTGCD() and gcdflag is set below.
* m1 / n1 is usually taken to be 3 / 4.
*/
{
unsigned int i, k, l, n, m, t, flag = 0, Flag = 0, gcdflag =
0;
unsigned int flagg, beta, K1 = 0, tau = 2, sigma = 0, rho;
MPI **D, *X, *Y, *Z, *H, *Tmp, *R, *M1, *N1;
MPMATI *C, *L, *B1ptr;
unsigned int norig;

Z = NULL;
m = Bptr->C;
n = Bptr->R;
norig = n;
B1ptr = COPYMATI(Bptr);
D = (MPI **)mmalloc((1 + n) * sizeof(MPI *));
D[0] = ONEI();
for (i = 1; i <= n; i++)
D[i] = ZEROI();
C = ZEROMNI(n, m);
L = ZEROMNI(n, n);

found:
n = B1ptr->R;
i = (K1 == 0) ? 1 : K1;
/* K1 = no. of consecutive rows of *B1ptr that don't need
updating
for the Gram Schmidt process. */
while (i <= B1ptr->R)
{
BASIS_UPDATE(i, m, &C, &L, B1ptr, D);
flag = 1;
for (t = 0; t < m; t++)
{
if (!EQZEROI(C->V[i - 1][t]))
{
flag = 0;
break;
}
}
if (flag)
break;
X = ZEROI();
for (t = 0; t < m; t++)
{
H = MULTI(C->V[i - 1][t], C->V[i - 1][t]);
Tmp = X;
X = ADDI(X, H);
FREEMPI(Tmp);
FREEMPI(H);
}
FREEMPI(D[i]);
D[i] = INT(X, D[i - 1]);
FREEMPI(X);
i++;
if (MLLLVERBOSE)
printf("i = %u\n", i);
}
beta = (flag) ? i : i - 1;
rho = K1 = i - 1;
if (MLLLVERBOSE)
printf("completed updating the basis\n");
/* Here K1 = no. of LI rows in *B1ptr found by Gram Schmidt process.
flag = 0 means all the rho = beta rows of *B1ptr are LI;
flag = 1 means that the first rho = beta - 1 rows of *B1ptr are
LI, but the
beta-th row is a LC of the preceding rows. So beta = number of
rows of *B1ptr
currently being examined by the LLL algorithm. */
k = tau;
M1 = CHANGE(m1);
N1 = CHANGE(n1);
while (k <= beta)
{
if (MLLLVERBOSE)
printf("beta - k = %u\n", beta -k);
l = k - 1;
Flag = STEP4000(k, l, &L, &B1ptr, D, N);
if (k >= norig && GCDFLAG)
{
gcdflag = 1;
goto FOUND;
}
if (Flag)/* STEP 9 of POHST. */
{
sigma++;
if (MLLLVERBOSE)
printf("relation vector number %u
found\n", sigma);
tau = k++;
goto found;
}
X = MULTI(D[k - 2], D[k]);
Y = MULTI(D[k - 1], D[k - 1]);
Tmp = Y;
Y = MULT_I(Y, m1);
FREEMPI(Tmp);
R = MULTI(L->V[k - 1][k - 2], L->V[k - 1][k - 2]);
Z = ADD0I(X, R);
Tmp = Z;
Z = MULT_I(Z, n1);
FREEMPI(Tmp);
if (RSV(Y, Z) == 1)/*& STEP 5 of POHST. */
{
flagg = 0;
if (EQZEROI(D[k]) && EQZEROI(R))
{/* CASE B=0 of STEP 7 of POHST. */
FREEMPI(D[k - 1]);
D[k - 1] = ZEROI();
STEP8000(k, &B1ptr, &L, N);
if (k - 1 < K1)
K1 = k - 1;
/* The swap may have changed 2nd last
row */
/* of *B1ptr. */
for (t = 0; t < m; t++)
{
FREEMPI(C->V[k - 2][t]);
C->V[k - 2][t] = ZEROI();
}
beta--;
flagg = 1;
if (k > 2)
k--;
FREEMPI(X);
FREEMPI(Y);
FREEMPI(R);
FREEMPI(Z);
continue;
}
if (flagg == 0)
{
for (i = k + 1; i <= beta; i++)
STEP7(i, k, &L, D);
}
STEP8000(k, &B1ptr, &L, N);
if (k - 2 < K1)
K1 = k - 2;
/* swap will change last two rows of *B1ptr.
*/
FREEMPI(R);
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
{ /* STEP 6 of POHST. */
FREEMPI(R);
FREEMPI(X);
FREEMPI(Y);
FOUND:
for (l = k - 2; l >= 1; l--)
{
Flag = STEP4000(k, l, &L, &B1ptr, D,
N);
if (Flag)
{
FREEMPI(Z);
sigma++;
if (MLLLVERBOSE)
printf("relation
vector number %u found\n", sigma);
tau = k++;
goto found; /* STEP 9 of
POHST. */
}
}
k++;
}
if (!gcdflag)
FREEMPI(Z);
else
break;
}
FREEMPI(M1);
FREEMPI(N1);
FREEMATI(C);
FREEMATI(L);
for (i = 0; i <= norig; i++)
FREEMPI(D[i]);
ffree((char *)D, (1 + norig) * sizeof(MPI *));
if (MLLLVERBOSE)
{
printf("number of basis vectors found = %u ;\n", rho);
printf("number of relation vectors found = %u .\n",
sigma);
}
return (B1ptr);
}

unsigned int STEP4000(USI k, USI l, MPMATI **Lptr, MPMATI **Bptr, MPI
*D[], MPI *N)
/*
* updates *Lptr, *Bptr.
* returns 1 if row k of *Bptr becomes zero, returns zero otherwise.
*/
{
unsigned int i, j, flag = 1, t, m, n;
MPI *X, *Y, *R, *Tmp, *Temp;
MPMATI *TmpMATI;

m = (*Bptr)->C;
Y = MULT_I((*Lptr)->V[k - 1][l - 1], 2);
if (RSV(Y, D[l]) == 1)
{
R = NEAREST_INTI((*Lptr)->V[k - 1][l - 1], D[l]);
X = MINUSI(R);
TmpMATI = *Bptr;
*Bptr = ADD_MULT_ROWI(l - 1, k - 1, X, *Bptr);
FREEMATI(TmpMATI);
n = (*Bptr)->R;
for (i = 0; i < n; i++)
{
for (j = n; j < m; j++)
{
Temp = POWERI(N, m - j);
(*Bptr)->V[i][j] = INTI((*Bptr)->V[i][j],
Temp);
FREEMPI(Temp);
}
}
if (HERMITEVERBOSE)
{
printf("Row %u -> Row %u + ", k,k);PRINTI(X);printf("
x Row %u\n", l);
PRINTMATI(0,(*Bptr)->R-1,0,(*Bptr)->C-1,*Bptr);
GetReturn();
}
for (i = 0; i < n; i++)
{
for (j = n; j < m; j++)
{
Temp = POWERI(N, m - j);
(*Bptr)->V[i][j] = MULTI((*Bptr)->V[i][j],
Temp);
FREEMPI(Temp);
}
}
/*
MAXI = UPDATEMAXI(MAXI, *Bptr);
*/
FREEMPI(X);
for (j = 1; j < l; j++)
{
X = MULTI((*Lptr)->V[l - 1][j - 1], R);
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], R);
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(R);
}
for (t = 0; t < m; t++)
{
if (!EQZEROI((*Bptr)->V[k - 1][t]))
{
flag = 0;
break;
}
}
if (flag)
{
TmpMATI = *Bptr;
*Bptr = DELETE_ROWI(k, *Bptr);
FREEMATI(TmpMATI);
}
FREEMPI(Y);
return (flag);
}
void STEP8000(USI k, MPMATI **B1ptr, MPMATI **Lptr, MPI *N)
{
MPI *T, *Temp;
USI i, j, n, m;

*B1ptr = SWAP_ROWSI1(k - 2, k - 1, *B1ptr);
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();

n = (*B1ptr)->R;
m = (*B1ptr)->C;
for (i = 0; i < n; i++)
{
for (j = n; j < m; j++)
{
Temp = POWERI(N, m - j);
(*B1ptr)->V[i][j] = INTI((*B1ptr)->V[i][j],
Temp);
FREEMPI(Temp);
}
}
if (HERMITEVERBOSE)
{
printf("Swapping Rows %u and %u\n", k-1,k);
PRINTMATI(0,(*B1ptr)->R-1,0,(*B1ptr)->C-1,*B1ptr);
GetReturn();
}
for (i = 0; i < n; i++)
{
for (j = n; j < m; j++)
{
Temp = POWERI(N, m - j);
(*B1ptr)->V[i][j] = MULTI((*B1ptr)->V[i][j],
Temp);
FREEMPI(Temp);
}
}

return;
}