##############################################################################
##
##  This is a collection of functions dealing with the ALGEBRA OF INVARIANTS
##  of a finite subgroup  G  of  GL(d,Z) acting in the canonical way on the 
##  GROUP ALGEBRA of the lattice  Z^d. The lattice elements are viewed as
##  rows and the matrices operate from the right.
##
##  The main functions are:
##
##  - Cl( G, n)  computes the class group of the algebra of invariants.
##               The input is  G  and  n = order of G . The output is a list
##               [q1,...,qr]  where the  qi  are prime powers. The class
##               group is isomorphic to the direct product of the groups
##               Z/qiZ.
##
##  - PicMat( G, n) does the same for the Picard group
##
##  - Reg( G, n) returns  `true`  if the invariant algebra is regular and
##               `false` otherwise. Input is the same as for  Cl( G, n).
##
##  The coefficient field  k  of the group algebra is assumed to have
##  characteristic not dividing the order of  G  (so the group algebra  kG
##  is semisimple). For simplicity, Cl( G, n) also assumes  k  to be large
##  enough so as to be a splitting field for the commutator factor group of
##  G/N  where  N  is the subgroup of  G  that is generated by the
##  reflections in  G  ( called  Refl(G)[1] below). No other information
##  about  k  is needed.
##
##############################################################################
##
##  Some functions dealing with the Hermite normal form of matrices
##  (contributed by Werner Nickel)
##
InfoHNF := Ignore;

IsDiagonalMat := function( M )
    
    local i, j;
    
    if Length(M) = 0 then return true; fi;
    if Length(M) <> Length(M[1]) then return false; fi;
    for i in [1..Length(M)] do
        for j in [i+1..Length(M[1])] do
            if M[i][j] <> 0 then return false; fi;
            if M[j][i] <> 0 then return false; fi;
        od;
    od;
    return true;
end;

RowEchelonForm := function( M )
    
    local a, i, j, k, m, n, r, Cleared;
    
    M := ShallowCopy( M );
    m := Length( M ); n := Length( M[1] );
    
    InfoHNF( "#I RowEchelonForm called with ", m, "x", n, "-matrix\n" );
    i := 1; j := 1;
    while i <= m and j <= n do
        k := i; while k <= m and M[k][j] = 0 do k := k+1; od;
        if k <= m then
            r := M[i]; M[i] := M[k]; M[k] := r;
            for k in [k+1..m] do
                a := AbsInt( M[k][j] );
                if a <> 0 and a < AbsInt( M[i][j] ) then
                    r := M[i]; M[i] := M[k]; M[k] := r;
                fi;
            od;
            if M[i][j] < 0 then M[i] := -1 * M[i]; fi;
            Cleared := true;
            for k in [i+1..m] do
                a := QuoInt(M[k][j],M[i][j]);
                if a <> 0 then  M[k] := M[k] - a * M[i]; fi;
                if M[k][j] <> 0 then Cleared := false; fi;
            od;
            if Cleared then i := i+1; j := j+1; fi;
        else 
            j := j+1;
        fi;
    od;
    return M{[1..i-1]};
end;

RowEchelonFormT := function( M, T )
    
    local a, i, j, k, m, n, r, Cleared;
    
    M := ShallowCopy( M );
    m := Length( M ); n := Length( M[1] );
    
    InfoHNF( "#I RowEchelonFormT called with ", m, "x", n, "-matrix\n" );
    i := 1; j := 1;
    while i <= m and j <= n do
        k := i; while k <= m and M[k][j] = 0 do k := k+1; od;
        if k <= m then
            r := M[i]; M[i] := M[k]; M[k] := r;
            r := T[i]; T[i] := T[k]; T[k] := r;
            for k in [k+1..m] do
                a := AbsInt( M[k][j] );
                if a <> 0 and a < AbsInt( M[i][j] ) then
                    r := M[i]; M[i] := M[k]; M[k] := r;
                    r := T[i]; T[i] := T[k]; T[k] := r;    
                fi;
            od;
            if M[i][j] < 0 then 
                M[i] := -1 * M[i]; 
                T[i] := -1 * T[i]; 
            fi;
            Cleared := true;
            for k in [i+1..m] do
                a := QuoInt(M[k][j],M[i][j]);
                if a <> 0 then 
                    M[k] := M[k] - a * M[i]; 
                    T[k] := T[k] - a * T[i];
                fi;
                if M[k][j] <> 0 then Cleared := false; fi;
            od;
            if Cleared then i := i+1; j := j+1; fi;
        else 
            j := j+1;
        fi;
    od;
    return M{[1..i-1]};
end;

DiagonalizedMat := function( M )
    
    M := RowEchelonForm(M);
    while not IsDiagonalMat(M) do
        M := TransposedMat(M);
        M := RowEchelonForm(M);
    od;
    
    return M;
end;

DiagonalizedMatT := function( M )
    local   T;
    
    T := IdentityMat( Length(M[1]) );
    M := RowEchelonForm(M);
    
    while not IsDiagonalMat(M) do
        M := TransposedMat(M);
        M := RowEchelonFormT(M,T);
        if not IsDiagonalMat(M) then
            M := TransposedMat(M);
            M := RowEchelonForm(M);
        fi;
    od;
    
    return [M, T];
end;

HermiteNormalForm := function( M )
    local   k,  h,  j;
    
    M := RowEchelonForm( M );
    
    for k in Reversed([1..Length(M)]) do
        h := 1; while M[k][h] = 0 do h := h+1; od;
        for j in [1..k-1] do
            M[j] := M[j] - Int( M[j][h]/M[k][h] ) * M[k];
        od;
    od;
    
    h := 1;
    for k in [1..Length(M)] do
        while M[k][h] = 0 do h := h+1; od;
        for j in [1..k-1] do
            M[j] := M[j] - Int( M[j][h]/M[k][h] ) * M[k];
            if M[j][h] < 0 then M[j] := M[j] + M[k]; fi;
        od;
    od;
    
    return M;
end;

SolveHomogenousEquations := function( M )
    local   Q,  solutions,  j,  v;
    
    Q := IdentityMat( Length(M[1]) );
    M := RowEchelonForm( M );
    
    while not IsDiagonalMat(M) do
        M := TransposedMat(M);
        M := RowEchelonFormT(M,Q);
        if not IsDiagonalMat(M) then
            M := TransposedMat(M);
            M := RowEchelonForm( M );
        fi;
    od;
    Q := TransposedMat( Q );

    #   Now we have  M * Q * X = b
    solutions := [];
    for j in [Length(M)+1..Length(Q)] do
        v := 0 * [1..Length(Q)]; v[j] := 1;
        Add( solutions, v );
    od;

    return HermiteNormalForm( List( solutions, w -> Q * w ) );
end;

SolveInhomogenousEquations := function( M, b )
    local   Q,  j,  lin,  v,  inhom;
    
    b := ShallowCopy(b);
    Q := IdentityMat( Length(M[1]) );
    M := RowEchelonFormT( M,b );
    
    while not IsDiagonalMat(M) do
        M := TransposedMat(M);
        M := RowEchelonFormT(M,Q);
        if not IsDiagonalMat(M) then
            M := TransposedMat(M);
            M := RowEchelonFormT(M,b);
        fi;
    od;
    Q := TransposedMat( Q );

    #   Now we have  M * Q * X = b
    for j in [Length(M)+1..Length(b)] do
        if b[j] <> 0 then
            return rec( lin := [],  inhom := [] );
        fi;
    od;

    lin := [];
    for j in [Length(M)+1..Length(Q)] do
        v := 0 * [1..Length(Q)]; v[j] := 1;
        Add( lin, v );
    od;

    v := 0 * [1..Length(Q)];
    for j in [1..Length(M)] do
        if b[j] mod M[j][j] <> 0 then
            return rec( lin := [],  inhom := [] );
        fi;
        v[j] := b[j] / M[j][j];
    od;

    return rec( lin := HermiteNormalForm(List( lin, w -> Q * w )),
                inhom := Q * v );
end;

ReduceVectorHNF := function( M, v )
    local   r,  h;
    h := 1;
    for r in M do
        while r[h] = 0 do h := h+1; od;
        v := v - Int( v[h] / r[h] ) * r;
        if v[h] < 0 then  v := v + r;  fi;
    od;
    return v;
end;

DecomposeVectorHNF := function( M, v )
    local   h,  w,  i,  r,  c;
    
    h := 1;
    w := [1..Length(M)] * 0;
    for i in [1..Length(M)] do
        r := M[i];
        while r[h] = 0 do h := h+1; od;
        c := v[h] / r[h];
        if not IsInt(c) then  Error(); fi;
        v := v - c * r;
        w[i] := c;
    od;
    return w;
end;
#############################################################################
##
#F  ZOneCohomologyMat( <G>, <|G|> )  . . . . . . . . . . . . . . . . H^1( G )
##  (contributed by Werner Nickel)
##
##  Compute the abelian invariants of H^1 for a subgroup G of GL( d, Z ) and
##  its action on the natural module. 
##
##  The second argument is the size of  G  because we have no good way
##  of computing the order of an integral matrix group.
##
ZOneCohomologyMat := function( G, n )
    local   d,  I,  S,  i,  j,  L;
    
    d := Length( G.generators[1] );
    I := IdentityMat( d );
    S := List( [1..d], i->[] );
    for i in [1..Length(G.generators)] do
        S{[1..d]}{[(i-1)*d+1..i*d]} := G.generators[i] - I;
    od;
    S := DiagonalizedMat( S );
    
    L := [];
    for i in [1..Length(S)] do
        d := GcdInt( S[i][i], n );
        if d <> 1 then Add( L, d ); fi;
    od;
    
    for i  in [1..Length(L)-1]  do
        for j  in [i+1..Length(L)]  do
            if L[j] mod L[i] <> 0  then
                d    := GcdInt( L[i], L[j] );
                L[j] := L[j] / d * L[i];
                L[i] := d;
            fi;
        od;
    od;
    return Filtered( L, x->x<>1 );
end;

##########################################################################
##
##  Inv( <G> ) . . . . . . . . . . . . . . . . . . . lattice invariants
##
##  Returns a Z-basis of the sublattice of  Z^d  consisting of the invariants
##  of a finite subgroup  G  of  GL(d,Z).
##
Inv:= function( G )

  local d, I, S, i;

  d := Length( G.generators[1] );
  I := IdentityMat( d );
  S := List( [1..d], i->[] );
  for i in [1..Length(G.generators)] do
      S{[1..d]}{[(i-1)*d+1..i*d]} := G.generators[i] - I;
  od;
  i:=[];
  if RankMat( S ) < d then
    i:= SolveHomogenousEquations( TransposedMat( S ) );
  fi;
  return i;
end;

##########################################################################
##
##  Refl( <G> ) . . . . . . . . . . . . . . . . . . . . .  reflections
##
##  Returns the normal subgroups of a finite subgroup  G  of  GL(d,Z)
##  that are generated by the reflections and the
##  diagonalizable (over the integers) reflections in  G, respectively.
##  The reflection subgroup is in Refl(G)[1], the diagonalizable reflection
##  subgroup in Refl(G)[2].
##
##  In order to distinguish between the two types of reflections, the
##  crucial observation is that, for  `d`  a diagonalizable reflection,
##  d-Id  is conjugate (over the integers) to  diag(-2,0,...,0) , and
##  for  `e` non-diagonalizable, e-Id  is conjugate to
##
##                 -2 1 0 ... 0
##                  0.........0
##                  ...........
##                  0.........0
##
##  In the first case, DiagonalizeMat(d-Id) returns  diag(2,0,...,0) and 
##  in the second case  DiagonalizeMat(e-Id) returns  diag(1,0,...,0).
##
Refl:= function( G )

  local refl, diagrefl, I, x, y;

  refl:=[]; diagrefl:=[];
  I:=IdentityMat( Length( G.generators[1] ) );
  for x in List(ConjugacyClasses(G), Representative) do
    if RankMat(x-I) <= 1 then Add(refl,x);fi;
  od;
  for x in refl do
    y:=x-I; DiagonalizeMat( y );
    if TraceMat( y ) = 2 then Add(diagrefl,x);fi;
  od;
  return [NormalClosure(G,Subgroup(G,refl)),
          NormalClosure(G,Subgroup(G,diagrefl)),refl,diagrefl];

end;

##############################################################################
##
##  Cl( <G>, <|G|> ) . . . . . . . . . . . . . . . . . class group
##
##  Returns the abelian invariants of the class group of the invariant algebra
##
Cl:= function( G, n )

  local  resgroup, classgroup, refl, U, V, mats, s;

  refl:=Refl(G);  s:=Size(refl[2]);
  classgroup:= AbelianInvariants(CommutatorFactorGroup(G/refl[1]));
  if s<n then
    if s=1 then Append( classgroup, ZOneCohomologyMat( G, n ) );
    else
      U:=Inv( refl[2] ); 
      if  U<>[] then
        V:=Inv( G );
        if V=[]  or  RankMat( U ) > RankMat( V )  then
          U:=RowSpace(Rationals, U);
          mats:=InducedActionSpaceMats( Basis(U), G.generators );
          resgroup:=Group( mats, IdentityMat( Length( mats[1] ) ) );
          Append( classgroup, ZOneCohomologyMat( resgroup, n ) );
        fi;
      fi;
    fi;
  fi;
  return classgroup;
end;

##############################################################################
##
##  Reg( <G>, <|G|> ) . . . . . . . . . . . . . . . . . regularity
##
##  Tests regularity of the invariant algebra
##
Reg:= function( G, n )

  local r;

  if n=Size(Refl(G)[1]) and Product(Cl( G, n ))=1 then r:=true;
  else r:=false;fi;
  return r;
end;

###############################################################################
##
##  PicMat( <G>, <|G|> ) . . . . . . . . . . . . . . . . . the Picard group
##
##  Returns the abelian invariants of the Picard group of the invariant algebra
##
PicMat:= function( G, n)
      local d, El, Ab, gen, D, g, j, C, L, c, pic, prod;

      d:=Length( G.generators[1] );
      El:=Filtered(GroupOps.Elements( G ) , x -> x <> G.identity);
      Ab:=AbelianGroup(AgWords, List([1..d], x->n) );
      gen:=Filtered(Ab.cgs, x -> Order(Ab,x)=n);
      D:=Ab;
      for  g in El  do
          C:=Subgroup( G, [g] );
          L:=[];
          for c in Inv(C) do
            prod:=gen[1]^c[1];
            for j in [2..d] do
              prod:=prod*gen[j]^c[j];
            od;
            Add(L, prod);
          od;
          D:=AgGroupOps.Intersection( D, Subgroup( Ab, L ) );
      od;
      L:=[];
      for c in Inv(G) do
        prod:=gen[1]^c[1];
        for j in [2..d] do
          prod:=prod*gen[j]^c[j];
        od;
        Add( L, prod );
      od;
      pic:=AgGroupOps.FactorGroup( D, Subgroup( Ab, L ) );
      return AgGroupOps.AbelianInvariants(pic);
end;