macro( HermiteNormalForm = `Matrix/HermiteNormalForm` );

#
#--> HermiteNormalForm
#
#  Compute the Hermite normal form of the n by m matrix A containing
#  integer entries.  See the comments in linalg[hermite].
#
#  Caveat: should return the multiplier matrix K such that K H = A .
#  See also: linalg[hermite], linalg[ismith] and linalg[gausselim]
#
#  MBM Dec/88  (modified for Gauss by TJS Mar 1989)
#

HermiteNormalForm := proc(B,M)

    local A, a, b, c, g, i, j, m, n, q, r, s, t, temp, ui, C;
  
    n := M['Rows'](B);
    m := M['Cols'](B);

    A := M['ArrayCoeffs'](B);

    C := M['CoefficientRing'];

    r := 1;
    for c to m while r <= n do

        #  Pivot selection from row k
        for i from r to n while C[`=`](A[i,c], C[0]) do  od;
        if i > n then next fi;

        #  Select the smallest non-zero entry as the pivot
        for j from i+1 to n do
            if C[`=`](A[j,c], C[0]) then next fi;
	    if C['SmallerEuclideanNorm'](A[j,c],A[i,c]) then i := j fi
        od;

        #  Pivot is A[i,c]: interchange row i with row r if necessary
        if i <> r then
            for j from c to m do
                t := A[i,j]; A[i,j] := A[r,j]; A[r,j] := t
            od
        fi;

        #  Ensure that the leading entry is unit normal
	ui := C['Unit'](A[r,c]);
        if C[`<>`](ui, C[1]) then
	    ui := C['Inv'](ui);
            for j from c to m do A[r,j] := C[`*`](ui,A[r,j]) od
        fi;

        #  Zero out column c from c+1 to n
        for i from r+1 to n do
            if C[`=`](A[i,c], C[0]) then next fi;
            g := C['Gcdex'](A[r,c], A[i,c], 's', 't');
            a := C['Quo'](A[r,c],g); b := C['Quo'](A[i,c],g);
            #
            #  We have  s A[r,c]/g + t A[i,c]/g = 1
            #
            #       [  s  t ]  [ A[r,c]  A[r,j] ]   [ g  ... ]
            #       [       ]  [                ] = [        ]
            #       [ -b  a ]  [ A[i,c]  A[i,j] ]   [ 0  ... ]
            #
            #       for j = c+1..m  where note  s a + t b = 1
            #
            for j from c+1 to m do
                temp := C[`+`](C[`*`](s,A[r,j]), C[`*`](t,A[i,j]));
                A[i,j] := C[`-`](C[`*`](a,A[i,j]), C[`*`](b,A[r,j]));
                A[r,j] := temp
            od;
            A[r,c] := g;
            A[i,c] := C[0];
        od;

        #  Reduce such that 0 <= A[i,c] < A[r,c] for i <> r
        for i to r-1 do
	    q := C['Quo']( A[i,c], A[r,c], 't');
            if C[`=`](q, C[0]) then next fi;
            for j from c+1 to m do 
		A[i,j] := C[`-`](A[i,j], C[`*`](q,A[r,j])) 
	    od;
            A[i,c] := t
        od;

        r := r+1

    od;

    M[Matrix](A,n,m);

end:

save `HNF.m`;
quit
