// FILENAME: canonical.m
// THIS ACCOMPANIES THE ARTICLE
// W. Castryck, F. Cools, "A minimal set of generators for the canonical ideal of a non-degenerate curve", preprint (2014)

//*********************************************************************************************************************

// THE MAIN FUNCTIONS IN THE CURRENT FILE ARE TorIdeal() AND NondegIdeal(). THE FIRST FUNCTION TAKES AS INPUT A LATTICE
// POLYGON AND A FIELD k, AND OUTPUTS A MINIMAL SET OF GENERATORS FOR THE CORRESPONDING TORIC IDEAL. THE SECOND
// FUNCTION TAKES AS INPUT A BIVARIATE POLYNOMIAL f TAKING COEFFICIENTS IN SOME FIELD k; IT IS ASSUMED TO BE 
// NON-DEGENERATE WITH RESPECT TO ITS NEWTON POLYGON (*). THE OUTPUT IS THEN A MINIMAL SET OF GENERATORS FOR THE 
// CANONICAL IDEAL OF THE FUNCTION FIELD DEFINED BY f.

// WE RELY ON A NUMBER OF FUNCTIONS FOR DEALING WITH LATTICE POLYGONS, THAT WERE WRITTEN AS PART OF THE PAPERS
// W. Castryck, "Moving out the edges of a lattice polygon", Disc. Comp. Geom. 47(3), pp. 496-518 (2012)
// W. Castryck, F. Cools, "Linear pencils encoded in the Newton polygon", preprint (2014) 
// THE FILE basic_commands.m THAT ACCOMPANIES THE LATTER PAPER IS LOADED BY THE COMMAND BELOW.

load "basic_commands.m";

// (*) IT IS IN FACT SUFFICIENT THAT THE ZERO LOCUS OF f INSIDE (k*)^2 COMPACTIFIES TO A SMOOTH CURVE IN THE TORIC 
// SURFACE ASSOCIATED TO THE NEWTON POLYGON OF f.
//
// IN BOTH CASES THE OUTPUT CONSISTS OF A SEQUENCE OF POLYNOMIALS IN A NUMBER OF VARIABLES. IN THE CASE OF TorIdeal()
// THE VARIABLES ARE INDEXED BY THE LATTICE POINTS IN THE INPUT POLYTOPE P. THE ORDER OF THE VARIABLES CORRESPONDS TO
// THE ORDER OF THE OUTPUT OF Sort(ConvexHull(P)). IN THE CASE OF NondegIdeal(), THE VARIABLES ARE INDEXED BY THE 
// LATTICE POINTS IN THE INTERIOR OF THE NEWTON POLYTOPE OF f. HERE THE ORDER OF THE VARIABLES CORRESPONDS TO THE ORDER
// OF THE OUTPUT OF Sort(InnerPoints(NewtonPolytope(f))).

//*********************************************************************************************************************

//* PairToRank() AND RankToPair() ARE FUNCTIONS MADE FOR ENUMERATING THE ENTRIES OF THE TRIANGLE
//*                                                       (N  , N  )
//*                                            (N-1, N-1) (N  , N-1)
//*                                 (N-2, N-2) (N-1, N-2) (N  , N-2)
//*                              .       .          .         .
//*                           .          .          .         .
//*                        .             .          .         .
//*            (2  , 2  )  .  .  .  (N-2, 2  ) (N-1, 2  ) (N  , 2  )
//* (1  , 1  ) (2  , 1  )  .  .  .  (N-2, 1  ) (N-1, 1  ) (N  , 1  )
//* BOTTOM-UP AND LEFT-TO-RIGHT.
//*
//* EXAMPLE:
//* > PairToRank(2,2,3);
//* 4
//* > RankToPair(4,3);
//* [ 2, 2]

function PairToRank(i1,i2,N)
 if i1 le i2 then a := i1; b := i2; else a := i2; b := i1; end if;
 return (a-1)*N + (b-a) + 1 - ((a-1)*(a-2) div 2);
end function;

function RankToPair(r,N)
 i1 := 0;
 i2 := r;
 while i2 gt (N - i1) do
  i2 := i2 - (N - i1);
  i1 +:= 1;
 end while;
 return [i1 + 1,i1+i2];
end function;

//* TorIdealLinAlg ACCEPTS AS INPUT A LATTICE POLYGON P (GIVEN AS A SEQUENCE OF VERTICES) AND A FIELD field, AND
//* RETURNS A LINEARLY INDEPENDENT GENERATING SET FOR THE SPACE OF QUADRATIC RELATIONS IN THE IDEAL OF THE TORIC
//* SURFACE OVER field ASSOCIATED TO P. EVERY ELEMENT OF THIS GENERATING SET IS A BINOMIAL ENCODED AS A COORDINATE 
//* VECTOR (0, 0, ..., 0, 1, 0, ..., 0, -1, 0, ..., 0, 0) WITH RESPECT TO THE STANDARD BASIS OF THE SPACE OF 
//* QUADRATIC POLYNOMIALS. THE FUNCTION ALSO OUTPUTS THE PARENT VECTOR SPACE.
//*
//* REMARK: THE QUADRATIC RELATIONS DO NOT NECESSARILY GENERATE THE WHOLE IDEAL. THIS IS TRUE IF AND ONLY IF THE
//* NUMBER OF LATTICE POINTS ON THE BOUNDARY OF P IS AT LEAST FOUR.
//*
//* EXAMPLE:
//* > P := [ <0, 0>, <2, 0>, <1, 1>, <0, 1> ];
//* > TorIdealLinAlg( P, Rationals() );
//* [
//*    ( 1  0  0  0  0  0  0  0 -1  0  0  0  0  0  0),
//*    ( 0  0  1  0  0  0  0 -1  0  0  0  0  0  0  0),
//*    ( 0  0  0  1  0  0  0  0  0  0  0 -1  0  0  0)
//* ]
//* Full Vector space of degree 15 over Rational Field


function TorIdealLinAlg(P,field)
  points := Sort(ConvexHull(P));
  N := #points;
  dim := N*(N+1) div 2;
  V := VectorSpace(field, dim);
  gens := [];
  for i1 in [1..N] do for i2 in [i1..N] do for i3 in [1..N] do for i4 in [(i3+1)..N] do
    if points[i1] + points[i2] eq points[i3] + points[i4] then 
      vec := V ! 0;
      pos1 := PairToRank(i1,i2,N); pos2 := PairToRank(i3,i4,N);
      if pos1 ne pos2 then vec[PairToRank(i1,i2,N)] := 1; vec[PairToRank(i3,i4,N)] := -1; end if;
      gens cat:=[vec];
    end if;
  end for; end for; end for; end for;
  U := sub<V | gens>;
  return Basis(U),V;
end function;

//* TorIdeal INVOKES TorIdealLinAlg FOR THE QUADRATIC RELATIONS, AND TRANSLATES THE OUTPUT OF THE LATTER
//* INTO ACTUAL POLYNOMIALS. THE TRANSLATION IS DONE BY THE AUXILIARY FUNCTION IdealFromLinAlg. IT THEN
//* MANUALLY ADDS CUBIC RELATIONS IF NEEDED, AS DESCRIBED IN OUR PAPER. THE FUNCTION FindVertexTriple IS
//* A TOOL FOR FINDING SUCH CUBIC RELATIONS.
//*
//* EXAMPLE:
//* > P := [ <0, 0>, <2, 0>, <1, 1>, <0, 1> ];
//* > TorIdeal(P, Rationals());
//* [
//*    $.1^2 - $.2*$.5,
//*    $.1*$.3 - $.2*$.4,
//*    $.1*$.4 - $.3*$.5
//* ]
//* > P := [ <0, 0>, <3, -1>, <2, 1> ];
//* > TorIdeal(P, Rationals());
//* [
//*    -$.4^2 + $.3*$.5,
//*    $.1*$.2*$.3 - $.4*$.5^2,
//*    $.1*$.2*$.4 - $.5^3
//* ]


function IdealFromLinAlg(bas,V)
  field := BaseField(V);
  dim := Dimension(V);
  N := Floor(Sqrt(2*dim));
  R := PolynomialRing(field,N,"grevlex");
  id := [];
  for b in bas do
    f := &+[ b[i]*R.(RankToPair(i,N)[1])*R.(RankToPair(i,N)[2]) : i in [1..dim]];
    id cat:= [f];
  end for;
  return id,R;
end function;

function FindVertexTriple(P,P1)
  vx := P[1][1] + P[2][1] + P[3][1]; vy := P[1][2] + P[2][2] + P[3][2];
  for triple in CartesianPower(P1,3) do
    if (vx eq triple[1][1] + triple[2][1] + triple[3][1]) and (vy eq triple[1][2] + triple[2][2] + triple[3][2]) then
      goodtriple := triple; break;
    end if;
  end for;
  return goodtriple;
end function;

function TorIdeal(P,field)
  bas,V := TorIdealLinAlg(P,field);
  I,R := IdealFromLinAlg(bas,V);
  points := Sort(ConvexHull(P));
  P1 := InnerPoints(P); g := #P1;
  if (#points - g) eq 3 and (g ge 1) then
    triple := FindVertexTriple(P,P1);
    I cat:= [ R.Index(points,P[1])*R.Index(points,P[2])*R.Index(points,P[3]) - 
              R.Index(points,triple[1])*R.Index(points,triple[2])*R.Index(points,triple[3]) ];
    // EXTRA RELATIONS IN THE HYPERELLIPTIC CASE
    if #LatticePolytope(P1) eq 2 then
      _,matrix,transl := AreEquivalent([<0,-1>,<g,0>,<-1,1>],P);
      vertex1 := < -matrix[1][2] + transl[1][1], -matrix[2][2] + transl[2][1]>;
      vertex2 := < matrix[1][2] - matrix[1][1] + transl[1][1], matrix[2][2] - matrix[2][1] + transl[2][1]>;
      origin := <transl[1][1],transl[2][1]>;
      for i in [1..(g-1)] do
        vertex3 := < (g-i)*matrix[1][1] + transl[1][1], (g-i)*matrix[2][1] + transl[2][1] >;
        otherpoint := < (g-i-1)*matrix[1][1] + transl[1][1], (g-i-1)*matrix[2][1] + transl[2][1] >;
        I cat:= [R.Index(points,vertex1)*R.Index(points,vertex2)*R.Index(points,vertex3) 
                 - R.Index(points,origin)^2*R.Index(points,otherpoint) ];
      end for;
    end if;
  end if;
  return I;
end function;

//* NondegIdealLinAlg TAKES AS INPUT A BIVARIATE LAURENT POLYNOMIAL f TAKING COEFFICIENTS IN ANY FIELD SUCH THAT
//* IT IS NON-DEGENERATE WITH RESPECT TO ITS NEWTON POLYGON P. WE ASSUME THAT THE INTERIOR P1 OF THE LATTER 
//* CONTAINS AT LEAST 3 LATTICE POINTS (I.E. THE GENUS OF THE CURVE DEFINED BY f IS AT LEAST 3). THE OUTPUT IS OF 
//* THE SAME TYPE AS IN TorIdealLinAlg, I.E. IT IS A LINEARLY INDEPENDENT GENERATING SET FOR THE SPACE OF QUADRATIC 
//* RELATIONS IN THE CANONICAL IDEAL THE CURVE DEFINED BY f, EACH OF WHICH IS ENCODED AS A COORDINATE VECTOR
//*
//* EXAMPLE:
//* > R<x,y> := PolynomialRing(Rationals(),2);
//* > f := 2*x^2 - 3*x^4*y^2 + 5*x^2*y^4 - 7*y^2;
//* > NondegIdealLinAlg(f);
//* [
//*    ( 0  0  1  0  0  0  0  0  0  0  0  0  0  0 -1),
//*    ( 0  0  0  0  0  0  0  1  0  0  0  0  0  0 -1),
//*    ( 2  0  0  0  0 -7  0  0  0  0  0  0 -3  0  5)
//* ]

function NondegIdealLinAlg(f);
  field := BaseRing(Parent(f));
  P := NewtonPolytope(f);
  points1 := Sort(InnerPoints(P));
  P1 := LatticePolytope(points1);
  bas := TorIdealLinAlg(P1,field);
  N := #points1;
  dim := N*(N+1) div 2;
  VSp := VectorSpace(field,dim);
  nulvec := VSp ! 0;


  points2 := InnerPoints(P1);
  P2 := LatticePolytope(points2); 

  T := Terms(f); C := Coefficients(f); numT := #T;  
  for basepoint in points2 do
    vec := nulvec;
    indexbasepoint := Index(points1,basepoint);
    for i in [1..numT] do
      exp := Exponents(T[i]); exp := <exp[1],exp[2]>;
      if exp in points1 then
        indexexp := Index(points1,exp);
        vec[PairToRank(indexexp,indexbasepoint,N)] := C[i];
      else
        for u in points1 do for v in points1 do
          if <u[1] - basepoint[1] + v[1] - basepoint[1], u[2] - basepoint[2] + v[2] - basepoint[2]> 
             eq <exp[1] - basepoint[1],exp[2] - basepoint[2]> 
             then U := u; V := v; break u; 
          end if;
        end for; end for;
        indexu := Index(points1,U);
        indexv := Index(points1,V);
        vec[PairToRank(indexu,indexv,N)] := C[i];
      end if;
    end for;
    bas cat:= [vec];
  end for;   
  return bas,VSp;
end function;

//* NondegIdeal DOES THE SAME, BUT NOW THE ONLY REQUIREMENT ON P1 IS THAT IT CONTAINS AT LEAST 3 LATTICE POINTS (I.E.
//* THE GENUS OF OUR CURVE IS AT LEAST 3), AND AGAIN OUTPUT IS NO LONGER ENCODED AS A SEQUENCE OF VECTORS, BUT REALLY AS 
//* A SEQUENCE OF HOMOGENEOUS POLYNOMIALS.
//*
//* EXAMPLE:
//* > P := [ <0, 0>, <2, 0>, <1, 1>, <0, 1> ];
//* > TorIdeal( P, Rationals() );
//* [
//*    $.1^2 - $.2*$.5,
//*    $.1*$.3 - $.2*$.4,
//*    $.1*$.4 - $.3*$.5
//* ]
//*
//* (THE ACTUAL TRANSLATION IS DONE BY THE AUXILIARY FUNCTION IdealFromLinAlg).


function NondegIdeal(f)
  bas,V := NondegIdealLinAlg(f);
  I,R := IdealFromLinAlg(bas,V);
  P := NewtonPolytope(f);
  points1 := Sort(InnerPoints(P));
  P1 := LatticePolytope(points1);  
  g := #points1;
  points2 := InnerPoints(P1);

  if #P1 gt 2 then

    // ADDING THE CUBIC GENERATOR OF THE TORIC SURFACE IF NEEDED

    if #points2 eq 1 and g eq 4 then
      I cat:= [ R.Index(points1,P1[1])*R.Index(points1,P1[2])*R.Index(points1,P1[3]) - R.Index(points1,points2[1])^3 ];
    end if;

    if #points2 eq 0 then

      T := Terms(f); C := Coefficients(f); numT := #T;

      // ADDING THE QUARTIC GENERATOR IN THE PLANE QUARTIC CASE
      if g eq 3 then
        _, matrix, transl := AreEquivalent( P1, [<1,1>,<2,1>,<1,2>] );
        polynomialtoadd := R ! 0;
        for i in [1..numT] do
          exp := Exponents(T[i]);
          newexp := < matrix[1][1]*exp[1] + matrix[1][2]*exp[2] + transl[1][1], matrix[2][1]*exp[1] + matrix[2][2]*exp[2] + transl[2][1] >;
          polynomialtoadd +:= C[i]*R.Index(points1,P1[3])^newexp[1]*R.Index(points1,P1[1])^newexp[2]*R.Index(points1,P1[2])^(4 - newexp[1] - newexp[2]);
        end for;
        I cat:= [polynomialtoadd];
      end if;

      // ADDING THE CUBIC GENERATORS IN THE CLIFFORD INDEX 1 CASE
      if g gt 3 then
        P1double := [];
        for point in P1 do P1double cat:= [ <2*point[1],2*point[2]> ]; end for;
        W2 := InnerPoints(P1double);
        for basepoint in W2 do
          polynomialtoadd := R ! 0;
          for i in [1..numT] do
            exp := Exponents(T[i]); exp := <exp[1], exp[2]>;
            for v1 in points1 do for v2 in points1 do for v3 in points1 do
              if exp eq <v1[1] + v2[1] + v3[1] - basepoint[1], v1[2] + v2[2] + v3[2] - basepoint[2]> then
                polynomialtoadd +:= C[i]*R.Index(points1,v1)*R.Index(points1,v2)*R.Index(points1,v3);
                break v1;
              end if;
            end for; end for; end for;
          end for;
          I cat:= [polynomialtoadd];
        end for;  
      end if;
    end if;

  end if;

  return I;
end function;
