A GCD based Hermite Normal Form algorithm, stabilized in the Chou Collins  or Lenstra Lenstra Lovász manner.

We present a rather simple code to illustrate that the principles that prevent coefficient explosion in the LLL algorithm also suffice in a naive GCD based Hermite Normal Form algorithm.This should not come as a surprise, because Chou and Collins have already shown this long ago, at least for an input matrix of full rank. 
Query: Does anyone know of a more recent treatment, not restricted to input of full rank?
A GCD based operation amounts to a sequence of ordinary row operations. We know from LLL that one does not need to length reduce everything all the time. From my analysis of an extended LLL algorithm it is clear that it suffices to do the following. Each time before a new row is allowed to enter the fray one must bring the matrix above it into (lower triangular) Hermite Normal Form. So one should not be satisfied with a lower triangular form.  To avoid coefficient explosion one should go all the way to Hermite Normal Form. One can prove an (extended) LLL like estimate for the given code in the manner of my analysis of extendedlll. Thus the analysis of Chou and Collins is very much related to that of LLL and extendedlll. In fact Chou and Collins covered the case that corresponds with the extended hnf call hnf[matrix,symbol] as introduced below. They did not cover the case of the call hnf[matrix], which is faster and allows dependent vectors. (They added an identity matrix to make the vectors independent).  Our main point is that it helps to think of their Hermite Normal Form algorithm as belonging to the same family as the LLL algorithm, or vice versa. 
We also like to stress that our simple code describes an algorithm that often outperforms the algorithm of P. Domich, R. Kannan and L. Trotter, as implemented by Victor Shoup in NTL. Recall that their algorithm counters coefficient explosion by using modular arithmetic. When no small multiple of the pivots exists (or when none is known), it seems much better to prevent the explosion than to try and cure it by computing modulo some large number.
Recall that one sees terrible explosions when the intermediate stabilizing steps are left out. 
Wilberd van der Kallen ,   November 2003 but converted to Mathematica 7 and a sign corrected in 2017.

The hnf code

SetAttributes[hnf,HoldAll]

hnf::usage:="{newmatrix,rank}=hnf[mat] gives a lower triangular Hermite Normal Form of the integral matrix mat, together with its rank.\n
{newmatrix,rank}=hnf[matrix,symbol] also assigns to symbol the matrix trans so that trans.mat=newmatrix."

hnf[bb_,uu_:False]:=
Module[
{u,m=Length[bb],n=Length[bb[[1]]],i,j,l,pp,bb0,rank=0,dimker=0,k=1,maxk=0,ans}
,
pp=0 Range[m];(* initialize as a row of zeroes *)
If[uu=!=False
,(* If the transformation is requested, 
put an identity matrix in front, to keep track of the transformation. *)
bb0=Transpose[Join[IdentityMatrix[m],Transpose[bb]]];
n=n+m
,
bb0=bb
];
For[i=1,i≤m,i++
,(* Record how far row i sticks out. *)
j=n;
While[(j>0)&&(bb0[[i,j]]==0),j--];
pp[[i]]=j;
If[(j>0)&&(bb0[[i,j]]<0),bb0[[i]]=-bb0[[i]]];(* Keep diagonal positive. *)
];
While[k≤m
,
maxk=k;(* number of new row that is to be processed *)
While[(k>dimker+1)&&(pp[[k]]≤pp[[k-1]])
,
swap[k,bb0,pp];(* uses extended Euclid *)
k--;
];
If[pp[[k]]==0,dimker++,rank++];
For[Null,k≤maxk,k++
,(* Reduce the columns, as required by Hermite. *)
For[l=k-1,l>dimker,l--
,(* We use the order of reduction that is recommended by LLL or Chou Collins *)
j=pp[[l]];
u=Quotient[bb0[[k,j]],bb0[[l,j]]];
bb0[[k]]=bb0[[k]]-u bb0[[l]]
]
];
k=maxk+1 (* Proceed to next row, now that part above it is in HNF form. *)
];(* end of main loop *)
If[uu=!=False
,(* Extract new matrix and transformation. *)
n=n-m;
uu=Take[#,m]&/@bb0;
ans=Take[#,-n]&/@bb0;
rank=0;
For[i=1,i≤m,i++,If[pp[[i]]>m,rank++]]
,
ans=bb0
];
{ans,rank}
]

SetAttributes[swap,HoldAll];

swap[kk_,bb_,pp_]:=(* Swaps row k with row k-1, using extended GCD if needed. *) 
Module[
{i,j,l,t1,t2,e,x,y,k=kk}
,
If[pp[[k]]<pp[[k-1]]
,(* Row k-1 sticks out further. *)
{bb[[k-1]],bb[[k]]}={bb[[k]],bb[[k-1]]};
{pp[[k-1]],pp[[k]]}={pp[[k]],pp[[k-1]]};
,(* They stick out equally. *)
j=pp[[k]];
{e,{x,y}}=ExtendedGCD[bb[[k-1,j]],bb[[k,j]]];
t1=Quotient[bb[[k,j]],e];
t2=-Quotient[bb[[k-1,j]],e];
{bb[[k-1]],bb[[k]]}={{t1,t2},{x,y}}.{bb[[k-1]],bb[[k]]};
For[i=0,True,i++,
If[(j>0)&&(bb[[k,j]]<0),bb[[k]]=-bb[[k]]];(* Keep diagonal positive. *)
If[i>0,Break[]];
k--;(* Look again how far it sticks out. *)
While[(j>0)&&(bb[[k,j]]==0),j--];
pp[[k]]=j
]
]
]

Examples

?hnf

ma={{990519,-4051508,-3802647,-43245,-369337,1967100},{1598245,-8116013,-6067138,-3323116,-2053036,3161362},{1786244,1069431,-465693,-227421,-2029817,-542818},{1084810,1537944,260218,-1057178,-1016978,-362744},{824048,-2096515,-2199740,-427760,-447464,1118810},
{-90848,5193882,1635132,6421956,2423012,-1147532}};

{new,rank}=hnf[ma,u];

u.ma==new

hnf[ma][[1]]==new

MatrixForm[ma]

MatrixForm[new]

MatrixForm[u]

The above matrix "new" is in lower triangular HNF form. A matrix is in lower triangular HNF form if it is either zero, or there is a "first pivot", in the example the entry 15915414528, which
is positive
has zero rows above it
has zeroes to the right of it in its row
has nonegative entries in its column, all other entries being smaller
has to its right a matrix in lower triangular HNF form. In the example this submatrix is:

MatrixForm[Take[#,-4]&/@new]

If you prefer a more traditional upper triangular (echelon) Hermite Normal Form, then you may proceed like this:

{new,rank}=hnf[Reverse[Reverse/@ma],u];new=Reverse[Reverse/@new];u=Reverse[Reverse/@u];u.ma==new

MatrixForm[ma]

MatrixForm[new]

MatrixForm[u]

Let us try a few random cases.

mat[n_,m_,d_]:=Table[RandomInteger[{-d,d}],{i,m},{j,n}]

First an example where the normal form itself is not so favorable. Then of course the algorithm can't keep the integers smaller than that.

{new,rank}=hnf[ma=mat[60,40,1000],u];

Max[Abs[new]]

Max[Abs[u]]

Next an example where the normal form itself is fine.

{new,rank}=hnf[ma=mat[40,60,1000],u];

Max[Abs[new]]

Max[Abs[u]]

You see that the entries of the transformation matrix get as bad as in the previous case. One has an estimate for the bitlength of entries of the transformation matrix that is roughly the number of rows of the input matrix times the the bitlength of its largest entry .

Log[N[Max[u]]]

Log[N[Max[ma]]]

Length[ma] Log[Length[ma]N[Max[ma]]]

To see how bad it gets when you are stupid and wait with the reductions until the end, let us consider the variant where one does not do any of these reductions.

SetAttributes[hnfstupid,HoldAll]

hnfstupid::usage:="{newmatrix,rank}=hnf[mat]  gives a lower triangular Normal Form of the integral matrix mat, together with its rank. 
{newmatrix,rank}=hnf[matrix,symbol] also assigns to symbol the matrix trans so that trans.mat=newmatrix."

hnfstupid[bb_,uu_:False]:=Module[{u,m=Length[bb],n=Length[bb[[1]]],i,j,l,pp,bb0,rank=0,dimker=0,k=1,maxk=0,ans},pp=0 Range[m];
If[uu=!=False,bb0=Transpose[Join[IdentityMatrix[m],Transpose[bb]]];
n=n+m,bb0=bb];
For[i=1,i≤m,i++,j=n;
While[(j>0)&&(bb0[[i,j]]==0),j--];
pp[[i]]=j;];
While[k≤m,maxk=k;
While[(k>dimker+1)&&(pp[[k]]≤pp[[k-1]]),swap[k,bb0,pp];
k--;];
If[pp[[k]]==0,dimker++,rank++];
For[Null,k≤maxk,k++,For[l=k-1,l>dimker,l--,j=pp[[l]];
u=Quotient[bb0[[k,j]],bb0[[l,j]]];
(* bb0[[k]]=bb0[[k]]-u bb0[[l]] *)
] ];
k=maxk+1];
If[uu=!=False,n=n-m;
uu=Take[#,m]&/@bb0;
ans=Take[#,-n]&/@bb0;
rank=0;
For[i=1,i≤m,i++,If[pp[[i]]>m,rank++]],ans=bb0];
{ans,rank}]

Now let us see how bad it gets.

{new,rank}=hnfstupid[ma=mat[10,20,1000],u];

Max[Abs[new]]

Max[Abs[u]]

Compare that with

{new,rank}=hnf[ma,u];

Max[Abs[new]]

Max[Abs[u]]

We leave it to you to explore larger cases.

Better transformation matrices.

Following Havas-Majewski-Matthews we can use LLL to improve on the transformation matrix when there are zero rows in the normal form. The difference with their approach is that we postpone this improving until our lightweight hnf has run its course. We believe this often makes sense. But we did not run tests to compare the two approaches. Note that we will be mixing the compiled code LatticeReduce with uncompiled code.

The hnfLLL code

SetAttributes[hnfLLL,HoldAll];

hnfLLL::usage:="{newmatrix,rank}=hnfLLL[mat,symbol] gives a lower triangular Hermite Normal Form of the integral matrix mat, together with its rank.\n
It also assigns to symbol an optimized matrix trans so that trans.mat=newmatrix.";

hnfLLL[bb_,uu_]:=
Module[{u,m=Length[bb],i,j,k,l,dd,lam,rank,dimker,ans},
dd=Range[m+1];
lam=0 IdentityMatrix[m];
{ans,rank}=hnf[bb,uu];
dimker=m-rank;
If[dimker>0,(* otherwise there is no point calling hnfLLL *)
uu=Join[LatticeReduce[Take[uu,dimker]],Take[uu,-rank]] (* LatticeReduce the nullspace. *)
];
For[k=1,k≤dimker,k++,(* Build the GramSchmidt table. *)
For[i=1,i≤k,i++,gramS[uu,dd,lam,i,k]];
For[l=k-1,l>0,l--,reduce[k,l,dd,lam,uu]] (* This should not be necessary. *)
];
For[k=dimker+1,k≤m,k++,(* Build more of the Gram-Schmidt table. *)
For[l=1,l≤dimker,l++,gramS[uu,dd,lam,l,k]];
For[l=dimker,l>0,l--,reduce[k,l,dd,lam,uu]] (* This should help a lot. *)
];
{ans,rank}
]

SetAttributes[gramS,HoldAll];

gramS[uu_,dd_,lam_,low_,up_]:=
Module[{u,t1,t2,i},(* Compute entry of Gram-Schmidt table. *)
u=uu[[low]].uu[[up]];
For[i=1,i<low,i++,
t1=u dd[[i+1]];
t2=lam[[up,i]] lam[[low,i]];
t1=t1-t2;
u=Quotient[t1,dd[[i]]];
];
If[low<up,lam[[up,low]]=u,dd[[up+1]]=u]
]

SetAttributes[reduce,HoldAll];

reduce[k_,l_,dd_,lam_,uu_]:=
Module[{t1,t2,u,i},(* Subtract multiple of row l from row k. *)
t2=dd[[l+1]];
t1=Abs[2 lam[[k,l]]];
If[t1>t2,
u=Quotient[2lam[[k,l]]+t2,2 t2];
uu[[k]]=uu[[k]]-u uu[[l]];
lam[[k,l]]=lam[[k,l]]-t2 u;
For[i=1,i<l,i++,
lam[[k,i]]=lam[[k,i]]-u lam[[l,i]]
]
]
]

Examples

?hnfLLL

ma={{51,78,-17},{-75,-63,48},{-92,87,-66},{-96,7,73},{41,17,3},{-78,6,58},{47,-34,-71},{34,28,-11},{13,-59,-94},{90,-32,-33},{76,-62,48},{-20,16,-87},{2,92,29},{14,-89,-32}};

hnf[ma,uu];Max[Abs[uu]]

hnfLLL[ma,uuLLL];Max[Abs[uuLLL]]

uuLLL.ma==uu.ma

Det[uuLLL]

MatrixForm[ma]

MatrixForm[uu]

MatrixForm[uuLLL]

ma=mat[1,18,100]

hnf[ma,uu]

Max[Abs[uu]]

hnfLLL[ma,uuLLL]

Max[Abs[uuLLL]]