// Functions used in paper 8:
// _Support varieties for modules_ by Jon F. Carlson

// Section 2: Notes on projectivity

IsProjective := function(M, B, n)
  V := VectorSpace(BaseRing(B[1]), Dimension(M));
  S := &+[ sub< V | RowSpace(x) > : x in B ];
  if (Dimension(M) - Dimension(S))*n eq Dimension(M) then
    return true;
  end if;
  return false;
end function;

LeftProjectiveStrip := function(M, G)
  MM := Restriction(M, G);
  repeat
    flag := true;
    for i := 1 to 5 do
      a := Random(MM);
      SM := sub< MM | a >;
      if Dimension(SM) eq #G then
        MM := quo< MM | SM >;
        flag := false;
        break;
      end if;
    end for;
  until flag;
  return MM;
end function;


// Section 4: Finding points on the variety

DimensionOfVariety := function(M, G, num)
  phi := Representation(M);
  n := Ngens(G);
  dim := Dimension(M);
  ff := ext< CoefficientRing(M) | 6 >;
  NM := ExtendField(M, ff);
  V := VectorSpace(ff, n);
  MA := MatrixAlgebra(ff, dim);
  newmatlist := [ MA!Representation(NM)(G.i) - Id(MA) : i in [1..n] ];
  for r := 1 to n do
    g := AbelianGroup(GrpPerm, [ 2 : i in [1..r] ]);
    RRL := [ ];
    flag := true;
    for j := 1 to num do
      a := [ Random(V) : i in [1..r] ];
      mmmlist := [ Id(MA) + &+[ a[i][j]*newmatlist[j] :
                   j in [1..n] ] : i in [1..r] ];
      mmm := GModule(g, mmmlist);
      str := LeftProjectiveStrip(mmm, g);
      if Dimension(str) eq 0 and r lt n then
        flag := false;
        break j;
      end if;
      if Dimension(str) eq 0 and r eq n then
        return g, [ <a, str> ], newmatlist;
      end if;
      RRL := Append(RRL, <a, str>);
      if #RRL eq 8 and r eq 1 then
        return g, RRL, newmatlist;
      end if;
    end for;
    if flag then
      return g, RRL, newmatlist;
    end if;
  end for;
end function;

StandardForm := function(M)
  r := Ngens(Group(M));
  id := MatrixAlgebra(BaseRing(M), Dimension(M))!1;
  LL := [ ActionGenerator(M, i) - id : i in [1..r] ];
  W := sub< M | &cat[ &cat[ Basis(RowSpace(LL[i]*LL[j])) :
               j in [1..i] ]: i in [1..r-1] ] >;
  MM := quo< M | W >;
  idd := MatrixAlgebra(BaseRing(MM), Dimension(MM)) ! 1;
  Rad := sub< MM | &cat[ Basis(RowSpace(ActionGenerator(MM, i) - idd)) :
                         i in [1..r-1] ] >;
  Gen := [ MM ! ExtendBasis(Rad, MM)[i] :
                i in [ Dimension(Rad)+1 .. Dimension(MM) ] ];
  Bas := Gen cat &cat[ [ x*ActionGenerator(MM, i) - x : x in Gen ] :
                       i in [1..r-1] ];
  b := Matrix(Dimension(MM), Dimension(MM), [ Vector(x) : x in Bas ]);
  NL := [ b*ActionGenerator(MM, i)*b^-1 : i in [1..r] ];
  return NL;
end function;

PointsOnVariety := function(M)
  G := Group(M);
  PL := [ ];
  BL := [ ];
  Npoints := 2^Ngens(G)*(Dimension(M) div #G + 1) + 20;
  GG, vl, ACT := DimensionOfVariety(M, G, Npoints);
  r := Ngens(GG);
  if r eq 1 then
   return [ ], 0, 0;
  end if;
  ff := CoefficientRing(vl[1][2]);
  if Dimension(vl[1][2]) eq 0 then
   return [0], Ngens(G), 0;
  end if;
  mini := Minimum([ Dimension(x[2]) : x in vl ]);
  w := mini div (2^(r-1));
  for x in vl do
    mm := x[2];
    if Dimension(mm) eq mini then
      bool := IsProjective(mm, [ ActionGenerator(mm, i)
              - MatrixAlgebra(ff, mini) ! 1 : i in [1..r-1] ], 2^(r-1));
      if bool then
        SF := StandardForm(mm);
        blocks := [ Submatrix(SF[r], 1, w*i + 1, w, w) : i in [1..r-1] ];
        minpolys := [ MinimalPolynomial(y) : y in blocks ];
        Append(~PL, <x[1], minpolys>);
        Append(~BL, blocks);
      end if;
    end if;
  end for;
  degrees := [ Degree(Factorization(PL[i][2][j])[k][1]) : k in
       [1..#Factorization(PL[i][2][j])], j in [1..#PL[i][2]], i in [1..#PL] ];
  xd := LCM(degrees);
  if xd eq 1 then
    RL := [ <PL[i][1], [ Roots(PL[i][2][j]) :
                         j in [1..#PL[i][2]] ]>: i in [1..#PL] ];
    NPL := PL;
    fff := ff;
  else
    NPL:= [ ];
    RL := [ ];
    fff, gamma := ext< ff | xd >;
    VV := VectorSpace(fff, Ngens(G));
    Q := PolynomialRing(fff);
    for i := 1 to #PL do
      npl := <[ VV ! a : a in PL[i][1] ], [ Q!f : f in PL[i][2] ]>;
      Append(~NPL, npl);
      rl := <npl[1], [ Roots(y): y in npl[2] ]>;
      Append(~RL, rl);
    end for;
  end if;
  VL1 := [ ];
  for w in RL do
    rt := w[2];
    rl := [ [ rt[j][k][1] : k in [1..#rt[j]] ] : j in [1..#rt] ];
    if #rl eq 1 then
      RL1 := [ [x,1] : x in rl[1] ];
    else
      RL1 := [ [x] : x in rl[1] ];
      for i := 2 to #rl do
        RL1 := [ Append(x, y) : x in RL1, y in rl[i] ];
      end for;
      RL1 := [ Append(x, 1) : x in RL1 ];
    end if;
    RL2 := [ &+[ w[1][j]*RL1[k][j] : j in [1..r] ] : k in [1..#RL1] ];
    VL1 := VL1 cat RL2;
  end for;
  VL := [ ]; m := Nrows(ACT[1]);
  if xd gt 1 then
    NACT := [ MatrixAlgebra(fff, m) ! x : x in ACT ];
  else
    NACT := ACT;
  end if;
  for w in VL1 do
    MAT := &+[ w[i]*NACT[i] : i in [1..Ngens(G)] ];
    if 2*Rank(MAT) lt m then
      Append(~VL, w);
    end if;
   end for;
   return VL, r-1, w;
end function;


// Section 5: Computing the variety from a set of points

VarietyPolynomials := function(P, pts, deg)
  V := Parent(pts[1]);
  F := CoefficientField(V);
  Mons := MonomialsOfDegree(P, deg);
  M := RMatrixSpace(F, #Mons, #pts);
  phi := [ hom< P -> F | ElementToSequence(pts[i]) > : i in [1..#pts] ];
  A := M ! &cat[ [ phi[i](Mons[j]) : i in [1..#pts], j in [1..#Mons] ] ];
  _, Null := Solution(A, Zero(VectorSpace(F, #pts)));
  return [ &+[ Basis(Null)[j][i] * Mons[i] :
       i in [1.. #Mons] ] : j in [1..Dimension(Null)] ];
end function;

SupportVariety := function(M)
  G := Group(M);
  MM := LeftProjectiveStrip(M, G);
  if Dimension(MM) eq 0 then
    P := PolynomialRing(CoefficientRing(M), Ngens(G));
    return [ P.i: i in [1..Rank(P)] ], P;
  end if;
  pts, codim, deg := PointsOnVariety(MM);
  if #pts eq 0 then
    P := PolynomialRing(CoefficientRing(M), Ngens(G));
    return [ ], P;
  end if;
  s := Ngens(G);
  ff := BaseRing(pts[1]);
  P := PolynomialRing(ff, s);
  POLS := &cat[ VarietyPolynomials(P, pts, t) : t in [1..deg] ];
  I := ideal< P | POLS >;
  Groebner(I);
  B := Basis(I);
  dd := [ TotalDegree(B[j]) : j in [1..#B] ];
  flag := true;
  MinGens := [ ];
  for i := 1 to #Basis(I) do
    a,b := Minimum(dd);
    if flag then
      II := ideal< P | MinGens >;
    end if;
    if B[b] notin II then
      Append(~MinGens, B[b]);
      flag := true;
    else
      flag := false;
    end if;
    Remove(~B, b);
    Remove(~dd, b);
  end for;
  variables := [ "z", "y", "x", "w", "v", "u", "t", "s" ];
  AssignNames(~P, [ variables[i] : i in [1..Rank(P)] ]);
  return MinGens, P;
end function;


// Section 6: Varieties of truncated syzygy modules

SyzygyModule := function(M, FM)
  r := Dimension(M) - Dimension(JacobsonRadical(M));
  F, _, prj := DirectSum([ FM : i in [1..r] ]);
  Hom := GHom(FM, M);
  flag := true;
  while flag do
    lf := [ Random(Hom) : i in [1..r] ];
    f := &+[ MapToMatrix(prj[i])*lf[i] : i in [1..r] ];
    if Rank(f) eq Dimension(M) then
      flag := false;
    end if;
  end while;
  OM := Kernel(f);
  return OM;
end function;