// THIS FILE CONTAINS A NUMBER OF FUNCTIONS FOR DEALING WITH 2-DIMENSIONAL LATTICE POLYTOPES AND NEWTON POLYTOPES
// FOR US, A LATTICE POLYTOPE IS AN OBJECT OF THE FORM [p1,p2,...,pn] WHERE THE pi ARE OF THE FORM <a,b> FOR INTEGERS
// a AND b. IT IS ASSUMED THAT THE pi ARE THE VERTICES OF THE LATTICE POLYTOPE, AND THAT IN THE 2-DIMENSIONAL CASE
// THEY ARE ENUMERATED COUNTERCLOCKWISE.

// NewtonPolytope OUTPUTS THE NEWTON POLYTOPE OF A BIVARIATE POLYNOMIAL
// E.G. > R<x,y> := PolynomialRing(Integers(),2);
//      > NewtonPolytope(1 + x + x^2 + y);
//      [ <0, 0>, <2, 0>, <0, 1> ]
//      > NewtonPolytope(R ! 1);
//      [ <0, 0> ]

NewtonPolytope := function(f)
  N := NewtonPolygon(f);
  V := AllVertices(N);
  V := [ < Integers() ! v[1], Integers() ! v[2] > : v in V ];
  return V;
end function;

// LatticePolytope OUTPUTS THE LATTICE POLYTOPE OBTAINED BY TAKING THE CONVEX HULL OF A BUNCH OF POINTS WITH
// INTEGER COORDINATES
// E.G. > LatticePolytope([<-1,-1>,<0,0>,<1,-1>,<1,1>]);
//      [ <-1, -1>, <1, -1>, <1, 1> ]

LatticePolytope := function(points)
  min := 0;
  for p in points do
    min := Minimum([min,p[1],p[2]]);
  end for;
  R<x,y> := PolynomialRing(Integers(),2);
  V := NewtonPolytope(&+[x^(p[1]-min)*y^(p[2]-min) : p in points]);
  return [ <v[1] + min, v[2] + min> : v in V];
end function;

// Edges OUTPUTS THE SEQUENCE OF 1-DIMENSIONAL EDGES OF THE LATTICE POLYTOPE. EACH EDGE IS REPRESENTED AS THE
// SEQUENCE OF LATTICE POINTS IT CONTAINS. WHEN THE INPUT IS 2-DIMENSIONAL, EVERYTHING IS GIVEN COUNTERCLOCKWISE 
// (BOTH THE SEQUENCE OF EDGES AND THE SEQUENCE OF LATTICE POINTS ON AN EDGE).
// E.G. > Edges([<-1,-1>,<1,-1>,<1,1>]);
//      [
//          [ <-1, -1>, <0, -1>, <1, -1> ],
//          [ <1, -1>, <1, 0>, <1, 1> ],
//          [ <1, 1>, <0, 0>, <-1, -1> ]
//      ]
//      > Edges([<-1,0>,<1,0>]);
//      [
//          [ <-1, 0>, <0, 0>, <1, 0> ]
//      ]
//      > Edges([<0,0>]);
//      []
 
Edges := function(V)           
      F := [];
      W := V cat [V[1]];
      if #V le 2 then r := #V-1; else r := #V; end if;
      for i in [1..r] do
        a := W[i+1][1] - W[i][1]; b := W[i+1][2] - W[i][2];
        g := GCD(a,b);
        F := F cat [ [ <W[i][1] + (a div g)*j, W[i][2] + (b div g)*j> : j in [0..g] ] ];
      end for;
      return F;
end function;

// NumBoundaryPoints OUTPUTS THE NUMBER OF LATTICE POINTS ON THE BOUNDARY OF THE LATTICE POLYTOPE. IN CASE 
// OF A 0-DIMENSIONAL OR 1-DIMENSIONAL INPUT, SIMPLY THE TOTAL NUMBER OF LATTICE POINTS IS RETURNED.
// BoundaryPoints OUTPUTS THE POINTS THEMSELVES.
// E.G. > NumBoundaryPoints([<-1,-1>]);
//      1
//      > BoundaryPoints([<-1,-1>]);
//      [ <-1, -1> ]

NumBoundaryPoints := function(V)
  F := Edges(V);
  if #F eq 0 then return #V;
  elif #F eq 1 then return #F[1]; 
  else return &+[#F[i] - 1 : i in [1..#F]];
  end if;
end function;

BoundaryPoints := function(V)
  F := Edges(V);
  if #F eq 1 then return F[1]; 
  elif #F eq 0 then return V; 
  else
  Boundary := [];
  for i in [1..#F] do
    for k in [2..#F[i]] do
      Boundary := Boundary cat [F[i][k]];
    end for;
  end for;
  return Boundary;
  end if;
end function; 

// FakeMod IS AN AUXILIARY FUNCTION. FakeMod(a,m) RETURNS THE NUMBER IN {1,...,m} CONGRUENT TO a MOD m.

FakeMod := function(a,m)
  f := a mod m;
  if f eq 0 then f := m; end if;
  return f;
end function;

// InnerPoints ENUMERATES THE LATTICE POINTS IN THE INTERIOR OF THE LATTICE POLYTOPE. THE INTERIOR OF A
// 0- OR 1-DIMENSIONAL LATTICE POLYTOPE IS CONSIDERED EMPTY.
// E.G. > N := LatticePolytope([<1,1>,<-2,3>,<-5,1>,<-1,2>,<2,2>]);
//      > N;
//      [ <-5, 1>, <1, 1>, <2, 2>, <-2, 3> ]
//      > InnerPoints(N);
//      [ <-3, 2>, <-2, 2>, <-1, 2>, <0, 2>, <1, 2> ]
// PolGenus GIVES THE NUMBER OF INTERIOR LATTICE POINTS

InnerPoints := function(V)
  F := Edges(V);
  if #F le 1 then return [];
  else
    r := #V;
    Min,Left := Minimum([v[2] : v in V]);
    Max := Maximum([v[2] : v in V]);
    Right := Left;
    if V[FakeMod(Left-1,r)][2] eq Min then Left := FakeMod(Left-1,r); end if;
    if V[FakeMod(Right+1,r)][2] eq Min then Right := FakeMod(Right+1,r); end if;
    Inner := []; Numerator := 1;
    for yvalue in [Min+1..Max-1] do
      if V[Left][2] lt yvalue then
        Lefter := FakeMod(Left-1,r);
        LeftInverseRico := (V[Lefter][1] - V[Left][1])/(V[Lefter][2] - V[Left][2]);
        Numerator := Maximum([Numerator,V[Lefter][2]-V[Left][2]]);
        Left := Lefter;
      end if;
      if V[Right][2] lt yvalue then
        Righter := FakeMod(Right+1,r);
        RightInverseRico := (V[Righter][1] - V[Right][1])/(V[Righter][2] - V[Right][2]);
        Numerator := Maximum([Numerator,V[Righter][2]-V[Right][2]]);
        Right := Righter;
      end if;
      for xvalue in [(V[Left][1] + Ceiling(LeftInverseRico*(yvalue-V[Left][2]) + 1/(2*Numerator)))..(V[Right][1] + Floor(RightInverseRico*(yvalue-V[Right][2]) - 1/(2*Numerator)))] do Inner := Inner cat [<xvalue,yvalue>]; end for;
    end for;
    return Inner;
  end if;
end function;

PolGenus := function(V)
  return #InnerPoints(V);
end function;

// PolytopeVolume RETURNS THE 2-DIMENSIONAL EUCLIDEAN VOLUME OF THE LATTICE POLYTOPE, USING PICK'S THEOREM
// E.G. > N := LatticePolytope([<0,0>,<1,0>,<0,1>]);
//      > PolytopeVolume(N);
//      1/2

PolytopeVolume := function(V)
  if #V le 2 then return 0; else
  I := #InnerPoints(V);
  R := NumBoundaryPoints(V);
  return I + R/2 - 1;
  end if;
end function;

// ConvexHull ACCEPTS THE SAME INPUT AS LatticePolytope, BUT RETURNS ALL OF THE LATTICE POINTS, INSTEAD OF JUST
// THE VERTICES.

ConvexHull := function(S)
  N := LatticePolytope(S);
  return (BoundaryPoints(N) cat InnerPoints(N));
end function;

// SupportingHyperplanes GIVES THE SEQUENCE OF SUPPORTING HYPERPLANE INEQUALITIES OF THE LATTICE POLYTOPE, GIVEN
// IT IS 2-DIMENSIONAL (OTHERWISE AN EMPTY SEQUENCE IS RETURNED). A HYPERPLANE INEQUALITY IS GIVEN AS
// [a,b,c], WHICH SHOULD BE READ AS aX + bY =< c WITH a,b COPRIME.
// E.G. > N := LatticePolytope([<1,0>,<0,0>,<0,1>]);
//      > N;
//      [ <0, 0>, <1, 0>, <0, 1> ]
//      > SupportingHyperplanes(N);
//      [
//          [ 0, -1, 0 ],
//          [ 1, 1, 1 ],
//          [ -1, 0, 0 ]
//      ]

SupportingHyperplanes := function(S)
  E := Edges(S);
  if #E le 1 then return []; else
    array := [];
    for e in E do
      a := e[2][2] - e[1][2]; b := e[1][1] - e[2][1]; c := e[2][1]*(e[2][2]-e[1][2]) + e[2][2]*(e[1][1]-e[2][1]);
      array cat:= [[a,b,c]];
    end for;
    return array;
  end if;
end function;

// IsInteriorPolytope RETURNS true IF THE POLYTOPE CAN BE REALIZED AS THE CONVEX HULL OF THE INTERIOR LATTICE
// POINTS OF ANOTHER LATTICE POLYTOPE, AND false ELSEWHERE. IF true, SUCH A LATTICE POLYTOPE IS RETURNED AS
// WELL, AND IF THE INPUT WAS 2-DIMENSIONAL IT CONCERNS THE UNIQUE MAXIMAL POLYTOPE WITH THAT INTERIOR

IsInteriorPolytope := function(S)
  case #S:
    when 1: 
      return true, LatticePolytope([<S[1][1]-1,S[1][2]-1>,<S[1][1]-1,S[1][2]+1>,<S[1][1]+1,S[1][2]+1>,<S[1][1]+1,S[1][2]-1>]);
    when 2:
      E := Edges(S)[1]; 
      genus := #E;
      endpoint1 := <2*E[1][1]-E[2][1],2*E[1][2]-E[2][2]>;
      endpoint2 := <2*E[genus][1]-E[genus-1][1],2*E[genus][2]-E[genus-1][2]>;
      points := [endpoint1,endpoint2];
      if E[2][1] ne E[1][1] then
        for e in E do
          points cat:= [<e[1],e[2]+1>,<e[1],e[2]-1>];
        end for;
        return true, LatticePolytope(points);
      else
        for e in E do
          points cat:= [<e[1]+1,e[2]>,<e[1]-1,e[2]>];
        end for;
        return true, LatticePolytope(points);
      end if;
    else:
      H := SupportingHyperplanes(S);
      numver := #H;
      okay := true;
      points := [];
      i := 1; 
      while okay and (i le numver) do
        j := FakeMod(i+1,numver);
        M := Matrix(Rationals(),2,2,[H[i][1],H[i][2],H[j][1],H[j][2]]);
        sols := M^(-1)* Matrix(Rationals(),2,1,[H[i][3]+1,H[j][3]+1]);
        if IsIntegral(sols[1][1]) and IsIntegral(sols[2][1]) then points cat:=[<Integers() ! sols[1][1], Integers() ! sols[2][1]>];
        else okay := false; end if;
        i +:=1;
      end while;
      // IF okay THEN WE KNOW THAT THE TRANSLATED EDGES INTERSECT IN LATTICE POINTS, BUT IT MIGHT HAPPEN THAT THE POLYTOPE
      // RESULTING FROM "CONNECTING THE DOTS" IS NOT CONVEX.
      if okay then 
        i := 1;
        while okay and (i le numver) do
          j := FakeMod(i+1,numver);
          k := FakeMod(i+2,numver);
          vector1 := <points[j][1]-points[i][1],points[j][2]-points[i][2]>;
          vector2 := <points[k][1]-points[i][1],points[k][2]-points[i][2]>;
          if vector1[1]*vector2[2] - vector1[2]*vector2[1] lt 0 then okay := false; end if;
          i +:=1;
        end while;
      end if;
      if okay then return true, LatticePolytope(points); else return false, "no polytope"; end if;
  end case; 
end function;

// AreEquivalent TAKES AS INPUT TWO LATTICE POLYTOPES S1 AND S2 AND OUTPUTS false IF THERE IS NO Z-AFFINE TRANSFORMATION
// CONNECTING BOTH. IF THERE IS, THE OUTPUT IS true AND SUCH A Z-AFFINE TRANSFORMATION IS GIVEN IN THE FORM OF AN
// INVERTIBLE 2x2 MATRIX (LINEAR PART) AND A 2x1 MATRIX (TRANSLATION PART) OVER THE INTEGERS.

AreEquivalent := function(S1,S2)
  if #S1 ne #S2 then return false;
  else
    E1 := Edges(S1);
    E2 := Edges(S2);
    if {* #e : e in E1 *} ne {* #e : e in E2 *} then return false;
    else
      case #S1:
      when 1:
        transfolinear := Matrix(Integers(),2,2,[1,0,0,1]);
        transfotranslate := Matrix(Integers(),2,1,[S2[1][1]-S1[1][1],S2[1][2]-S1[1][2]]);
        return true, transfolinear, transfotranslate;
      when 2:
        vertex1 := E1[1][1]; point1 := E1[1][2]; vector1 := <point1[1]-vertex1[1],point1[2]-vertex1[2]>;
        vertex2 := E2[1][1]; point2 := E2[1][2]; vector2 := <point2[1]-vertex2[1],point2[2]-vertex2[2]>;
        _1,a1,b1 := XGCD(vector1[1],vector1[2]);
        _2,a2,b2 := XGCD(vector2[1],vector2[2]);
        transfolinear := Matrix(Integers(),2,2,[a2,b2,-vector2[2],vector2[1]])^(-1)*Matrix(Integers(),2,2,[a1,b1,-vector1[2],vector1[1]]);
        transfotranslate := Matrix(Integers(),2,1,[vertex2[1],vertex2[2]])-transfolinear*Matrix(Integers(),2,1,[vertex1[1],vertex1[2]]);
        return true, transfolinear, transfotranslate;
      else:
        n := #S1;
        foundtransfo := false;
        i := 1;
        vertex := E1[2][1]; point1 := E1[2][2]; point2 := E1[1][#E1[1]-1];
        matrix := Matrix(Rationals(),2,2,[point1[1] - vertex[1], point2[1]-vertex[1],point1[2]-vertex[2],point2[2]-vertex[2]]);
        minv := matrix^(-1);
        vertexmatrix := Matrix(Rationals(),2,1,[vertex[1],vertex[2]]);
        det2 := Determinant(matrix)^2;
        posorientation := true;
        while (not foundtransfo) and i le n do
          candvertex := E2[i][1]; 
          candpoint1 := E2[i][2]; candpoint2 := E2[FakeMod(i-1,n)][#E2[FakeMod(i-1,n)]-1];
          if not posorientation then backup := candpoint1; candpoint1 := candpoint2; candpoint2 := backup; end if;
          candmatrix := Matrix(Rationals(),2,2,[candpoint1[1]-candvertex[1], candpoint2[1]-candvertex[1],candpoint1[2]-candvertex[2], candpoint2[2]-candvertex[2]]);
          if det2 eq Determinant(candmatrix)^2 then
            transfolinear := candmatrix*minv;
            if IsIntegral(transfolinear[1][1]) and IsIntegral(transfolinear[1][2]) and IsIntegral(transfolinear[2][1]) and IsIntegral(transfolinear[2][2]) then
              candvertexmatrix := Matrix(Rationals(),2,1,[candvertex[1],candvertex[2]]);
              transfotranslate := candvertexmatrix - transfolinear*vertexmatrix;
              transformedpolytope := [];
              for p in S1 do
                transformedpolytope cat:= [ transfolinear*Matrix(Rationals(),2,1,[p[1],p[2]]) + transfotranslate ];
              end for;
              transformedpolytope := [ <Integers()! p[1][1], Integers() ! p[2][1]> : p in transformedpolytope ];
              if LatticePolytope(transformedpolytope) eq S2 then foundtransfo := true; end if;
            end if;
          end if;
          if posorientation then posorientation := false; else posorientation := true; i +:= 1; end if;
        end while;
        if not foundtransfo then return false; else 
          transfolinear := Matrix(Integers(), 2,2,[transfolinear[1][1],transfolinear[1][2],transfolinear[2][1],transfolinear[2][2]]);
          transfotranslate := Matrix(Integers(), 2,1,[transfotranslate[1][1],transfotranslate[2][1]]);
          return true, transfolinear, transfotranslate; 
        end if;        
      end case;
    end if;
  end if;
  //  return "under construction";
end function;

// LatticeWidth(V) COMPUTES THE LATTICE WIDTH OF A POLYTOPE AND RETURNS A PRIMITIVE VECTOR THAT
// SHOULD BE MAPPED TO (1,0) TO PUT THE POLYTOPE INTO A HORIZONTAL STRIP OF THAT WIDTH

MinimalInCone := function(v1,v2,v)
  g1 := GCD(v1[1],v1[2]);
  g2 := GCD(v2[1],v2[2]);
  w1 := <v1[1] div g1, v1[2] div g1>;
  w2 := <v2[1] div g2, v2[2] div g2>;
  P := ConvexHull([<0,0>,w1,w2]);
  P := Remove(P,Index(P,<0,0>));
  min := Abs(v[1]*P[1][1] + v[2]*P[1][2]);
  vec := P[1];
  for i in [2..#P] do
    a := Abs(v[1]*P[i][1] + v[2]*P[i][2]);
    if a lt min then min := a; vec := P[i]; end if;
  end for;
  return min, vec;
end function;

// RETURN INDEX OF SMALLEST ANGLE

CompareAngles := function(c1,c2);
  scal1 := c1[1][1]*c1[2][1] + c1[1][2]*c1[2][2];
  scal2 := c2[1][1]*c2[2][1] + c2[1][2]*c2[2][2];
  if (scal1 le 0) and (scal2 ge 0) then return 2;
  else
    if (scal2 le 0) and (scal1 ge 0) then return 1;
    else
      normprodsquared1 := (c1[1][1]^2 + c1[1][2]^2)*(c1[2][1]^2 + c1[2][2]^2);
      normprodsquared2 := (c2[1][1]^2 + c2[1][2]^2)*(c2[2][1]^2 + c2[2][2]^2);
      cosanglesquared1 := scal1^2/normprodsquared1;
      cosanglesquared2 := scal2^2/normprodsquared2;
      if Sign(scal1)*cosanglesquared1 ge Sign(scal2)*cosanglesquared2 then return 1; else return 2; end if;
    end if;
  end if;
end function;

LatticeWidth := function(V)
if #V eq 0 then return -1;
else if #V le 2 then return 0;
     else
       min := Maximum([v[2] : v in V]) - Minimum([v[2] : v in V])+1;
       vec1 := <1,0>;
       n := #V;
       step := 1;
       main := 1;
       mainstart := <V[1][1] - V[n][1], V[1][2] - V[n][2]>;
       // FIND OPPOSITE VERTEX
       S := SupportingHyperplanes(V);
       _,opp := Minimum([S[n][1]*v[1] + S[n][2]*v[2] : v in V]);
       while step le 2*n do
           oppend := <V[FakeMod(opp+1,n)][1] - V[opp][1], V[FakeMod(opp+1,n)][2] - V[opp][2]>;
	   if oppend[1]*mainstart[2] - oppend[2]*mainstart[1] eq 0 then
	     opp := FakeMod(opp + 1,n);
	   else
             mainend := <V[FakeMod(main+1,n)][1] - V[main][1], V[FakeMod(main+1,n)][2] - V[main][2]>;
             oppstart := <-mainstart[1],-mainstart[2]>;
	     if CompareAngles([mainstart,mainend],[oppstart,oppend]) eq 1 then
	       // print main, opp, "Cone 1:", mainstart, mainend, " Cone 2:", oppstart, oppend, " Cone 1 smaller than Cone 2";
               candidatemin,vec := MinimalInCone(mainstart,mainend,<-V[opp][2]+V[main][2],V[opp][1]-V[main][1]>);
	       mainnew := opp;
	       mainstart := <-mainend[1],-mainend[2]>;
	       opp := FakeMod(main + 1,n);
	       main := mainnew;
	     else
	       // print main, opp, "Cone 1:", mainstart, mainend, " Cone 2:", oppstart, oppend, " Cone 2 smaller than Cone 1";
               candidatemin,vec := MinimalInCone(oppstart,oppend,<-V[main][2]+V[opp][2],V[main][1]-V[opp][1]>);
               mainstart := <-oppend[1],-oppend[2]>;
	       opp := FakeMod(opp + 1,n);
	     end if;
	     if candidatemin lt min then min := candidatemin; vec1 := vec; end if;
	   end if;
	 step +:=1;
       end while;
     end if;
end if;
return min, vec1;
end function;


// EnumerateHyperellipticPolytopes LISTS, FOR INPUT g AT LEAST 2, POLYTOPES REPRESENTING ALL EQUIVALENCE
// CLASSES OF POLYTOPES HAVING g INTERIOR LATTICE POINTS ON A LINE. THE CONSTRUCTION FOLLOWS
// A PROOF OF THE NUMBER OF SUCH EQUIVALENCE CLASSES DUE TO KOELMAN.

EnumerateHyperellipticPolytopes := function(g)
if g lt 2 then return "Genus should be at least 2...";
else
  L := [];
  // TYPE 0
  for k in [g..2*g] do L cat:= [LatticePolytope([<0,-1>,<1,1>,<1+2*g-k,1>,<k,-1>])]; end for;
  // TYPE 1
  for k in [0..g] do
    for u in [0..k] do
      L cat:= [LatticePolytope([<0,-1>,<1,1>,<1+u,1>,<g+1,0>,<k,-1>])];
    end for;
  end for;
  for k in [(g+1)..(2*g+1)] do
    for u in [0..(2*g-k+1)] do
      L cat:= [LatticePolytope([<0,-1>,<1,1>,<1+u,1>,<g+1,0>,<k,-1>])];
    end for;
  end for;
  // TYPE 2
  for s in [0..(g+1)] do
    for k in [0..(g+1-s)] do
      for u in [0..k] do
        L cat:= [LatticePolytope([<0,0>,<0,-1>,<s,1>,<k,-1>,<g+1,0>,<u+s,1>])];
      end for;
    end for;
    for k in [(g+2-s)..(2*g+2-2*s)] do
      for u in [0..(2*g-k-2*s+2)] do
        L cat:= [LatticePolytope([<0,0>,<0,-1>,<s,1>,<k,-1>,<g+1,0>,<u+s,1>])];
      end for;
    end for;
  end for;
  return L;
end if;
end function;

// EnumerateGenusZeroPolytopes LISTS, FOR INPUT R AT LEAST 3, POLYTOPES REPRESENTING ALL EQUIVALENCE CLASSES
// OF 2-DIMENSIONAL LATTICE POLYTOPES HAVING NO INTERIOR LATTICE POINTS AND HAVING R LATTICE POINTS ON THE BOUNDARY.
// E.G. > EnumerateGenusZeroPolytopes(6);
//      [
//          [ <0, 0>, <2, 0>, <0, 2> ],
//          [ <0, 0>, <4, 0>, <0, 1> ],
//          [ <0, 0>, <3, 0>, <1, 1>, <0, 1> ],
//          [ <0, 0>, <2, 0>, <2, 1>, <0, 1> ]
//      ]

EnumerateGenusZeroPolytopes := function(R)
  if R eq 6 then list := [[<0,0>,<2,0>,<0,2>]]; else list := []; end if;
  for i in [1..(R div 2)] do
    list cat:=[LatticePolytope([<0,0>,<R-i-1,0>,<i-1,1>,<0,1>])];
  end for;
  return list;
end function;

EnumeratePolytopesWithInterior := function(V)
  // CURRENTLY ONLY COVERS THE NON-HYPERELLIPTIC CASE -- TO BE CONTINUED
  isinterior,maxpol := IsInteriorPolytope(V);
  if isinterior then
    genus := PolGenus(maxpol);
    E := Edges(maxpol); n := #E;
    listofsubsets := [];
    for e in E do
      ee := SequenceToSet(Remove(e,1));
      listofsubsets cat:=[Subsets(ee,1) join Subsets(ee,2) join {{e[1]}}];
    end for;
    firstlist := [];
    for c in CartesianProduct(listofsubsets) do
      s := &join[c[i] : i in [1..n]];
      L := LatticePolytope(s);
      if PolGenus(L) eq genus then firstlist cat:= [L]; end if;
    end for;
    cleanedlist := [firstlist[1]];
    for f in firstlist do
      add := true;
      i := 1; n := #cleanedlist;
      while add and i le n do
        if AreEquivalent(f,cleanedlist[i]) then add:=false; end if;
        i +:=1;
      end while;
      if add then cleanedlist cat:= [f]; end if;
    end for;
    return cleanedlist;
  else return [];
  end if;
end function;

EnumeratePolytopesOfGenus := function(g);
  tim := Cputime();
  //GENUS 1 BY HAND
  genuslist := [ [ [], [], [ [ <0, 0>, <2, 1>, <1, 2> ] ],
  [ [ <0, 0>, <2, 0>, <1, 2> ], [ <0, 1>, <1, 0>, <2, 0>, <1, 2> ],
    [ <0, 1>, <1, 0>, <2, 1>, <1, 2> ] ],
  [ [ <0, 0>, <2, 0>, <2, 1>, <1, 2> ],
    [ <0, 0>, <1, 0>, <2, 1>, <1, 2>, <0, 1> ] ],
  [ [ <0, 0>, <3, 0>, <0, 2> ], [ <0, 0>, <2, 0>, <1, 2>, <0, 2> ],
    [ <0, 0>, <2, 0>, <2, 1>, <1, 2>, <0, 1> ],
    [ <0, 1>, <1, 0>, <2, 0>, <2, 1>, <1, 2>, <0, 2> ] ],
  [ [ <0, 0>, <3, 0>, <2, 1>, <0, 2> ],
    [ <0, 0>, <2, 0>, <2, 1>, <1, 2>, <0, 2> ] ],
  [ [ <0, 0>, <4, 0>, <0, 2> ],
    [ <0, 0>, <3, 0>, <1, 2>, <0, 2> ],
    [ <0, 0>, <2, 0>, <2, 2>, <0, 2> ] ],
  [ [ <0, 0>, <3, 0>, <0, 3> ] ]
    ] ];
  print "Genus 1 done, 16 found. Time needed so far:", Cputime(tim);
  // GENUS 2 QUASI BY HAND
  tobetreated := EnumerateHyperellipticPolytopes(2);
  list := [ [], [], [], [], [], [], [], [], [], [] ];
  for V in tobetreated do
    list[NumBoundaryPoints(V)] cat:= [V];
  end for;
  print "Genus 2 done, 45 found. Time needed so far:", Cputime(tim);
  genuslist[2] := list;
  for i in [3..g] do
  list := []; for j in [1..(2*i+6)] do list cat:=[ [] ]; end for;  
    // HYPERELLIPTICS
    tobetreated := EnumerateHyperellipticPolytopes(i);
    // GENUS ZERO INTERIOR
    for V in EnumerateGenusZeroPolytopes(i) do tobetreated cat:=EnumeratePolytopesWithInterior(V); end for;
    // LOWER GENUS INTERIOR
    for j in [1..(i-1)] do
      if #genuslist[j] ge i-j then
        for V in genuslist[j][i-j] do tobetreated cat:= EnumeratePolytopesWithInterior(V); end for;
      end if;
    end for;
  print "Genus",i,"done,", #tobetreated, "found. Time needed so far:", Cputime(tim);
  for V in tobetreated do
    list[NumBoundaryPoints(V)] cat:= [V];
  end for;
  genuslist[i] := list;
  end for;
  // FINALIZE TO EASE COUNTING
  genlist := [];
  for i in [1..g] do
    genlist[i] := &cat(genuslist[i]);
  end for;
  return genlist;
end function;
