Attributes[hnfLLL] = {HoldAll}
 
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, uu = Join[LatticeReduce[Take[uu, dimker]], 
         Take[uu, -rank]]]; For[k = 1, k <= dimker, k++, 
       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]]]; 
      For[k = dimker + 1, k <= m, k++, For[l = 1, l <= dimker, l++, 
         gramS[uu, dd, lam, l, k]]; For[l = dimker, l > 0, l--, 
         reduce[k, l, dd, lam, uu]]]; {ans, rank}]
 
hnfLLL /: 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."
 
Attributes[hnf] = {HoldAll}
 
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]; 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[#1, m] & ) /@ bb0; 
        ans = (Take[#1, -n] & ) /@ bb0; rank = 0; For[i = 1, i <= m, i++, 
         If[pp[[i]] > m, rank++]], ans = bb0]; {ans, rank}]
 
hnf /: 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."
 
Attributes[swap] = {HoldAll}
 
swap[kk_, bb_, pp_] := Module[{i, j, l, t1, t2, e, x, y, k = kk}, 
     If[pp[[k]] < pp[[k - 1]], {bb[[k - 1]], bb[[k]]} = 
        {bb[[k]], bb[[k - 1]]}; {pp[[k - 1]], pp[[k]]} = 
        {pp[[k]], pp[[k - 1]]}; , 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]]}; 
       If[j > 0 && bb[[k,j]] < 0, bb[[k]] = -bb[[k]]]; k--; 
       While[j > 0 && bb[[k,j]] == 0, j--]; pp[[k]] = j; 
       If[j > 0 && bb[[k,j]] < 0, bb[[k]] = -bb[[k]]]; ]]
 
Attributes[gramS] = {HoldAll}
 
gramS[uu_, dd_, lam_, low_, up_] := Module[{u, t1, t2, i}, 
     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]]
 
Attributes[reduce] = {HoldAll}
 
reduce[k_, l_, dd_, lam_, uu_] := Module[{t1, t2, u, i}, 
     t2 = dd[[l + 1]]; t1 = Abs[2*lam[[k,l]]]; 
      If[t1 > t2, u = Quotient[2*lam[[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]]]]]