gap-libs 4r6p5-3 / usr / share / gap / lib / matint.gi

#############################################################################
##
#W  matint.gi                GAP library                        A. Storjohann
#W                                                              R. Wainwright
#W                                                                 F. Gähler
#W                                                                    D. Holt
#W                                                                  A. Hulpke
##
##
#Y  Copyright (C) 2003 The GAP Group
##
##  This file contains functions that compute Hermite and Smith normal forms
##  of integer matrices, with or without the HNF/SNF  expressed as the linear
##  combination of the input.
##

########################################################
##
##      auxiliary + main code for all in one function
##
##  MATINTsplit
##  MATINTrgcd
##  MATINTmgcdex
##  MATINTbezout
##  SNFofREF
##  NormalFormIntMat
##

########################################################
#
# MATINTsplit(<N>,<a>) - returns product of prime factors of N which are not factors of a.
#
BindGlobal("MATINTsplit",function(N,a)
local x,t;

  x:=a;
  t:=N;
  while x<>1 do
    x:=GcdInt(x,t);
    t:=QuoInt(t,x);
  od;
  return t;
end);

################################################
#
# MATINTrgcd(<N>,<a>) - Returns smallest nonnegative c such that gcd(N,a+c) = 1
#
BindGlobal("MATINTrgcd",function(N,a)
local k,r,d,i,c,g,q;

  if N=1 then return 0; fi;
  k := 1;
  r:=[(a-1) mod N];
  d:=[N];
  c:=0;
  while true do
    for i in [1..k] do r[i]:=(r[i]+1) mod d[i]; od;
    i:=PositionProperty(r,x->x<=0);
    if i=fail then
      g:=1;i:=0;
      while g=1 and i<k do
        i:=i+1;
        g:=GcdInt(r[i],d[i]);
      od;
      if g=1 then return c; fi;
      q:=MATINTsplit(QuoInt(d[i],g),g);
      if q>1 then
        k:=k+1;
        r[k]:=r[i] mod q;
        d[k]:=q;
      fi;
      r[i]:=0;
      d[i]:=g;
    fi;
    c:=c+1;
  od;

end);

#######################################################
#
#  MATINTmgcdex(<N>,<a>,<v>) - Returns c[1],c[2],...c[k] such that
#
#   gcd(N,a+c[1]*b[1]+...+c[n]*b[k]) = gcd(N,a,b[1],b[2],...,b[k])
#
BindGlobal("MATINTmgcdex", function(N,a,v)
local h,g,M,c,i,d,b,l;
  l:=Length(v);
  c:=[]; M:=[];
  h := N;
  for i in [1..l] do
    g := h;
    h:=GcdInt(g,v[i]);
    M[i]:=QuoInt(g,h);
  od;
  h:=GcdInt(a,h);
  g:=QuoInt(a,h);

  for i in [l,l-1..1] do
    b:=QuoInt(v[i],h);
    d:=MATINTsplit(M[i],b);
    if d=1 then
      c[i]:=0;
    else
      c[i]:=MATINTrgcd(d,g/b mod d);
      g:=g+c[i]*b;
    fi;
  od;

  return c;

end);




#####################################################
#
#  MATINTbezout(a,b,c,d) - returns row transform , P, to transform, A, to hnf :
#
#  PA=H;
#
#  [ s  t ] [ a  b ]   [ e  f ]
#  [      ] [      ] = [       ]
#  [ u  v ] [ c  d ]   [    g ]
#
BindGlobal("MATINTbezout", function(a,b,c,d)
local e,f,g,q;

  e := Gcdex(a,c);
  f := e.coeff1*b+e.coeff2*d;
  g := e.coeff3*b+e.coeff4*d;
  if g<0 then
    e.coeff3 := -e.coeff3;
    e.coeff4 := -e.coeff4;
    g := -g;
  fi;
  if g>0 then
    q := QuoInt(f-(f mod g),g);
    e.coeff1 := e.coeff1-q*e.coeff3;
    e.coeff2 := e.coeff2-q*e.coeff4;
  fi;
  return e;

end);

#####################################################
##
## SNFofREF - fast SNF of REF matrix
##
##
InstallGlobalFunction(SNFofREF , function(R,destroy)
local k,g,b,ii,m1,m2,t,tt,si,n,m,i,j,r,jj,piv,d,gt,tmp,A,T,TT,kk;

  Info(InfoMatInt,1,"SNFofREF - initializing work matrix");
  n := Length(R);
  m := Length(R[1]);

  piv := List(R,x->PositionProperty(x,y->y<>0));
  r := PositionProperty(piv,x->x=fail);
  if r=fail then
    r := Length(piv);
  else
    r := r-1;
    piv := piv{[1..r]};
  fi;
  Append(piv,Difference([1..m],piv));

  if destroy then
    T:=R;
    ##  Need to be careful: we are trying to permute the cols in place
    for i in [1..r] do
      T[i]{[1..m]} := T[i]{piv};
    od;
  else
    T := NullMat(n,m);
    for j in [1..m] do
      for i in [1..Minimum(r,j)] do
        T[i][j]:=R[i][piv[j]];
      od;
    od;
  fi;

  si := 1;
  A := [];
  d := 2;
  for k in [1..m] do
    Info(InfoMatInt,2,"SNFofREF - working on column ",k);
    if k<=r then
      d := d*AbsInt(T[k][k]);
      Apply(T[k],x->x mod (2*d));
    fi;

    t := Minimum(k,r);
    for i in [t-1,t-2..si] do
      t := MATINTmgcdex(A[i],T[i][k],[T[i+1][k]])[1];
      if t<>0 then
        AddRowVector(T[i],T[i+1],t);
        Apply(T[i],x->x mod A[i]);
      fi;
    od;

    for i in [si..Minimum(k-1,r)] do
      g := Gcdex(A[i],T[i][k]);
      T[i][k] := 0;
      if g.gcd<>A[i] then
        b := QuoInt(A[i],g.gcd);
        A[i] := g.gcd;
        for ii in [i+1..Minimum(k-1,r)] do
          AddRowVector(T[ii],T[i],-g.coeff2*QuoInt(T[ii][k],A[i]) mod A[ii]);
          T[ii][k] := b*T[ii][k];

          Apply(T[ii],x->x mod A[ii]);
        od;
        if k<=r then
          t := g.coeff2*QuoInt(T[k][k],g.gcd);
          AddRowVector(T[k],T[i],-t);
          T[k][k]:=b*T[k][k];
        fi;
        Apply(T[i],x->x mod A[i]);
        if A[i]=1 then si := i+1; fi;
      fi;
    od;

    if k<=r then
      A[k] := AbsInt(T[k][k]);
      Apply(T[k],x->x mod A[k]);
    fi;

  od;

  for i in [1..r] do T[i][i] := A[i]; od;

  return T;

end);

BindGlobal("BITLISTS_NFIM",
  [ [ false, false, false, false, false ], [ true, false, false, false, false ],
    [ false, true, false, false, false ], [ true, true, false, false, false ],
    [ false, false, true, false, false ], [ true, false, true, false, false ],
    [ false, true, true, false, false ], [ true, true, true, false, false ],
    [ false, false, false, true, false ], [ true, false, false, true, false ],
    [ false, true, false, true, false ], [ true, true, false, true, false ],
    [ false, false, true, true, false ], [ true, false, true, true, false ],
    [ false, true, true, true, false ], [ true, true, true, true, false ],
    [ false, false, false, false, true ], [ true, false, false, false, true ],
    [ false, true, false, false, true ], [ true, true, false, false, true ],
    [ false, false, true, false, true ], [ true, false, true, false, true ],
    [ false, true, true, false, true ], [ true, true, true, false, true ],
    [ false, false, false, true, true ], [ true, false, false, true, true ],
    [ false, true, false, true, true ], [ true, true, false, true, true ],
    [ false, false, true, true, true ], [ true, false, true, true, true ],
    [ false, true, true, true, true ], [ true, true, true, true, true ] ] );


###########################################################
#
# DoNFIM(<mat>,<options>)
#
# Options bit values:
#
# 1  - Triangular / Smith
# 2  - No / Yes  Reduce off diag entries
# 4  - No / Yes  Row Transforms
# 8  - No / Yes  Col Transforms
# 16 - change original matrix in place (The rows still change) -- save memory
#
# Compute a Triangular, Hermite or Smith form of the n x m
# integer input matrix A.  Optionally, compute n x n / m x m
# unimodular transforming matrices which satisfy Q C A = H
# or  Q C A B P = S.
#
# Triangular / Hermite :
#
# Let I be the min(r+1,n) x min(r+1,n) identity matrix with r = rank(A).
# Then Q and C can be written using a block decomposition as
#
#             [ Q1 |   ]  [ C1 | C2 ]
#             [----+---]  [----+----]  A  =  H.
#             [ Q2 | I ]  [    | I  ]
#
# Smith :
#
#  [ Q1 |   ]  [ C1 | C2 ]     [ B1 |   ]  [ P1 | P2 ]
#  [----+---]  [----+----]  A  [----+---]  [----+----] = S.
#  [ Q2 | I ]  [    | I  ]     [ B2 | I ]  [ *  | I  ]
#
# * - possible non-zero entry in upper right corner...
#
#
BindGlobal("DoNFIM", function(arg)
local opt, sig, n, m, A, C, Q, B, P, r, c2, rp, c1, j, k, N, L, b, a, g, c,
      t, tmp, i, q, R, rank, signdet;

  if not Length(arg)=2
     or not (IsMatrix(arg[1])
         or (IsList(arg[1]) and Length(arg[1])=1
             and IsList(arg[1][1]) and Length(arg[1][1])=0))
     or not IsInt(arg[2]) then
    Error("syntax is DoNFIM(<matrix>,<options>)");
  fi;

  #Parse options
  opt := BITLISTS_NFIM[arg[2]+1];
  #List(CoefficientsQadic(arg[2],2),x->x=1);
  #if Length(opt)<4 then
  #  opt{[Length(opt)+1..4]} := List([Length(opt)+1..4],x->false);
  #fi;

  sig:=1;

  #Embed arg[1] in 2 larger "id" matrix
  n := Length(arg[1])+2;
  m := Length(arg[1][1])+2;
  k:=ListWithIdenticalEntries(m,0);
  if opt[5] then
    # change the matrix
    A:=arg[1];
    for i in [n-1,n-2..2] do
      A[i]:=ShallowCopy(k);
      A[i]{[2..m-1]}:=A[i-1];
    od;
  else
    A := [];
    for i in [2..n-1] do
      #A[i] := [0];
      #Append(A[i],arg[1][i-1]);
      #A[i][m] := 0;
      A[i]:=ShallowCopy(k);
      A[i]{[2..m-1]}:=arg[1][i-1];
    od;
  fi;
  A[1]:=ShallowCopy(k);
  A[n]:=k;
  A[1][1] := 1;
  A[n][m] := 1;

  if opt[3] then
    C := IdentityMat(n);
    Q := NullMat(n,n);
    Q[1][1] := 1;
  fi;

  if opt[1] and opt[4] then
    B := IdentityMat(m);
    P := IdentityMat(m);
  fi;

  r := 0;
  c2 := 1;
  rp := [];
  while m>c2 do
    Info(InfoMatInt,2,"DoNFIM - reached column ",c2," of ",m);
    r := r+1;
    c1 := c2;
    rp[r] := c1;
    if opt[3] then Q[r+1][r+1] := 1; fi;

    j := c1+1;
    while j<=m do
      k := r+1;
      while k<=n and A[r][c1]*A[k][j]=A[k][c1]*A[r][j] do k := k+1; od;
      if k<=n then c2 := j; j := m; fi;
      j := j+1;
    od;
    #Smith with some transforms..
    if opt[1] and (opt[4] or opt[3]) and c2<m then
      N := Gcd(Flat(A{[r..n]}[c2]));
      L := [c1+1..c2-1];
      Append(L,[c2+1..m-1]);
      Add(L,c2);
      for j in L do
        if j=c2 then
          b:=A[r][c2];a:=A[r][c1];
          for i in [r+1..n] do
            if b<>1 then
              g:=Gcdex(b,A[i][c2]);
              b:=g.gcd;
              a:=g.coeff1*a+g.coeff2*A[i][c1];
            fi;
          od;
          N:=0;
          for i in [r..n] do
            if N<>1 then N:=GcdInt(N,A[i][c1]-QuoInt(A[i][c2],b)*a);fi;
          od;
        else
          c := MATINTmgcdex(N,A[r][j],A{[r+1..n]}[j]);
          b := A[r][j]+c*A{[r+1..n]}[j];
          a := A[r][c1]+c*A{[r+1..n]}[c1];
        fi;
        t := MATINTmgcdex(N,a,[b])[1];
        tmp := A[r][c1]+t*A[r][j];
        while tmp=0 or tmp*A[k][c2]=(A[k][c1]+t*A[k][j])*A[r][c2] do
          t := t+1+MATINTmgcdex(N,a+t*b+b,[b])[1];
          tmp := A[r][c1]+t*A[r][j];
        od;
        if t>0 then
          for i in [1..n] do A[i][c1] := A[i][c1]+t*A[i][j]; od;
          if opt[4] then B[j][c1] := B[j][c1]+t; fi;
        fi;
      od;
      if A[r][c1]*A[k][c1+1]=A[k][c1]*A[r][c1+1] then
        for i in [1..n] do A[i][c1+1] := A[i][c1+1]+A[i][c2]; od;
        if opt[4] then B[c2][c1+1] := 1; fi;
      fi;
      c2 := c1+1;
    fi;

    c := MATINTmgcdex(AbsInt(A[r][c1]),A[r+1][c1],A{[r+2..n]}[c1]);
    for i in [r+2..n] do
      if c[i-r-1]<>0 then
        AddRowVector(A[r+1],A[i],c[i-r-1]);
        if opt[3] then
          C[r+1][i] := c[i-r-1];
          AddRowVector(Q[r+1],Q[i],c[i-r-1]);
        fi;
      fi;
    od;

    i := r+1;
    while A[r][c1]*A[i][c2]=A[i][c1]*A[r][c2] do i := i+1; od;
    if i>r+1 then
      c := MATINTmgcdex(AbsInt(A[r][c1]),A[r+1][c1]+A[i][c1],[A[i][c1]])[1]+1;;
      AddRowVector(A[r+1],A[i],c);
      if opt[3] then
        C[r+1][i] := C[r+1][i]+c;
        AddRowVector(Q[r+1],Q[i],c);
      fi;
    fi;

    g := MATINTbezout(A[r][c1],A[r][c2],A[r+1][c1],A[r+1][c2]);
    sig:=sig*SignInt(A[r][c1]*A[r+1][c2]-A[r][c2]*A[r+1][c1]);
    A{[r,r+1]} := [[g.coeff1,g.coeff2],[g.coeff3,g.coeff4]]*A{[r,r+1]};
    if opt[3] then
      Q{[r,r+1]} := [[g.coeff1,g.coeff2],[g.coeff3,g.coeff4]]*Q{[r,r+1]};
    fi;

    for i in [r+2..n] do
      q := QuoInt(A[i][c1],A[r][c1]);
      AddRowVector(A[i],A[r],-q);
      if opt[3] then AddRowVector(Q[i],Q[r],-q); fi;
      q := QuoInt(A[i][c2],A[r+1][c2]);
      AddRowVector(A[i],A[r+1],-q);
      if opt[3] then AddRowVector(Q[i],Q[r+1],-q); fi;
    od;

  od;
  rp[r+1] := m;
  Info(InfoMatInt,2,"DoNFIM - r,m,n=",r,m,n);
  if n=m and r+1<n then sig:=0;fi;

  #smith w/ NO transforms - farm the work out...
  if opt[1] and not (opt[3] or opt[4]) then
    #R:=rec(normal:=SNFofREF(A{[2..n-1]}{[2..m-1]}),rank:=r-1);
    for i in [2..n-1] do
      A[i-1]:=A[i]{[2..m-1]};
    od;
    Unbind(A[n-1]);
    Unbind(A[n]);
    R:=rec(normal:=SNFofREF(A,opt[5]),rank:=r-1);
    if n=m then R. signdet:=sig;fi;
    return R;
  fi;

  # hermite or (smith w/ column transforms)
  if (not opt[1] and opt[2]) or (opt[1] and opt[4]) then
    for i in [r, r-1 .. 1] do
      Info(InfoMatInt,2,"DoNFIM - reducing row ",i);
      for j in [i+1 .. r+1] do
        q := QuoInt(A[i][rp[j]]-(A[i][rp[j]] mod A[j][rp[j]]),A[j][rp[j]]);
        AddRowVector(A[i],A[j],-q);
        if opt[3] then AddRowVector(Q[i],Q[j],-q); fi;
      od;
      if opt[1] and i<r then
        for j in [i+1..m] do
          q := QuoInt(A[i][j],A[i][i]);
          for k in [1..i] do A[k][j] := A[k][j]-q*A[k][i]; od;
          if opt[4] then P[i][j] := -q; fi;
        od;
      fi;
    od;
  fi;

  #Smith w/ row but not col transforms
  if opt[1] and opt[3] and not opt[4] then
    for i in [1..r-1] do
      t := A[i][i];
      A[i] := List([1..m],x->0);
      A[i][i] := t;
    od;
    for j in [r+1..m-1] do
      A[r][r] := GcdInt(A[r][r],A[r][j]);
      A[r][j] := 0;
    od;
  fi;

  #smith w/ col transforms
  if opt[1] and opt[4] and r<m-1 then
    c := MATINTmgcdex(A[r][r],A[r][r+1],A[r]{[r+2..m-1]});
    for j in [r+2..m-1] do
      A[r][r+1] := A[r][r+1]+c[j-r-1]*A[r][j];
      B[j][r+1] := c[j-r-1];
      for i in [1..r] do P[i][r+1] := P[i][r+1]+c[j-r-1]*P[i][j]; od;
    od;
    P[r+1] := List([1..m],x->0);
    P[r+1][r+1] := 1;
    g := Gcdex(A[r][r],A[r][r+1]);
    A[r][r] := g.gcd;
    A[r][r+1] := 0;
    for i in [1..r+1] do
      t := P[i][r];
      P[i][r] := P[i][r]*g.coeff1+P[i][r+1]*g.coeff2;
      P[i][r+1] := t*g.coeff3+P[i][r+1]*g.coeff4;
    od;
    for j in [r+2..m-1] do
      q := QuoInt(A[r][j],A[r][r]);
      for i in [1..r+1] do P[i][j] := P[i][j]-q*P[i][r]; od;
      A[r][j] := 0;
    od;
    for i in [r+2..m-1] do
      P[i] := List([1..m],x->0);
      P[i][i] := 1;
    od;
  fi;

  #row transforms finisher
  if opt[3] then for i in [r+2..n] do Q[i][i]:= 1; od; fi;

  for i in [2..n-1] do
    A[i-1]:=A[i]{[2..m-1]};
  od;
  Unbind(A[n-1]);
  Unbind(A[n]);
  R:=rec(normal:=A);

  if opt[3] then
    R.rowC:=C{[2..n-1]}{[2..n-1]};
    R.rowQ:=Q{[2..n-1]}{[2..n-1]};
  fi;

  if opt[1] and opt[4] then
    R.colC:=B{[2..m-1]}{[2..m-1]};
    R.colQ:=P{[2..m-1]}{[2..m-1]};
  fi;

  R.rank:=r-1;
  if n=m then R.signdet:=sig;fi;
  return R;

end);

#############################################################################
##
#F  NormalFormIntMat(<mat>,<options>)
##
InstallGlobalFunction(NormalFormIntMat,
function(mat,options)
  local r,opt;
  r:=DoNFIM(mat,options);
  opt := BITLISTS_NFIM[options+1];
  #opt := List(CoefficientsQadic(options,2),x->x=1);
  #if Length(opt)<4 then
  #  opt{[Length(opt)+1..4]} := List([Length(opt)+1..4],x->false);
  #fi;

  if opt[3] then
    r.rowtrans:=r.rowQ*r.rowC;
    #Unbind(r.rowQ);
    #Unbind(r.rowC);
  fi;

  if opt[1] and opt[4] then
    r.coltrans:=r.colC*r.colQ;
    #Unbind(r.colQ);
    #Unbind(r.colC);
  fi;
  return r;
end);

#############################################################################
##
#O  TriangulizedIntegerMat(<mat>);
##
InstallMethod(TriangulizedIntegerMat,"dispatch",true,[IsMatrix],0,
function(mat)
  return DoNFIM(mat,0).normal;
end);

InstallOtherMethod(TriangulizedIntegerMat,"empty matrix",true,[IsList],0,
function(mat)
  if Length(mat)<>1 or (not IsList(mat[1])) or Length(mat[1])<>0 then
    TryNextMethod();
  fi;
  return DoNFIM(mat,0).normal;
end);
InstallOtherMethod(TriangulizedIntegerMat,"empty",true,[IsEmpty],0,Immutable);

#############################################################################
##
#O  TriangulizedIntegerMatTransform(<mat>);
##
InstallMethod(TriangulizedIntegerMatTransform,"dispatch",true,[IsMatrix],0,
function(mat)
  return NormalFormIntMat(mat,4);
end);

InstallOtherMethod(TriangulizedIntegerMatTransform,"empty matrix",true,[IsList],0,
function(mat)
  if Length(mat)<>1 or (not IsList(mat[1])) or Length(mat[1])<>0 then
    TryNextMethod();
  fi;
  return NormalFormIntMat(mat,4);
end);
InstallOtherMethod(TriangulizedIntegerMatTransform,"empty",true,[IsEmpty],0,
function(mat)
  return rec(normal:=Immutable(mat),rowtrans:=Immutable([[1]]));
end);

#############################################################################
##
#O  TriangulizeIntegerMat(<mat>);
##
InstallMethod(TriangulizeIntegerMat,"dispatch",true,[IsMatrix and IsMutable],0,
function(mat)
  DoNFIM(mat,16);
end);

InstallOtherMethod(TriangulizeIntegerMat,"empty",true,[IsEmpty],0,Immutable);

#############################################################################
##
#O  HermiteNormalFormIntegerMat(<mat>);
##
InstallMethod(HermiteNormalFormIntegerMat,"dispatch",true,[IsMatrix],0,
function(mat)
  return DoNFIM(mat,2).normal;
end);

InstallOtherMethod(HermiteNormalFormIntegerMat,"empty matrix",true,[IsList],0,
function(mat)
  if Length(mat)<>1 or (not IsList(mat[1])) or Length(mat[1])<>0 then
    TryNextMethod();
  fi;
  return DoNFIM(mat,2).normal;
end);
InstallOtherMethod(HermiteNormalFormIntegerMat,"empty",true,[IsEmpty],0,
  Immutable);

#############################################################################
##
#O  HermiteNormalFormIntegerMatTransform(<mat>);
##
InstallMethod(HermiteNormalFormIntegerMatTransform,"dispatch",true,[IsMatrix],0,
function(mat)
  return NormalFormIntMat(mat,6);
end);

InstallOtherMethod(HermiteNormalFormIntegerMatTransform,"empty matrix",
  true,[IsList],0,
function(mat)
  if Length(mat)<>1 or (not IsList(mat[1])) or Length(mat[1])<>0 then
    TryNextMethod();
  fi;
  return NormalFormIntMat(mat,6);
end);
InstallOtherMethod(HermiteNormalFormIntegerMatTransform,"empty",true,
  [IsEmpty],0,
function(mat)
  return rec(normal:=Immutable(mat),rowtrans:=Immutable([[1]]));
end);

#############################################################################
##
#O  SmithNormalFormIntegerMat(<mat>);
##
InstallMethod(SmithNormalFormIntegerMat,"dispatch",true,[IsMatrix],0,
function(mat)
  return DoNFIM(mat,1).normal;
end);

InstallOtherMethod(SmithNormalFormIntegerMat,"empty matrix",true,[IsList],0,
function(mat)
  if Length(mat)<>1 or (not IsList(mat[1])) or Length(mat[1])<>0 then
    TryNextMethod();
  fi;
  return DoNFIM(mat,1).normal;
end);
InstallOtherMethod(SmithNormalFormIntegerMat,"empty",true,[IsEmpty],0,
  Immutable);

#############################################################################
##
#O  SmithNormalFormIntegerMatTransforms(<mat>);
##
InstallMethod(SmithNormalFormIntegerMatTransforms,"dispatch",true,[IsMatrix],0,
function(mat)
  return NormalFormIntMat(mat,13);
end);

InstallOtherMethod(SmithNormalFormIntegerMatTransforms,"empty matrix",
  true,[IsList],0,
function(mat)
  if Length(mat)<>1 or (not IsList(mat[1])) or Length(mat[1])<>0 then
    TryNextMethod();
  fi;
  return NormalFormIntMat(mat,13);
end);
InstallOtherMethod(SmithNormalFormIntegerMatTransforms,"empty",true,
  [IsEmpty],0,
function(mat)
  return
  rec(normal:=Immutable(mat),rowtrans:=Immutable([[1]]),
      coltrans:=Immutable([[1]]));
end);

InstallGlobalFunction( DiagonalizeIntMat, function ( mat )
  DoNFIM(mat,17);
end);

#############################################################################
##
#M  DiagonalizeMat(<integers>,<mat>)
##
InstallMethod( DiagonalizeMat, "over the integers",
  [ IsIntegers, IsMatrix and IsMutable ],
function(I,mat)
  DiagonalizeIntMat(mat);
end );

#############################################################################
##
#M  ElementaryDivisorsTransformationsMat(<integers>,<mat>)
##
InstallMethod( ElementaryDivisorsTransformationsMat, "over the integers",
  [ IsIntegers, IsMatrix and IsMutable ],
function(I,mat)
  return SmithNormalFormIntegerMatTransforms(mat);
end );

#############################################################################
##
#M  BaseIntMat(<mat>)
##
InstallMethod(BaseIntMat,"use HNF",true,
  [IsMatrix and IsCyclotomicCollColl],0,
function( mat )
local norm;
  norm := NormalFormIntMat( mat, 2 );
  return norm.normal{[1..norm.rank]};
end);

InstallOtherMethod(BaseIntMat,"empty",true,
  [IsEmpty],0,Immutable);

#############################################################################
##
#M  BaseIntersectionIntMats(<m1>,<m2>)
##
InstallMethod(BaseIntersectionIntMats,"use HNF",true,
  [IsMatrix and IsCyclotomicCollColl,IsMatrix and IsCyclotomicCollColl],0,
function( M1, M2 )
local M, Q, r, T;
  M := Concatenation( M1, M2 );
  r := NormalFormIntMat( M, 4 );
  T := r.rowtrans{[r.rank+1..Length(M)]}{[1..Length(M1)]};
  if not IsEmpty( T ) then T := T * M1; fi;
  return BaseIntMat( T );
end);

InstallOtherMethod(BaseIntersectionIntMats,"emptyL",true,
  [IsEmpty,IsObject],0,
function(L,R)
  return Immutable(L);
end);

InstallOtherMethod(BaseIntersectionIntMats,"emptyR",true,
  [IsObject,IsEmpty],0,
function(L,R)
  return Immutable(R);
end);

#############################################################################
##
#M  ComplementIntMat(<m1>,<m2>)
##
InstallMethod(ComplementIntMat,"use HNF and SNF",true,
  [IsMatrix and IsCyclotomicCollColl,IsMatrix and IsCyclotomicCollColl],0,
function( full,sub )
local F, S, M, r, T, R;
  F := BaseIntMat( full );
  if IsEmpty( sub ) or IsZero( sub ) then
    return rec( complement := F, sub := [], moduli := [] );
  fi;
  S := BaseIntersectionIntMats( F, sub );
  if S <> BaseIntMat( sub ) then
    Error( "sub must be submodule of full" );
  fi;
  # find T such that BaseIntMat(T*F) = S
  M := Concatenation( F, S );
  T := NormalFormIntMat( M, 4 ).rowtrans^-1;
  T := T{[Length(F)+1..Length(T)]}{[1..Length(F)]};

  # r.rowtrans * T * F = r.normal * r.coltrans^-1 * F
  r := NormalFormIntMat( T, 13 );
  M := r.coltrans^-1 * F;
  R := rec( complement := BaseIntMat( M{[1+r.rank..Length(M)]} ),
            sub := r.rowtrans * T * F,
            moduli := List( [1..r.rank], i -> r.normal[i][i] ) );
  return R;
end);

InstallOtherMethod(ComplementIntMat,"empty submodule",true,
  [IsMatrix and IsCyclotomicCollColl,IsList and IsEmpty],0,
function( full, sub )
  return rec( complement := BaseIntMat( full ), sub := [], moduli := [] );
end );

#############################################################################
##
#M  NullspaceIntMat(<mat>)
##
InstallMethod(NullspaceIntMat,"use HNF",true,
  [IsMatrix and IsCyclotomicCollColl],0,
function( mat )
local norm, kern;
  norm := NormalFormIntMat( mat, 4 );
  kern := norm.rowtrans{[norm.rank+1..Length(mat)]};
  return BaseIntMat( kern );
end);

#############################################################################
##
#M  SolutionIntMat(<mat>,<vec>)
##
InstallMethod(SolutionIntMat,"use HNF",true,
  [IsMatrix and IsCyclotomicCollColl,
   IsList and IsCyclotomicCollection],0,
function( mat,v )
local norm, rs, t, M, r;
  if IsZero(mat) then
    if IsZero(v) then
      return ListWithIdenticalEntries( Length(mat), 0 );
    else
      return fail;
    fi;
  fi;
  norm := NormalFormIntMat( mat, 4 );
  t := norm.rowtrans;
  rs :=  norm.normal{[1..norm.rank]};
  M := Concatenation( rs, [v] );
  r := NormalFormIntMat( M, 4 );
  if r.rank = Length(r.normal) or
     r.rowtrans[Length(M)][Length(M)] <> 1 then
    return fail;
  fi;
  return -r.rowtrans[Length(M)]{[1..r.rank]} * t{[1..r.rank]};
end);

InstallOtherMethod(SolutionIntMat,"empty",true,[IsEmpty,IsObject],0,
function(mat,v)
  return fail;
end);

#############################################################################
##
#M  SolutionNullspaceIntMat(<mat>,<vec>)
##
InstallMethod(SolutionNullspaceIntMat,"use HNF",true,
  [IsMatrix and IsCyclotomicCollColl,
   IsList and IsCyclotomicCollection],0,
function( mat,v )
local norm, rs, t, M, r, kern, len;
  if IsZero(mat) then
    len := Length(mat);
    if IsZero(v) then
      return [ListWithIdenticalEntries(len,0), IdentityMat(len)];
    else
      return [fail, IdentityMat(len)];
    fi;
  fi;
  norm := NormalFormIntMat( mat, 4 );
  kern := norm.rowtrans{[norm.rank+1..Length(mat)]};
  kern := BaseIntMat( kern );
  t := norm.rowtrans;
  rs :=  norm.normal{[1..norm.rank]};
  M := Concatenation( rs, [v] );
  r := NormalFormIntMat( M, 4 );
  if r.rank = Length(r.normal) or
     r.rowtrans[Length(M)][Length(M)] <> 1 then
    return [fail,kern];
  fi;
  return [-r.rowtrans[Length(M)]{[1..r.rank]} * t{[1..r.rank]}, kern];
end);


#############################################################################
##
#F  DeterminantIntMat(<mat>)
##
InstallGlobalFunction(DeterminantIntMat,function(mat)
local sig, n, m, A, r, c2, c1, j, k, c, i, g, q;

  sig:=1;

  #Embed mat in 2 larger "id" matrix
  n := Length(mat)+2;
  # Crossover point roughly 20x20 matrices, so farm the work if smaller..
  if n<22 then return DeterminantMat(mat);fi;
  m := Length(mat[1])+2;

  if not n=m then
    Error( "DeterminantIntMat: <mat> must be a square matrix" );
  fi;

  A := [List([1..m],x->0)];
  for i in [2..n-1] do
    A[i] := [0];
    Append(A[i],mat[i-1]);
    A[i][m] := 0;
  od;
  A[n] := List([1..m],x->0);
  A[1][1] := 1;      A[n][m] := 1;

  r := 0;    c2 := 1;
  while m>c2 do
    Info(InfoMatInt,2,"DeterminantIntMat - reached column ",c2," of ",m);
    r := r+1;
    c1 := c2;

    j := c1+1;
    while j<=m do
      k := r+1;
      while k<=n and A[r][c1]*A[k][j]=A[k][c1]*A[r][j] do k := k+1; od;
      if k<=n then c2 := j; j := m; fi;
      j := j+1;
    od;

    c := MATINTmgcdex(AbsInt(A[r][c1]),A[r+1][c1],A{[r+2..n]}[c1]);
    for i in [r+2..n] do
      if c[i-r-1]<>0 then
        AddRowVector(A[r+1],A[i],c[i-r-1]);
      fi;
    od;

    i := r+1;
    while A[r][c1]*A[i][c2]=A[i][c1]*A[r][c2] do
      i := i+1;
    od;

    if i>r+1 then
      c := MATINTmgcdex(AbsInt(A[r][c1]),A[r+1][c1]+A[i][c1],[A[i][c1]])[1]+1;;
      AddRowVector(A[r+1],A[i],c);
    fi;

    g := MATINTbezout(A[r][c1],A[r][c2],A[r+1][c1],A[r+1][c2]);
    sig:=sig*SignInt(A[r][c1]*A[r+1][c2]-A[r][c2]*A[r+1][c1]);
    if sig=0 then return 0;fi;
    A{[r,r+1]} := [[g.coeff1,g.coeff2],[g.coeff3,g.coeff4]]*A{[r,r+1]};

    for i in [r+2..n] do
      q := QuoInt(A[i][c1],A[r][c1]);
      AddRowVector(A[i],A[r],-q);
      q := QuoInt(A[i][c2],A[r+1][c2]);
      AddRowVector(A[i],A[r+1],-q);
    od;
  od;

  for i in [2..r+1] do
    sig:=sig*A[i][i];
  od;

  return sig;

end);

#############################################################################
##
#M  AbelianInvariantsOfList( <list> ) . . . . .  abelian invariants of a list
##
InstallMethod( AbelianInvariantsOfList,
    [ IsCyclotomicCollection ],
function ( list )
    local   invs, elm;

    invs := [];
    for elm  in list  do
        if elm = 0  then
            Add( invs, 0 );
        elif 1 < elm  then
            Append( invs, List( Collected(Factors(elm)), x->x[1]^x[2] ) );
        elif elm < -1 then
            Append( invs, List( Collected(Factors(-elm)), x->x[1]^x[2] ) );
        fi;
    od;
    Sort(invs);
    return invs;
end );

InstallOtherMethod( AbelianInvariantsOfList,
    [ IsList and IsEmpty ],
    list -> [] );