needsPackage "Polyhedra"
---------------------------------------------------------------------------
-- PURPOSE : Computations of implicit equations of toric rational surfaces embedded in a projective space, by means of approximation complexes 
-- PROGRAMMER : Nicolás Botbol, Marc Dohm, Manuel Dubinsky
-- UPDATE HISTORY : August 2010
---------------------------------------------------------------------------
newPackage("Implicitization",
   Version => "1.2.2",
   Date => "24 August 2010",
   Authors => {
	 {  Name => "Nicolás Botbol, Marc Dohm, Manuel Dubinsky",
	    Email => "nbotbol@dm.uba.ar",
	    HomePage => "http://mate.dm.uba.ar/~nbotbol/"}
         },
   Headline => "Implicit equations of toric rational surfaces embedded in a projective space, by means of approximation complexes",
   DebuggingMode => true
   )

export {
	"isGoodDegree",		-- tested
	"degreeImplicitEq", 	-- tested
	"ToricEmbedding",
	"newToricEmbedding", 	-- tested
	"teToricRing", 		-- tested
	"getPolytope",
	"teToricRationalMap", 	-- tested
	"representationMatrix", -- tested
	"implicitEq",		-- tested
	"polynomialsToPolytope",-- tested
	"relationNandP", 		-- tested
	"MaxMinor"
}

needsPackage "Polyhedra"

---------------------------------------------------------------------------
-- DEFINING NEW TYPES
---------------------------------------------------------------------------

-- Defining the new type ToricEmbedding
ToricEmbedding = new Type of MutableHashTable
globalAssignment ToricEmbedding

-- Modifying the standard output for a ToricEmbedding to give an overview of its characteristica
net ToricEmbedding := te -> ( horizontalJoin flatten ( 
	  "{", 
	  -- prints the parts vertically
	  stack (horizontalJoin \ sort apply({userRing,
							tePolyhedron,
							teLattice,
							teLatticeLen,
							teRing}, key -> (net key, " => ", net te#key))),
	  "}" ))

---------------------------------------------------------------
---------------------------------------------------------------
-- PURPOSE : Auxiliary method that computes a new list 'l3' with only the elements of 'l1' that are not on 'l2'
diffLists = method()

--   INPUT : '(l1, l2)', two lists 'l1' and 'l2' of elements of the same type
--  OUTPUT : 'l3' a new list with only the elements of 'l1' that are not on 'l2' 
diffLists (List, List) := (l1, l2) -> (
	 return select(l1, i -> toString(position(l2, j -> toString(j) == toString(i) )) == "null" )	
)

---------------------------------------------------------------
-- FUNCTIONS TO CONSTRUCT TORIC EMBEDDINGS
---------------------------------------------------------------

-- PURPOSE : Method that build the toric embedding of the affine space of dimension 2 from the ample divisor given by a ploytope
newToricEmbedding = method()

--   INPUT : '(aPolyhedron, aRing)', a two-dimensional polyhedron 'aPolyhedron' and the two-variable PolynomialRing 
--  OUTPUT : A 'ToricEmbedding' corresponding to the toric embedding of spec(aRing) given by the ample divisor associated to 'aPolyhedron'
newToricEmbedding (Polyhedron, Ring) := (aPolyhedron, aRing) -> (
	lat := latticePoints aPolyhedron;
	latLen := length(lat) - 1; -- for indexing
	ringSymbols := diffLists({u,v,w,x,y},{aRing_0,aRing_1});
	new ToricEmbedding from {
		userRing => aRing,
		tePolyhedron => aPolyhedron,
		teLattice => lat,
		teLatticeLen => latLen,
		teRing => QQ[aRing_0,aRing_1,ringSymbols_0,ringSymbols_1,ringSymbols_2,T_0..T_latLen,MonomialOrder=>Eliminate 5]
	}
)

--   INPUT : 'polynomialList', a list of polynomials
--  OUTPUT : A 'ToricEmbedding' corresponding to the toric embedding of spec(aRing) given by the ample divisor associated to 'aPolyhedron'
newToricEmbedding (List) := (polynomialList) -> (
	newToricEmbedding(polynomialsToPolytope (polynomialList), ring polynomialList_0)
)

---------------------------------------------------------------------------
---------------------------------------------------------------------------
-- PURPOSE : Return the convexHull of the union of the Newton polytope of all the polynomial in the list
polynomialsToPolytope = method();

--   INPUT : 'polynomialList' a list of polynomials
--  OUTPUT : 'N' the Newton polytope of all polynomials in 'polynomialList'
polynomialsToPolytope (List) := (polynomialList) -> (
	aPolytope := newtonPolytope polynomialList_0;
	scan(polynomialList, i -> aPolytope = convexHull(aPolytope, newtonPolytope i));
	return aPolytope;
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------
-- PURPOSE : Auxiliary method that computes the maximal 'x' and 'y' coordinate of the Newton polytope associated to the toric embedding
teMaxDegree = method()

--   INPUT : '(te, indexVar)' a toric embedding 'te' and 0 for the 'x' coordinate and 1 for the 'y' coordinate
--  OUTPUT : 'ZZ'  an integer with the maximum degree on the corresponding coordinate 'x'or 'y'
teMaxDegree (ToricEmbedding, ZZ) := (te, indexVar) -> (
 	max apply(0..te#teLatticeLen, i -> (te#teLattice_i_0_indexVar))
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------
-- PURPOSE : Auxiliary method that bihomogenize all the monomials with exponents in the newton polytope of 'te'.
teBigradedMonomials = method()

--   INPUT : '(te)' a toric embedding.
--  OUTPUT : 'bigradedMonomials' is a list that contains the biomogenization of all the monomials with exponents in the newton polytope of 'te'.
teBigradedMonomials (ToricEmbedding) := (te) -> (
	use te#teRing;
	maxDegreeS:= teMaxDegree(te, 0); -- max degree of all polynomials in S_0
	maxDegreeT:= teMaxDegree(te, 1); -- idem S_1	LPP:= latticePoints P;
	bigradedMonomials := {};
	for i from 0 to te#teLatticeLen do (
		bigradedMonomials = bigradedMonomials | {((te#teRing_0)^(te#teLattice_i_0_0)) * ((te#teRing_1)^(te#teLattice_i_0_1)) * ((te#teRing_2)^(maxDegreeS - (te#teLattice_i_0_0))) * ((te#teRing_3)^(maxDegreeT - (te#teLattice_i_0_1)))};
	);
	bigradedMonomials
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------
-- PURPOSE : Method that computes the toric coordinate ring from a given polytope
teToricRing = method()

--   INPUT : 'te' a Toric Embedding that contains the information of the toric embedding
--  OUTPUT : 'A' a QuotienRing which is  the toric coordinate ring of 'te'
teToricRing (ToricEmbedding) := (te) -> (
	if ( not te#?teToricRing) then {
		use te#teRing;
		st := teBigradedMonomials(te);
		prodT := 1;
		scan((0..te#teLatticeLen), i -> prodT=prodT*T_i); -- Product of all Tt_i
		J:=ideal(1 - (te#teRing_4) * prodT);
		scan((0..te#teLatticeLen), i -> J=J+ideal (T_i - st_i)); -- ?
		J= selectInSubring(1,gens gb J);
		QofT:=QQ[T_0..T_(te#teLatticeLen)];
		J=sub(J,QofT);
		te#teToricRing = QofT/ideal(J);
	};

	return te#teToricRing;
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------
-- PURPOSE : Auxiliary method that gets the polytope from a ToricEmbedding
getPolytope = method();

--   INPUT : 'te' a ToricEmbedding
--  OUTPUT : 'P' the polytope of the embedding 'te'
getPolytope(ToricEmbedding) := (te) -> (
	return te#tePolyhedron;
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------
-- PURPOSE : Method that rewrites the lattice points of N in terms of some homothety of a different polytope P
-- COMMENT : see Pennies, Nickels, Dimes and Quarters (Computations in algebraic geometry with Macaulay 2 Editors: D. Eisenbud, D. Grayson,M. Stillman, and B. Sturmfels)
relationNandP = method();

--   INPUT : '(N,P)' two polytopes 'N' and 'P'. 'P' is spected to be smaller, but in such a way that some homothety of 'P' covers 'N'.
--  OUTPUT : 'Mat' a matrix that rewrites the lattice points of 'N' in terms of some homothety of 'P'
relationNandP(Polyhedron, Polyhedron) := (N,P) -> (
	latticeP := latticePoints P;
	latticePLen = length latticeP;

	ones := {}; scan((0..latticePLen - 1), i -> ones = ones | {1});

	R := QQ[T_0..T_(latticePLen - 1), Degrees => transpose ({ones} | matrixToList matrix{latticeP})];

	homothetyNAndP := homothety(N,P);

	latticeN := latticePoints N;
	coordinatesNinTs := basis({homothetyNAndP,latticeN_0_0_0,latticeN_0_0_1}, R);
	
	latticeNLen = length latticeN;

	for i from 1 to latticeNLen - 1 do (
		coordinatesNinTs = coordinatesNinTs || matrix {{(basis({homothetyNAndP,latticeN_i_0_0,latticeN_i_0_1}, R))_0_0}};
	);

	return coordinatesNinTs;
); 

---------------------------------------------------------------------------
---------------------------------------------------------------------------
-- PURPOSE : Auxiliary method that computs the smallest integer 'h' such that 'N' is included in 'h.P' 
homothety = method();

--   INPUT : '(N,P)' two polytopes 'N' and 'P'. 'P' is spected to be smaller, but in such a way that some homothety of 'P' covers 'N'.
--  OUTPUT : 'h' is the smallest integer such that 'N' is included in 'h.P' 
homothety (Polyhedron, Polyhedron) := (N,P) -> (
	latticeP := latticePoints P;
	h := 1;
	while not contains(convexHull(h * latticeP), N) do (
		h = h + 1;
	);
	return h;
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------
-- PURPOSE : Computes the rational map G defined over the toric coordinate ring given by {f0,f1,f2,f3}
teToricRationalMap = method();

--   INPUT : '{f0,f1,f2,f3}' a list of polynomials in two variables
--  OUTPUT : 'G' is a 1x4-matrix that represents the rational map G defined over the toric coordinate ring 'te#teRing' of the toric variety associated to the Newton polytope of '{f0,f1,f2,f3}'
teToricRationalMap(List) := (polynomialList) -> (
	te := newToricEmbedding(polynomialList);

	use te#teRing;
	torVar := teToricRing(te);
	-- substitute monomials by variables Ti interpreted in the ring te#teRing
	(M,C) := coefficients( sub( matrix {polynomialList}, te#teRing), Variables=>{te#teRing_0, te#teRing_1}, Monomials=>teBigradedMonomials( te ) );
	Ts := sub(matrix {gens torVar},te#teRing);
	G := Ts*C; 
	G = matrix{{G_(0,0),G_(0,1),G_(0,2),G_(0,3)}};

	return sub (G,torVar);
);

--   INPUT : '({f0,f1,f2,f3}, P)' where '{f0,f1,f2,f3}' a list of polynomials in two variables and P is a Polyhedron.
--  OUTPUT : 'G' is a 1x4-matrix that represents the rational map G defined over the toric coordinate ring 'te#teRing' of the toric variety associated to P
teToricRationalMap(List, Polyhedron) := (polynomialList, P) -> (
	te := newToricEmbedding(P, ring polynomialList_0);
	N := polynomialsToPolytope polynomialList;
	latticeN := latticePoints N;
	latticeNLen := length latticeN;
	teBigrMonNInTeRing := {sub((teBigradedMonomials( newToricEmbedding(polynomialList) ))_0,te#teRing)};
	scan((1..latticeNLen - 1), i -> teBigrMonNInTeRing = teBigrMonNInTeRing | {sub((teBigradedMonomials( newToricEmbedding(polynomialList) ))_i,te#teRing)});

	use te#teRing;
	torVar := teToricRing(te);
	-- substitute monomials by variables Ti interpreted in the ring te#teRing
	(M,C) := coefficients( sub( matrix {polynomialList}, te#teRing), Variables=>{te#teRing_0, te#teRing_1}, Monomials=>teBigrMonNInTeRing );
	Ts := transpose sub(relationNandP(polynomialsToPolytope(polynomialList), P),te#teRing);
	G := Ts*C; 
	G = matrix{{G_(0,0),G_(0,1),G_(0,2),G_(0,3)}};

	return sub (G,torVar);
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------
-- PURPOSE : Method that converts a matrix to a list starting from the up-left corner and ending on the down-right corner
matrixToList = method();

--   INPUT : 'm' a matrix
--  OUTPUT : 'l' a list which has the elements of 'm' listed from the up-left corner to the down-right corner
matrixToList(Matrix) := (m) -> (
	l := {};
	for i from 0 to (rank target m) - 1 do (
		row := {};
		for j from 0 to (rank source m) - 1 do (
			row = row | {m_j_i};
		);
		l = l | {row};
	);
	return l;
)

---------------------------------------------------------------------------
---------------------------------------------------------------------------
-- PURPOSE :Method that computes the right-most map of the Z-complex in degree nu
representationMatrix = method();

--   INPUT : (G, nu) where 'G' is a 1x4-matrix that contains the four polynomials that parametrisez the rational map, and 'nu' is an integer that represents the degree in the grading of the toric ring.
--  OUTPUT : 'RM' a matrix in 'X_0,X_1,X_2,X_3' whose rank drops when X_i = G_0_i
representationMatrix(Matrix,ZZ) := (G,nu) -> (
	Z1 := kernel koszul(1,G);
	toricRing := ring G;
	toricRingAndXs := toricRing[X_0..X_((rank source G) - 1)];

	-- create the nu strand of Z1 = kernel koszul(1,G)
	d :=  (degree G_0_0)_0; -- polynomials' degree
	Z1nu := super basis(nu+d, Z1); -- add d to compensate the natural shift in Z1
	Xnu := matrix{{X_0, X_1, X_2, X_3}} * substitute(Z1nu, toricRingAndXs); 

	-- matrix of the map Z1_nu -> A_nu[X]
	(m,MatInAXs) := coefficients(Xnu,Variables=>(matrixToList(sub(matrix{gens toricRing}, toricRingAndXs)))_0,Monomials=>substitute(basis(nu,toricRing),toricRingAndXs));

	QofXs := QQ[X_0..X_3];
	ListofTand0 := {X_0,X_1,X_2,X_3};
	for i from 0 to (length gens toricRing) - 1 do { ListofTand0=append(ListofTand0,0) };
	ppp := map(QofXs,toricRingAndXs,ListofTand0);

	M := ppp(MatInAXs);

	if (rank target M != rank M) then {
		print "WARNING: representation matrix is not of full rank (rank target M != rank M).";
	};

	return M;
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------
-- PURPOSE : Method that computes the gcd of the right-most map of the Z-complex in degree nu,
implicitEq = method();

--   INPUT : '(polynomialList,nu)', where polynomialList with the polynomials  {f0,f1,f2,f3} and 'nu' is an integer that represents the degree in the grading of the toric ring.
--  OUTPUT : A polynomial in 'X_0,X_1,X_2,X_3' that vanishes when is evaluated on X_i = fi
implicitEq(List, ZZ) := (polynomialList, nu) -> (
	N := MaxMinor(representationMatrix(teToricRationalMap (polynomialList),nu)); 
	-- A multiple of the implicit equation (taking several minors we erase extra factors)
	return det(N)
);

--   INPUT : '(polynomialList,nu,P)', where polynomialList with the polynomials  {f0,f1,f2,f3}, 'nu' is an integer that represents the degree in the grading of the toric ring and 'P' is a polyhedron that determines the toric ring.
--  OUTPUT : A polynomial in 'X_0,X_1,X_2,X_3' that vanishes when is evaluated on X_i = fi
implicitEq(List, ZZ, Polyhedron) := (polynomialList, nu, aPolyhedra) -> (
	N := MaxMinor(representationMatrix(teToricRationalMap (polynomialList, aPolyhedra),nu)); 
	-- A multiple of the implicit equation (taking several minors we erase extra factors)
	return det(N)
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------
-- PURPOSE : Method that verifies if the Euler Chrasteristric of the Z-complex is zero in the given degree
isGoodDegree = method();

--   INPUT : '(polynomialList,nu)', where polynomialList with the polynomials  {f0,f1,f2,f3} and 'nu' is an integer that represents the degree in the grading of the toric ring.
--  OUTPUT : 'True/False' a Boolean that says if the Euler Chrasteristric of the Z-complex is zero in the given degree 'nu'
isGoodDegree(Matrix, ZZ) := (polynomialList, nu) -> (
	F := polynomialList;
	Z0 := (ring F)^1;
	Z1 := kernel koszul(1,F);
	Z2 := kernel koszul(2,F);
	Z3 := kernel koszul(3,F);
	d := (degree (polynomialList_0))_0;

	eulerChar := hilbertFunction(nu,Z0)-hilbertFunction(nu+d,Z1)+hilbertFunction(nu+2*d,Z2)-hilbertFunction(nu+3*d,Z3);

	return (eulerChar == 0)
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------
-- PURPOSE :Method that computes the degree of the MacRae's invariant = degree of implicit equation = degree of det((Z)_nu)
degreeImplicitEq = method();

--   INPUT : '(polynomialList,nu)', where polynomialList is a matrix with the polynomials  {f0,f1,f2,f3} and 'nu' is an integer that represents the degree in the grading of the toric ring.
--  OUTPUT : 'd' an integer which is the degree of H^{deg(f)}.G, where H is the implicit equation, deg(f) is the degree of f={f0,f1,f2,f3}, and G is an extra factor that may appear.
degreeImplicitEq(Matrix, ZZ) := (polynomialList, nu) -> (
	F := polynomialList;
	Z0 := (ring F)^1;
	Z1 := kernel koszul(1,F);
	Z2 := kernel koszul(2,F);
	Z3 := kernel koszul(3,F);
	d := (degree (polynomialList_0_0))_0;

	return (hilbertFunction(nu+d,Z1)-2*hilbertFunction(nu+2*d,Z2)+3*hilbertFunction(nu+3*d,Z3))
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------
-- PURPOSE : Auxiliary method that takes from a given m x n - Matrix of rank r, a new full rank m x r - Matrix
MaxCol = method();

--   INPUT : 'M' an m x n - Matrix of rank r
--  OUTPUT : 'N' an m x r - Matrix of rank r
MaxCol(Matrix) := (M) -> {
     	--M is matrix in a Polynomial Ring
     	--MaxCol returns a submatrix of M with rank M col, of rank M,
     	--using generic substitution.
	rm:=rank(M); 
	i:=0;
	c:={};
	while #c<rm do ( 
		if rank M_(c|{i}) == #(c|{i}) then c=c|{i};
		i=i+1;
	);
	return(M_c)
     	--c is the list of columns to keep
};

---------------------------------------------------------------------------
---------------------------------------------------------------------------
-- PURPOSE : Method that takes from a given m x n - Matrix of rank r, a new r x r full rank Matrix"
MaxMinor = method();

--   INPUT : 'M' an m x n - Matrix of rank r
--  OUTPUT : 'N' an r x r - Matrix of rank r
MaxMinor (Matrix) := (M) -> {
	--Apply two times MaxCol to obtain a maximal minor
	--This is not very efficient...
	return (transpose MaxCol(transpose MaxCol(M)) )
};

----------------------------------------------------
--------------- DOCUMENTATION ----------------------
----------------------------------------------------

beginDocumentation()

---------------------------------------------------------
---------------------------------------------------------
-- Simple Doc information
---------------------------------------------------------
---------------------------------------------------------

--*******************************************************
-- DOCUMENTATION FOR PACKAGE
--*******************************************************

doc ///
	Key 
		Implicitization
	Headline
		A package for computing implicit equations of bigraded rational surfaces by means of approximation complexes
	Description
		Text
			@EM "Implicitization"@ is a package for computing implicit equations of toric rational surfaces 
			embedded in a projective space, by means of approximation complexes.

			We provide a method for computing a matrix representation (see @TO representationMatrix @) and the 
			implicit equation (see @TO implicitEq @ ) by means of the method developed in [BDD09] and [Bot10]. 
			As it is probably the most interesting case from a practical point of view, we restrict our 
			computations to parametrizations from toric surfaces embedded in some projective space. 
			This implementation allows to compute small examples for the better understanding of the theory 
			developed in [BDD09] and [Bot10].

			[BDD09] N. Botbol, A. Dickenstein, M. Dohm. "Matrix representations for toric parametrizations". e-prints server. Computer Aided Geometric Design. Vol 26, Issue 7 (2009), 757-771.
			[Bot10] N. Botbol. "Compactifications of rational maps and the implicit equations of their images". arXiv:0910.1340v2. To appear in Journal of Pure and Applied Algebra. (2010). 

///
--	isGoodDegree		OK
--	degreeImplicitEq	OK	
--	ToricEmbedding,		OK
--	newToricEmbedding 	OK
--	teToricRing 		OK
--	getPolytope,		OK
--	teToricRationalMap 	OK
--	representationMatrix	OK
--	implicitEq		OK
--	MaxMinor, 		OK
--	polynomialsToPolytope, 	OK
--	relationNandP 		OK


------------------------ \ documentation isGoodDegree / ------------------------

document {
     	Key => {isGoodDegree},
	Headline => "verifies if the Euler Chrasteristric of the Z-complex is zero in the given degree",
	Usage => " isGood = isGoodDegree(polynomialList,nu)",

	Inputs => {
		"polynomialList" => Matrix,
		"nu" => ZZ
	},
	Outputs => {
		"Bool" => Boolean
	},
	PARA {}, "The integer nu needs to be a good degree for the first parameter {f0,f1,f2,f3} that can be verified by doing isGoodDegree(polinomialList,nu).",
	PARA {}, TT "isGoodDegree", " verifies if the approximation complex Z associated to the polynomials given is acyclic in degree nu.",
	PARA {}, "Given a list of polynomials {f0,f1,f2,f3}, the approximation complex of cycles is bigraded. This funtion verifies if the strand of degree 'nu' is acyclic.",
	PARA {}, "Precisely, it computes the Euler characteristic of the nu-strand of the Z-complex, by computing the alternate sum of (-1)^i * hilbertFunction(nu+i*d,Z_i).",
	EXAMPLE {" S = QQ[s,t]; ",
		" f0 = s^2*t^2; ",
		" f1 = t^2+s^2*t; ",
		" f2 = t^3+s; ",
		" f3 = s^2+s^2*t^3; ",
		"isGoodDegree (teToricRationalMap {f0,f1,f2,f3},2)" 
		},
     
	SeeAlso => {"Polyhedra::Polyhedra", teToricRationalMap}
}

------------------------ \ documentation degreeImplicitEq / ------------------------

document {
	Key => {degreeImplicitEq},
	Headline => "computes the degree of det((Z)_nu)",
	Usage => " isGood = isGoodDegree(polynomialList,nu)",
	Inputs => {
		"polynomialList" => Matrix,
		"nu" => ZZ
	},
	Outputs => {
		"degH" => ZZ
	},
	PARA {}, " The integer 'nu' needs to be a 'good degree' for the first parameter '{f0,f1,f2,f3}' that can be verified by doing isGoodDegree(polinomialList,nu) ",
	PARA {}, TT "degreeImplicitEq", " computes the degree of the gcd of the maximal minors of the matrix representation of {f0,f1,f2,f3} in degree 'nu', equivalently, the degree of determinant of the Z-complex associated to {f0,f1,f2,f3} in degree 'nu'. This is:",
	PARA {}, " degreeImplicitEq(teToricRationalMap {f0,f1,f2,f3}, nu) = degree(representationMatrix({f0,f1,f2,f3},nu) that computes 'deg(det((Z(f0,f1,f2,f3))_nu))'",
	PARA {}, " Given a list of polynomials {f0,f1,f2,f3}, it computes the degree of the gcd of the maximal minors of the right-most map of the strand of degree 'nu'.",
	PARA {}, " Precisely, 'degreeImplicitEq' returns an integer which is the degree of H^{deg(f)}.G, where H is the implicit equation, deg(f) is the degree of f={f0,f1,f2,f3}, and G is an extra factor that may appear.",
	EXAMPLE {" S = QQ[s,t]; ",
		" f0 = s^2*t^2; ",
		" f1 = t^2+s^2*t; ",
		" f2 = t^3+s; ",
		" f3 = s^2+s^2*t^3; ",
		"degreeImplicitEq (teToricRationalMap {f0,f1,f2,f3},2)" },
	SeeAlso => {representationMatrix,teToricRationalMap, isGoodDegree, degree}
}

------------------------ \ documentation ToricEmbedding / ------------------------

document {
	Key => {ToricEmbedding},
	Headline => "Embedding of a variety with ambient space A^2 in an n-dimensional projective espace.",
	PARA {}, " Class that builds an embedding of the affine plane into an n-dimensional projective espace."
}


------------------------ \ documentation newToricEmbedding / ------------------------

document {
	Key => {(newToricEmbedding,List)},
	Headline => "constructs a ToricEmbedding",
	Usage => " TE = newToricEmbedding(polynomialList)",
	Inputs => {
		"polynomialList" => List => {"with polynomials {f0,f1,f2,f3}"},
	},
	
	Outputs => {
		"degH" => ToricEmbedding
	},
	PARA {}, " Given a list of four polynomials {f0,f1,f2,f3} in a polynomial ring in two variables, namely QQ[s,t] we can compute its Newton polytope, namely, the", TO convexHull, " of the union of the Newton polytope of all the polynomial in the list, with the method", TO polynomialsToPolytope, ". Take N= polynomialsToPolytope ({f0,f1,f2,f3}). Then", TO ToricEmbedding, " of {f0,f1,f2,f3} can be defined as", TO ToricEmbedding, " of N.",
	PARA {}, "From the polynedron N we construc the ambient ring for building the coordinate ring of the toric variety associated to N. Precisely, assume N has p+1 lattice points, this is length latticePoints N = p+1, and denote with T_0 ... T_p the lattice points of N. Consider the ring QQ[T_0..T_p]. Let J be the toric ideal of N, then newToricEmbedding provides the ring 'QQ[s,t,u,v,w,T_0..T_p,MonomialOrder=>Eliminate 5]' which has all the information of the coordinate ring QQ[T_0..T_p]/J.",
	PARA {}, " METHODS: teLattice, teLatticeLen,teRing",
	PARA {}, "teLattice: returns the list of lattice points of the polytope associated to the toric embedding. Precisely, when tEmb = newToricEmbedding(P,R), tEmb#teLattice = latticePoints P.",
	PARA {}, "teLatticeLen: returns the length list ''of lattice points of the polytope associated to the toric embedding. Precisely, when tEmb = newToricEmbedding(P,R), tEmb#teLattice = length latticePoints P.",
	PARA {}, "teRing: returns the ring 'QQ[s,t,u,v,w,T_0..T_tEmb#teLatticeLen,MonomialOrder=>Eliminate 5]'.",
	EXAMPLE {" S = QQ[s,t]; ",
		" f0 = s^2*t^2; ",
		" f1 = t^2+s^2*t; ",
		" f2 = t^3+s; ",
		" f3 = s^2+s^2*t^3; ",
		" l = {f0,f1,f2,f3};",
		" teFromPolyList = newToricEmbedding(l);"  },
	SeeAlso => {teToricRing, teToricRationalMap, polynomialsToPolytope}
}

document {
	Key => {newToricEmbedding,(newToricEmbedding,Polyhedron,Ring)},
	Headline => "constructs a ToricEmbedding",
	Usage => " TE = newToricEmbedding(P,R)",
	Inputs => {
		"P" => Polyhedron ,
		"R" => Ring
	},
	
	Outputs => {
		"degH" => ToricEmbedding
	},
	PARA {}, " Given a polynedron P and an affine ring R, we construc the ambient ring for building the coordinate ring of the toric variety associated to P. Precisely, assume P has p+1 lattice points, this is length latticePoints P = p+1, and denote with T_0 ... T_p the lattice points of P. Consider the ring QQ[T_0..T_p]. Let J be the toric ideal of P, then newToricEmbedding provides the ring 'QQ[s,t,u,v,w,T_0..T_p,MonomialOrder=>Eliminate 5]' which has all the information of the coordinate ring QQ[T_0..T_p]/J.",
	PARA {}, " METHODS: teLattice, teLatticeLen,teRing",
	PARA {}, "teLattice: returns the list of lattice points of the polytope associated to the toric embedding. Precisely, when tEmb = newToricEmbedding(P), tEmb#teLattice = latticePoints P.",
	PARA {}, "teLatticeLen: returns the length list ''of lattice points of the polytope associated to the toric embedding. Precisely, when tEmb = newToricEmbedding(P,R), tEmb#teLattice = length latticePoints P.",
	PARA {}, "teRing: returns the ring 'QQ[s,t,u,v,w,T_0..T_tEmb#teLatticeLen,MonomialOrder=>Eliminate 5]'.",
	EXAMPLE {" R = QQ[s,t]; ",
		"  P = convexHull matrix{{0,0,1,1},{0,1,0,2}};",
		" teFromPolyList = newToricEmbedding(P,R);"  },
	SeeAlso => {teToricRing, teToricRationalMap, polynomialsToPolytope}
}

------------------------ \ documentation teToricRing / ------------------------

document {
	Key => {teToricRing},
	Headline => "Returns the coordinate ring of the toric variety",
	Usage => " A = teToricRing(TE)",
	Inputs => {
		"TE" => ToricEmbedding
	},
	Outputs => {
		"degH" => QuotientRing
	},
	PARA {}, TT "teToricRing", " constructs from a",TO ToricEmbedding, " the associated coordinate ring.",   
	PARA {}, " Given a polyhedron P we construct the coordinate ring of the toric variety associated to P. Precisely, assume P has p+1 lattice points, this is", TO length, TO latticePoints, " P = p+1, and denote with T_0 ... T_p the lattice points of P. Consider the ring QQ[T_0..T_p]. Let J be the toric ideal of P, then teToricRing gives coordinate ring QQ[T_0..T_p]/J.",
	EXAMPLE {" S = QQ[s,t]; ",
		" f0 = s^2*t^2; ",
		" f1 = t^2+s^2*t; ",
		" f2 = t^3+s; ",
		" f3 = s^2+s^2*t^3; ",
		" l = {f0,f1,f2,f3};",
		" TE = newToricEmbedding(l);",
		" V = teToricRing(TE);"},
	SeeAlso => {newToricEmbedding, teToricRationalMap }
}

------------------------ \ documentation getPolytope / ------------------------

document {
	Key => {getPolytope},
	Headline => "returns the Polytope associated to the ToricEmbedding",
	Usage => " P = getPolytope(TE)",
	Inputs => {
		"TE" => ToricEmbedding
	},
	Outputs => {
		"P" => Polyhedron,
	},
	PARA {}, TT "getPolytope", "returns the polytope used for building the", TO ToricEmbedding ," TE.", 
	EXAMPLE {" S = QQ[s,t]; ",
		" f0 = s^2*t^2; ",
		" f1 = t^2+s^2*t; ",
		" f2 = t^3+s; ",
		" f3 = s^2+s^2*t^3; ",
		" l = {f0,f1,f2,f3}; ",
		" TE = newToricEmbedding(l); ",
		" getPolytope TE "
	},
}

------------------------ \ documentation teToricRationalMap / ------------------------
document {
	Key => "teToricRationalMap, (teToricRationalMap, Polyhedron, List)",
	Headline => "Computes the rational map G defined over the toric coordinate ring given by {f0,f1,f2,f3}",
	Usage => " G = teToricRationalMap(polynomialList,P)",
	Inputs => {
		"polynomialList" => List => {"with polynomials {f0,f1,f2,f3}"},
		"P" => Polyhedron
	},
	Outputs => {
		"G" => Matrix
	},
	PARA {}, " Given a polynedron P we can associate a Toric variety Tv_P with coordinate ring A. Precisely A can be computed as 'A = teToricRing (newToricEmbedding {f0,...,f3})'. Assume P has p+1 lattice points, this is length latticePoints P = p+1, and denote with T_0 ... T_p the lattice points of P. Consider the ring QQ[T_0..T_p]. Let J be the toric ideal of P, then teToricRing gives coordinate ring QQ[T_0..T_p]/J",
	PARA {}, "The polynomials f0,...,f3 induces a map homogeneous g=(g_0:...:g_3):Tv_P --> P^3",
	PARA {}, "teToricRationalMap ({f0,...,f3},P) gives the map g=(g_0:...:g_3):Tv_P --> P^3 induced by the polynomials f0,...,f3",
	EXAMPLE {" S = QQ[s,t]; ",
		" f0 = s^2*t^2; ",
		" f1 = t^2+s^2*t; ",
		" f2 = t^3+s; ",
		" f3 = s^2+s^2*t^3; ",
		"  l = {f0,f1,f2,f3}; ",
		"  P = convexHull matrix{{0,0,1,1},{0,1,0,2}};",
		"  G = teToricRationalMap(l,P); ",
	},
	SeeAlso => {newToricEmbedding, teToricRing, polynomialsToPolytope}
}

document {
	Key => "(teToricRationalMap, List)",
	Headline => "Computes the rational map G defined over the toric coordinate ring given by {f0,f1,f2,f3}",
	Usage => " G = teToricRationalMap(polynomialList)",
	Inputs => {
		"polynomialList" => List => {"with polynomials {f0,f1,f2,f3}"}
	},
	Outputs => {
		"G" => Matrix
	},
	PARA {}, "Given a list of polynomials f0,...,f3 induces a map homogeneous g=(g_0:...:g_3):Tv_N --> P^3, where Tv_N is the toric variety associates to the Newton polytope of the polynomials {f0,f1,f2,f3}",
	PARA {}, "teToricRationalMap ({f0,...,f3}) gives the map g=(g_0:...:g_3):Tv_N --> P^3 induced by the polynomials f0,...,f3, where N = polynomialsToPolytope({f0,...,f3}).",
	EXAMPLE {" S = QQ[s,t]; ",
		" f0 = s^2*t^2; ",
		" f1 = t^2+s^2*t; ",
		" f2 = t^3+s; ",
		" f3 = s^2+s^2*t^3; ",
		" l = {f0,f1,f2,f3}; ",
		" G = teToricRationalMap(l); ",
	},
	SeeAlso => {newToricEmbedding, teToricRing, polynomialsToPolytope}
}

------------------------ \ documentation representationMatrix / ------------------------

document {
	Key => "representationMatrix",
	Headline => "computes the right-most map of the Z-complex in degree nu)",
	Usage => " RM = representationMatrix(G,nu)",
	Inputs => {
		"G" => Matrix,
		"nu" => ZZ
	},
	Outputs => {
		"G" => Matrix
	},
	PARA {}, " NOTE: the first imput need to be a list or a matrix of homogenous polynomials. It is usually composed with 'teToricRationalMap'",
	PARA {}, " The integer 'nu' needs to be a 'good degree' for the first parameter '{f0,f1,f2,f3}' that can be verified by doing isGoodDegree(teToricRationalMap{f0,f1,f2,f3},nu) ",
	PARA {}, "representationMatrix(teToricRationalMap{f0,f1,f2,f3},nu) computes the right-most map of the Z-complex in degree nu. Its determinant vanishes on the implicit equation of the image of the map given by (f0:...:fn) in P^n. The output matrix in 'X_0,X_1,X_2,X_3' whose rank drops when X_i = fi",
	EXAMPLE {" S = QQ[s,t]; ",
		" f0 = s^2*t^2; ",
		" f1 = t^2+s^2*t; ",
		" f2 = t^3+s; ",
		" f3 = s^2+s^2*t^3; ",
		" l = {f0,f1,f2,f3}; ",
		"representationMatrix (teToricRationalMap(l),2)" },
	SeeAlso => { implicitEq, isGoodDegree}
}

------------------------ \ documentation implicitEq / ------------------------

document {
	Key => "implicitEq, (implicitEq, List, ZZ,Polyhedron)",
	Headline => "computes the gcd of the right-most map of the Z-complex in degree nu)",
	Usage => " iEq = implicitEq(polynomialList,nu,P)",
	Inputs => {
		"polynomialList" => List => {"with polynomials {f0,f1,f2,f3}"},
		"nu" => ZZ,
		"P" => Polyhedron
	},
	Outputs => {
		"G" => RingElement
	},
	PARA {}, " The integer 'nu' needs to be a 'good degree' for the first parameter '{f0,f1,f2,f3}' that can be verified by doing isGoodDegree(teToricRationalMap ({f0,f1,f2,f3},P),nu) ",
	PARA {},  "implicitEq({f0,f1,f2,f3},nu,P) computes the determinant of the maximal minors of representationMatrix({f0,f1,f2,f3},nu,P). Equivalently, it computes the determinant of the Z-complex in degree 'nu'. This is:",
	PARA {},  "implicitEq({f0,f1,f2,f3}, nu)=det((Z(f0,f1,f2,f3))_nu)",
	PARA {},  "We also have that: degree(implicitEq({f0,f1,f2,f3}, nu))=degreeImplicitEq(teToricRationalMap {f0,f1,f2,f3}, nu)",
	EXAMPLE {" S = QQ[s,t]; ",
		" f0 = s^2*t^2; ",
		" f1 = t^2+s^2*t; ",
		" f2 = t^3+s; ",
		" f3 = s^2+s^2*t^3; ",
		"  P = convexHull matrix{{0,0,1,1},{0,1,0,2}};",
		"implicitEq ({f0,f1,f2,f3},2,P)" },
	PARA {},  "For the programmer: it computes the determinant of a maximal minor of 'representationMatrix (teToricRationalMap{f0,f1,f2,f3},2,P)' using the function 'MaxMinor'",
	SeeAlso => { representationMatrix, degreeImplicitEq, MaxMinor, isGoodDegree}
}

document {
	Key => "(implicitEq, List, ZZ)",
	Headline => "computes the gcd of the right-most map of the Z-complex in degree nu)",
	Usage => " iEq = implicitEq(polynomialList,nu)",
	Inputs => {
		"polynomialList" => List => {"with polynomials {f0,f1,f2,f3}"},
		"nu" => ZZ
	},
	Outputs => {
		"G" => RingElement
	},
	PARA {}, " The integer 'nu' needs to be a 'good degree' for the first parameter '{f0,f1,f2,f3}' that can be verified by doing isGoodDegree(teToricRationalMap{f0,f1,f2,f3},nu) ",
	PARA {},  "implicitEq({f0,f1,f2,f3},nu) computes the determinant of the maximal minors of representationMatrix({f0,f1,f2,f3},nu). Equivalently, it computes the determinant of the Z-complex in degree 'nu'. This is:",
	PARA {},  "implicitEq({f0,f1,f2,f3}, nu)=det((Z(f0,f1,f2,f3))_nu)",
	PARA {},  "We also have that: degree(implicitEq({f0,f1,f2,f3}, nu))=degreeImplicitEq(teToricRationalMap {f0,f1,f2,f3}, nu)",
	EXAMPLE {" S = QQ[s,t]; ",
		" f0 = s^2*t^2; ",
		" f1 = t^2+s^2*t; ",
		" f2 = t^3+s; ",
		" f3 = s^2+s^2*t^3; ",
		"implicitEq ({f0,f1,f2,f3},2)" },
	PARA {},  "For the programmer: it computes the determinant of a maximal minor of 'representationMatrix (teToricRationalMap{f0,f1,f2,f3},2)' using the function 'MaxMinor'",
	SeeAlso => { representationMatrix, degreeImplicitEq, MaxMinor, isGoodDegree}
}
		
------------------------ \ documentation MaxMinor / ------------------------

document {
	Key => "MaxMinor",
	Headline => "Returns a maximal minor of the matrix of full rank.",
	Usage => " MM = MaxMinor(Mat)",
	Inputs => {
		"Mat" => Matrix
	},
	Outputs => {
		"MM" => Matrix
	},
	PARA {}, " From a given m x n - Matrix of rank r,", TO MaxMinor, " returns an r x r full rank Matrix",
	EXAMPLE { 
		" M = matrix {{1,2,3},{1,2,3},{4,5,6},{4,5,6}}",
		" MaxMinor M;"},
}
------------------------ \ documentation polynomialsToPolytope / ------------------------

document {
	Key => "polynomialsToPolytope",
	Headline => "Return the convexHull of the union of the Newton polytope of all the polynomial in the list",
	Usage => " N = polynomialsToPolytope(polynomialList)",
	Inputs => {
		"polynomialList" => List => {"with polynomials {f0,f1,f2,f3}"}
	},
	Outputs => {
		"N" => Polyhedron
	},
	PARA {}, " Given a list l={f_0..f_n} of polynomials in two variables this method computes the convexHull of the union of the Newton polytope of all the polynomial in the list. This is, if N_i= newtonPolytope(f_i), then polynomialsToPolytope(l)= convexHull(N_0,...,N_n).",
	EXAMPLE {" S = QQ[s,t]; ",
		" f0 = s^2*t^2; ",
		" f1 = t^2+s^2*t; ",
		" f2 = t^3+s; ",
		" f3 = s^2+s^2*t^3; ",
		" l = {f0,f1,f2,f3};",
		" N = polynomialsToPolytope(l);" },
	SeeAlso => {teToricRing, teToricRationalMap, representationMatrix}
}

------------------------ \ documentation relationNandP / ------------------------

document {
	Key => "relationNandP",
	Headline => "Returns the list of coordinates of N based on the smallest homothety of P containing N.",
	Usage => " N = polynomialsToPolytope(polynomialList)",
	Inputs => {
		"N" => Polyhedron,
		"P" => Polyhedron
	},
	Outputs => {
		"Mat" => Matrix
	},
	PARA {}, " Assume we are given two polyhedron, N and P, such that there exists a positive integer k such that kP contains N. Take the minimun k with this prpoerty. Asumme P has p+1 lattice points, this is length latticePoints P = p+1. Now denote with T_0 ... T_p the lattice points of P, and consider the ring QQ[T_0..T_p]. The method 'relationNandP(N,P)' gives a 1-column matrix with entries homogeneous elements of QQ[T_0..T_p] of degree k. The i-th row has the information of how the ih point in N can be written as a product of k T's, this is, as a sum of k points of P.",
	EXAMPLE {" S = QQ[s,t]; ",
		" f0 = s^2*t^2; ",
		" f1 = t^2+s^2*t; ",
		" f2 = t^3+s; ",
		" f3 = s^2+s^2*t^3; ",
		" l = {f0,f1,f2,f3};",
		" P = convexHull(matrix{{0,0,1},{0,1,0}});",
		" rNP = relationNandP(polynomialsToPolytope(l),P)" },
	PARA {}, "For the programmer: we use a simiular algorithm of that used in 'Pennies, Nickels, Dimes and Quarters (Computations in algebraic geometry with Macaulay 2 Editors: D. Eisenbud, D. Grayson,M. Stillman, and B. Sturmfels)'",
	SeeAlso => {teToricRing, teToricRationalMap, representationMatrix, polynomialsToPolytope}
}


---------------------------------------
-- TESTS
---------------------------------------

-- Test 0
-- Checking polynomialsToPolytope
TEST ///
S=QQ[s,t]; 
f0=s*t^6+2; 
f1=s*t^5-3*s*t^3; 
f2=s*t^4+5*s^2*t^6; 
f3=2+s^2*t^6; 
l = {f0,f1,f2,f3};

LATTICE = latticePoints polynomialsToPolytope(l);

assert(LATTICE_0 == matrix{{0},{0}})
assert(LATTICE_1 == matrix{{1},{3}})
assert(LATTICE_2 == matrix{{1},{4}})
assert(LATTICE_3 == matrix{{1},{5}})
assert(LATTICE_4 == matrix{{1},{6}})
assert(LATTICE_5 == matrix{{2},{6}})
///

-- Test 1
-- Checking lattice points of the Newton polytope
TEST ///
S=QQ[s,t]; 
f0=s*t^6+2; 
f1=s*t^5-3*s*t^3; 
f2=s*t^4+5*s^2*t^6; 
f3=2+s^2*t^6; 
l = {f0,f1,f2,f3};

TE = newToricEmbedding(l);

assert(length(latticePoints getPolytope TE) == 6)
///

-- Test 2
-- Checking that rings with different vars coincide
TEST ///
S=QQ[s,t]; 
f0=s*t^6+2; 
f1=s*t^5-3*s*t^3; 
f2=s*t^4+5*s^2*t^6; 
f3=2+s^2*t^6; 
l1 = {f0,f1,f2,f3};

TE1 = newToricEmbedding(l1);

R=QQ[a,b]; 
g0=a*b^6+2; 
g1=a*b^5-3*a*b^3; 
g2=a*b^4+5*a^2*b^6; 
g3=2+a^2*b^6; 
l2 = {g0,g1,g2,g3};

TE2 = newToricEmbedding(l2);

assert(latticePoints getPolytope TE1 == latticePoints getPolytope TE2)
///

-- Test 3
-- Test relationNandP
TEST ///
S=QQ[s,t]; 
f0=s*t^6+2; 
f1=s*t^5-3*s*t^3; 
f2=s*t^4+5*s^2*t^6; 
f3=2+s^2*t^6; 
l1 = {f0,f1,f2,f3};

N = polynomialsToPolytope(l1)

P = convexHull(matrix{{0,0,1},{0,2,2}});

matNandP = relationNandP(N,P);

expected = matrix{{T_0^3},{T_0*T_1*T_3},{T_0*T_2*T_3},{T_1*T_2*T_3},{T_2^2*T_3},{T_2*T_3^2}}

-- verifies by hand the construction of the monomials written on s,t and rewritten on T_0,T_1,T_2,T_3, by interpreting the points of N in 3*P.
assert(matNandP_0_0 == expected_0_0)
assert(matNandP_0_1 == expected_0_1)
assert(matNandP_0_2 == expected_0_2)
assert(matNandP_0_3 == expected_0_3)
assert(matNandP_0_4 == expected_0_4)
assert(matNandP_0_5 == expected_0_5)
///

-- Test 4
-- Test teToricRing
TEST ///
S=QQ[s,t]; 
f0=s*t^6+2; 
f1=s*t^5-3*s*t^3; 
f2=s*t^4+5*s^2*t^6; 
f3=2+s^2*t^6; 
l = {f0,f1,f2,f3};

TE = newToricEmbedding(l);

A = teToricRing(TE);

-- verifies by hand that the variables and the relations are the correct ones. Thus, if the quotient is correct.
gensSet = gens gb ideal A;
assert(gensSet_0_0 == T_3^2-T_2*T_4)
assert(gensSet_1_0 == T_2*T_3-T_1*T_4)
assert(gensSet_2_0 == T_2^2-T_1*T_3)
assert(gensSet_3_0 == T_1^2-T_0*T_5)

ambientVars = vars ambient A
assert(ambientVars_0_0 == T_0)
assert(ambientVars_1_0 == T_1)
assert(ambientVars_2_0 == T_2)
assert(ambientVars_3_0 == T_3)
assert(ambientVars_4_0 == T_4)
assert(ambientVars_5_0 == T_5)
///

-- Test 5
-- Test teToricRationalMap
TEST ///
S=QQ[s,t]; 
f0=s*t^6+2; 
f1=s*t^5-3*s*t^3; 
f2=s*t^4+5*s^2*t^6; 
f3=2+s^2*t^6; 
l = {f0,f1,f2,f3};

P := convexHull(matrix{{0,0,1},{0,2,2}});

myG := teToricRationalMap(l,P);

-- the number d is the relation between the sizes of the polytope P defining the toric variety and the Newton polytope of f0,f1,f2,f3.
d := 3 --This number is computed as: homothety (polynomialsToPolytope(l), P)
-- verifies that the degree of the rational map G defined over the toric coordinate ring is given by homogeneous polynomials of degree d.
assert(rank source myG == 4)
assert((degree myG_0_0)_0 == d)
assert((degree myG_1_0)_0 == d)
assert((degree myG_2_0)_0 == d)
assert((degree myG_3_0)_0 == d)

assert(isHomogeneous myG_0_0)
assert(isHomogeneous myG_1_0)
assert(isHomogeneous myG_2_0)
assert(isHomogeneous myG_3_0)
///

-- Test 6
-- Test teToricRationalMap
TEST ///
S=QQ[s,t]; 
f0=s*t^6+2; 
f1=s*t^5-3*s*t^3; 
f2=s*t^4+5*s^2*t^6; 
f3=2+s^2*t^6; 
l = {f0,f1,f2,f3};

myG := teToricRationalMap(l);

-- verifies that the Representation matrix RM drops its rank when is evaluated on f0,f1,f2,f3.
RM=representationMatrix (myG,2);

R=QQ[s,t,X_0..X_3];

FinRM=sub(sub(RM,R), {X_0=>(sub(f0,R)), X_1=>(sub(f1,R)), X_2=>(sub(f2,R)), X_3=>(sub(f3,R))});

assert(rank FinRM < rank RM)

-- verifies if the number of rows corresponds to the dimension of the toric coordinate ring in degree 2, which coincides with the number of lattice points in 2 times the newton polytope.
assert(rank target RM == length latticePoints(2*polynomialsToPolytope(l)))
///

-- Test 7
-- Test implicitEq
TEST ///
S=QQ[s,t]; 
f0=s*t^6+2; 
f1=s*t^5-3*s*t^3; 
f2=s*t^4+5*s^2*t^6; 
f3=2+s^2*t^6; 
l = {f0,f1,f2,f3};

EQ=implicitEq(l,2);

R=QQ[s,t,X_0..X_3];

FinEQ=sub(sub(EQ,R), {X_0=>(sub(f0,R)), X_1=>(sub(f1,R)), X_2=>(sub(f2,R)), X_3=>(sub(f3,R))});

assert(FinEQ == 0)
assert(isHomogeneous EQ)

strEQ = "24118074000*X_0^4*X_1^8*X_2^5+1064681172000*X_0^2*X_1^10*X_2^5-48236148000*X_0^5*X_1^6*X_2^6+573682176000*X_0^3*X_1^8*X_2^6+24118074000*X_0^6*X_1^4*X_2^7-434279880000*X_0^4*X_1^6*X_2^7+746775936000*X_0^2*X_1^8*X_2^7+1820232000*X_0^5*X_1^4*X_2^8-122007060000*X_0^3*X_1^6*X_2^8+26427708000*X_0^4*X_1^4*X_2^9+117044352000*X_0^2*X_1^6*X_2^9+995976000*X_0^3*X_1^4*X_2^10+7220826000*X_0^2*X_1^4*X_2^11+120590370000*X_0^5*X_1^6*X_2^5*X_3-1458323514000*X_0^3*X_1^8*X_2^5*X_3-120590370000*X_0^6*X_1^4*X_2^6*X_3+2243753622000*X_0^4*X_1^6*X_2^6*X_3-4020720768000*X_0^2*X_1^8*X_2^6*X_3-61887888000*X_0^5*X_1^4*X_2^7*X_3+1349332830000*X_0^3*X_1^6*X_2^7*X_3-267007428000*X_0^4*X_1^4*X_2^8*X_3-1109439990000*X_0^2*X_1^6*X_2^8*X_3-38877408000*X_0^3*X_1^4*X_2^9*X_3-108810378000*X_0^2*X_1^4*X_2^10*X_3-143343270000*X_0^4*X_1^6*X_2^5*X_3^2+1456185600000*X_0^2*X_1^8*X_2^5*X_3^2+263933640000*X_0^5*X_1^4*X_2^6*X_3^2-3721035024000*X_0^3*X_1^6*X_2^6*X_3^2+707778324000*X_0^4*X_1^4*X_2^7*X_3^2+2316030570000*X_0^2*X_1^6*X_2^7*X_3^2+315071856000*X_0^3*X_1^4*X_2^8*X_3^2+549040356000*X_0^2*X_1^4*X_2^9*X_3^2+23826150000*X_0^3*X_1^6*X_2^5*X_3^3-167169420000*X_0^4*X_1^4*X_2^6*X_3^3+1525517550000*X_0^2*X_1^6*X_2^6*X_3^3-732900960000*X_0^3*X_1^4*X_2^7*X_3^3-940081140000*X_0^2*X_1^4*X_2^8*X_3^3-1073250000*X_0^2*X_1^6*X_2^5*X_3^4+24899400000*X_0^3*X_1^4*X_2^6*X_3^4+62892450000*X_0^2*X_1^4*X_2^7*X_3^4-1073250000*X_0^2*X_1^4*X_2^6*X_3^5"

assert(toString EQ == strEQ)
///

-- Test 8
-- Test isGoodDegree
TEST ///
S=QQ[s,t]; 
f0=s*t^6+2; 
f1=s*t^5-3*s*t^3; 
f2=s*t^4+5*s^2*t^6; 
f3=2+s^2*t^6; 
l = {f0,f1,f2,f3};

myG := teToricRationalMap(l);

assert(not isGoodDegree(myG,0))
assert(not isGoodDegree(myG,1))
assert(isGoodDegree(myG,2))
assert(isGoodDegree(myG,3))
///

-- Test 9
-- Test degreeImplicitEq
TEST ///
S=QQ[s,t]; 
f0=s*t^6+2; 
f1=s*t^5-3*s*t^3; 
f2=s*t^4+5*s^2*t^6; 
f3=2+s^2*t^6; 
l = {f0,f1,f2,f3};

myG := teToricRationalMap(l);

assert(degreeImplicitEq(myG, 0) == 0)
assert(degreeImplicitEq(myG, 1) == 2)
assert(degreeImplicitEq(myG, 2) == 6)
assert(degreeImplicitEq(myG, 5) == 6)
assert(degreeImplicitEq(myG, 10) == 6)
///

---------------- TESTS with EXAMPLE 8-bis

-- Test 10
-- Checking polynomialsToPolytope
TEST ///
 S = QQ[s,t];
f0 = 3*s^2*t^2-2*s*t^3-s^2*t+s*t^2-3*s*t-t^2+4*t-1;
f1 = 3*s^2*t^2-s^2*t-3*s*t^2-s*t+t^2+1+t^3+s^2;
f2 = 2*s^2*t^3-3*s^2*t^2-s^2*t+s*t^2+3*s*t-3*t^2+2*t-1;
f3 = 2*s^2*t^3-3*s^2*t^2-2*s*t^3+s^2*t+5*s*t^2-3*s*t-3*t^2+4*t-1;
l = {f0,f1,f2,f3};

LATTICE = latticePoints polynomialsToPolytope(l);

assert(LATTICE_0 == matrix{{0},{0}})
assert(LATTICE_1 == matrix{{0},{1}})
assert(LATTICE_2 == matrix{{0},{2}})
assert(LATTICE_3 == matrix{{0},{3}})
assert(LATTICE_4 == matrix{{1},{0}})
assert(LATTICE_5 == matrix{{1},{1}})
assert(LATTICE_6 == matrix{{1},{2}})
assert(LATTICE_7 == matrix{{1},{3}})
assert(LATTICE_8 == matrix{{2},{0}})
assert(LATTICE_9 == matrix{{2},{1}})
assert(LATTICE_10 == matrix{{2},{2}})
assert(LATTICE_11 == matrix{{2},{3}})
///

-- Test 11
-- Checking lattice points of the Newton polytope
TEST ///
 S = QQ[s,t];
f0 = 3*s^2*t^2-2*s*t^3-s^2*t+s*t^2-3*s*t-t^2+4*t-1;
f1 = 3*s^2*t^2-s^2*t-3*s*t^2-s*t+t^2+1+t^3+s^2;
f2 = 2*s^2*t^3-3*s^2*t^2-s^2*t+s*t^2+3*s*t-3*t^2+2*t-1;
f3 = 2*s^2*t^3-3*s^2*t^2-2*s*t^3+s^2*t+5*s*t^2-3*s*t-3*t^2+4*t-1;
l = {f0,f1,f2,f3};

TE = newToricEmbedding(l);

assert(length(latticePoints getPolytope TE) == 12)
///

-- Test 12
-- Checking that rings with different vars coincide
TEST ///
 S = QQ[s,t];
f0 = 3*s^2*t^2-2*s*t^3-s^2*t+s*t^2-3*s*t-t^2+4*t-1;
f1 = 3*s^2*t^2-s^2*t-3*s*t^2-s*t+t^2+1+t^3+s^2;
f2 = 2*s^2*t^3-3*s^2*t^2-s^2*t+s*t^2+3*s*t-3*t^2+2*t-1;
f3 = 2*s^2*t^3-3*s^2*t^2-2*s*t^3+s^2*t+5*s*t^2-3*s*t-3*t^2+4*t-1;
l1 = {f0,f1,f2,f3};

TE1 = newToricEmbedding(l1);

 S = QQ[a,b];
g0 = 3*a^2*b^2-2*a*b^3-a^2*b+a*b^2-3*a*b-b^2+4*b-1;
g1 = 3*a^2*b^2-a^2*b-3*a*b^2-a*b+b^2+1+b^3+a^2;
g2 = 2*a^2*b^3-3*a^2*b^2-a^2*b+a*b^2+3*a*b-3*b^2+2*b-1;
g3 = 2*a^2*b^3-3*a^2*b^2-2*a*b^3+a^2*b+5*a*b^2-3*a*b-3*b^2+4*b-1;
l2 = {g0,g1,g2,g3};

TE2 = newToricEmbedding(l2);

assert(latticePoints getPolytope TE1 == latticePoints getPolytope TE2)
///

-- Test 13
-- Test relationNandP
TEST ///
 S = QQ[s,t];
f0 = 3*s^2*t^2-2*s*t^3-s^2*t+s*t^2-3*s*t-t^2+4*t-1;
f1 = 3*s^2*t^2-s^2*t-3*s*t^2-s*t+t^2+1+t^3+s^2;
f2 = 2*s^2*t^3-3*s^2*t^2-s^2*t+s*t^2+3*s*t-3*t^2+2*t-1;
f3 = 2*s^2*t^3-3*s^2*t^2-2*s*t^3+s^2*t+5*s*t^2-3*s*t-3*t^2+4*t-1;
l1 = {f0,f1,f2,f3};

N = polynomialsToPolytope(l1)

P = convexHull(matrix{{0,0,1,1},{0,1,0,1}});

matNandP = relationNandP(N,P);

expected = matrix{{T_0^3},{T_0^2*T_1},{T_0*T_1^2},{T_1^3},{T_0^2*T_2},{T_0^2*T_3},{T_0*T_1*T_3},{T_1^2*T_3},{T_0*T_2^2},{T_0*T_2*T_3},{T_0*T_3^2},{T_1*T_3^2}}

-- verifies by hand the construction of the monomials written on s,t and rewritten on T_0,T_1,T_2,T_3, by interpreting the points of N in 3*P.
assert(matNandP_0 == expected_0)
///

-- Test 14
-- Test teToricRing
TEST ///
 S = QQ[s,t];
f0 = 3*s^2*t^2-2*s*t^3-s^2*t+s*t^2-3*s*t-t^2+4*t-1;
f1 = 3*s^2*t^2-s^2*t-3*s*t^2-s*t+t^2+1+t^3+s^2;
f2 = 2*s^2*t^3-3*s^2*t^2-s^2*t+s*t^2+3*s*t-3*t^2+2*t-1;
f3 = 2*s^2*t^3-3*s^2*t^2-2*s*t^3+s^2*t+5*s*t^2-3*s*t-3*t^2+4*t-1;
l = {f0,f1,f2,f3};

TE = newToricEmbedding(l);

A = teToricRing(TE);

-- verifies by hand that the variables and the relations are the correct ones. Thus, if the quotient is correct.
gensSet = gens gb ideal A;
    
assert(gensSet_0_0 == T_10^2-T_9*T_11)
assert(gensSet_1_0 == T_9*T_10-T_8*T_11)
assert(gensSet_2_0 == T_7*T_10-T_6*T_11)
assert(gensSet_3_0 == T_6*T_10-T_5*T_11)
assert(gensSet_4_0 == T_5*T_10-T_4*T_11)
assert(gensSet_5_0 == T_3*T_10-T_2*T_11)
assert(gensSet_6_0 == T_2*T_10-T_1*T_11)
assert(gensSet_7_0 == T_1*T_10-T_0*T_11)
assert(gensSet_8_0 == T_9^2-T_8*T_10)
assert(gensSet_9_0 == T_7*T_9-T_5*T_11)
assert(gensSet_10_0 == T_6*T_9-T_4*T_11)
assert(gensSet_11_0 == T_5*T_9-T_4*T_10)
assert(gensSet_12_0 == T_3*T_9-T_1*T_11)
assert(gensSet_13_0 == T_2*T_9-T_0*T_11)
assert(gensSet_14_0 == T_1*T_9-T_0*T_10)
assert(gensSet_15_0 == T_7*T_8-T_4*T_11)
assert(gensSet_16_0 == T_6*T_8-T_4*T_10)
assert(gensSet_17_0 == T_5*T_8-T_4*T_9)
assert(gensSet_18_0 == T_3*T_8-T_0*T_11)
assert(gensSet_19_0 == T_2*T_8-T_0*T_10)
assert(gensSet_20_0 == T_1*T_8-T_0*T_9)
assert(gensSet_21_0 == T_7^2-T_3*T_11)
assert(gensSet_22_0 == T_6*T_7-T_2*T_11)
assert(gensSet_23_0 == T_5*T_7-T_1*T_11)
assert(gensSet_24_0 == T_4*T_7-T_0*T_11)
assert(gensSet_25_0 == T_6^2-T_1*T_11)
assert(gensSet_26_0 == T_5*T_6-T_0*T_11)
assert(gensSet_27_0 == T_4*T_6-T_0*T_10)
assert(gensSet_28_0 == T_3*T_6-T_2*T_7)
assert(gensSet_29_0 == T_2*T_6-T_1*T_7)
assert(gensSet_30_0 == T_1*T_6-T_0*T_7)
assert(gensSet_31_0 == T_5^2-T_0*T_10)
assert(gensSet_32_0 == T_4*T_5-T_0*T_9)
assert(gensSet_33_0 == T_3*T_5-T_1*T_7)
assert(gensSet_34_0 == T_2*T_5-T_0*T_7)
assert(gensSet_35_0 == T_1*T_5-T_0*T_6)
assert(gensSet_36_0 == T_4^2-T_0*T_8)
assert(gensSet_37_0 == T_3*T_4-T_0*T_7)
assert(gensSet_38_0 == T_2*T_4-T_0*T_6)
assert(gensSet_39_0 == T_1*T_4-T_0*T_5)
assert(gensSet_40_0 == T_2^2-T_1*T_3)
assert(gensSet_41_0 == T_1*T_2-T_0*T_3)
assert(gensSet_42_0 == T_1^2-T_0*T_2)

ambientVars = vars ambient A
assert(ambientVars_0_0 == T_0)
assert(ambientVars_1_0 == T_1)
assert(ambientVars_2_0 == T_2)
assert(ambientVars_3_0 == T_3)
assert(ambientVars_4_0 == T_4)
assert(ambientVars_5_0 == T_5)
assert(ambientVars_6_0 == T_6)
assert(ambientVars_7_0 == T_7)
assert(ambientVars_8_0 == T_8)
assert(ambientVars_9_0 == T_9)
assert(ambientVars_10_0 == T_10)
assert(ambientVars_11_0 == T_11)
///

-- Test 15
-- Test teToricRationalMap
TEST ///
 S = QQ[s,t];
f0 = 3*s^2*t^2-2*s*t^3-s^2*t+s*t^2-3*s*t-t^2+4*t-1;
f1 = 3*s^2*t^2-s^2*t-3*s*t^2-s*t+t^2+1+t^3+s^2;
f2 = 2*s^2*t^3-3*s^2*t^2-s^2*t+s*t^2+3*s*t-3*t^2+2*t-1;
f3 = 2*s^2*t^3-3*s^2*t^2-2*s*t^3+s^2*t+5*s*t^2-3*s*t-3*t^2+4*t-1;
l = {f0,f1,f2,f3};

P = convexHull(matrix{{0,0,1,1},{0,1,0,1}});

myG := teToricRationalMap(l,P);

-- the number d is the relation between the sizes of the polytope P defining the toric variety and the Newton polytope of f0,f1,f2,f3.
d := 3 --This number is computed as: homothety (polynomialsToPolytope(l), P)
-- verifies that the degree of the rational map G defined over the toric coordinate ring is given by homogeneous polynomials of degree d.
assert(rank source myG == 4)
assert((degree myG_0_0)_0 == d)
assert((degree myG_1_0)_0 == d)
assert((degree myG_2_0)_0 == d)
assert((degree myG_3_0)_0 == d)

assert(isHomogeneous myG_0_0)
assert(isHomogeneous myG_1_0)
assert(isHomogeneous myG_2_0)
assert(isHomogeneous myG_3_0)
///

-- ACA
-- Test 16
-- Test teToricRationalMap
TEST ///
 S = QQ[s,t];
f0 = 3*s^2*t^2-2*s*t^3-s^2*t+s*t^2-3*s*t-t^2+4*t-1;
f1 = 3*s^2*t^2-s^2*t-3*s*t^2-s*t+t^2+1+t^3+s^2;
f2 = 2*s^2*t^3-3*s^2*t^2-s^2*t+s*t^2+3*s*t-3*t^2+2*t-1;
f3 = 2*s^2*t^3-3*s^2*t^2-2*s*t^3+s^2*t+5*s*t^2-3*s*t-3*t^2+4*t-1;
l = {f0,f1,f2,f3};

myG := teToricRationalMap(l);

-- verifies that the Representation matrix RM drops its rank when is evaluated on f0,f1,f2,f3.
RM=representationMatrix (myG,2);

R=QQ[s,t,X_0..X_3];

FinRM=sub(sub(RM,R), {X_0=>(sub(f0,R)), X_1=>(sub(f1,R)), X_2=>(sub(f2,R)), X_3=>(sub(f3,R))});

assert(rank FinRM < rank RM)

-- verifies if the number of rows corresponds to the dimension of the toric coordinate ring in degree 2, which coincides with the number of lattice points in 2 times the newton polytope.
assert(rank target RM == length latticePoints(2*polynomialsToPolytope(l)))
///

---------------- TESTS with EXAMPLE 3

-- Test 17
-- Test implicitEq
TEST ///

 S = QQ[s,t];
f0 = t^2+3*s*t^2
f1 = t^2+2*t
f2 = t^3-s*t
f3 = s+s*t^3
l = {f0,f1,f2,f3};

EQ=implicitEq(l,2);

R=QQ[s,t,X_0..X_3];

FinEQ=sub(sub(EQ,R), {X_0=>(sub(f0,R)), X_1=>(sub(f1,R)), X_2=>(sub(f2,R)), X_3=>(sub(f3,R))});

assert(FinEQ == 0)
assert(isHomogeneous EQ)

strEQ = "36283613184*X_0^17*X_1+93300719616*X_0^16*X_1^2-233251799040*X_0^15*X_1^3-77750599680*X_0^14*X_1^4+155501199360*X_0^13*X_1^5+171051319296*X_0^12*X_1^6-145134452736*X_0^11*X_1^7+54425419776*X_0^16*X_1*X_2+318777458688*X_0^15*X_1^2*X_2-1422835974144*X_0^14*X_1^3*X_2+1632762593280*X_0^13*X_1^4*X_2-583129497600*X_0^12*X_1^5*X_2+256576978944*X_0^11*X_1^6*X_2-474278658048*X_0^10*X_1^7*X_2+217701679104*X_0^9*X_1^8*X_2+27212709888*X_0^15*X_1*X_2^2+1764938612736*X_0^14*X_1^2*X_2^2-6235598094336*X_0^13*X_1^3*X_2^2+8058849656832*X_0^12*X_1^4*X_2^2-4839974830080*X_0^11*X_1^5*X_2^2+1586112233472*X_0^10*X_1^6*X_2^2-598679617536*X_0^9*X_1^7*X_2^2+345990168576*X_0^8*X_1^8*X_2^2-108850839552*X_0^7*X_1^9*X_2^2-975122104320*X_0^14*X_1*X_2^3+5211233943552*X_0^13*X_1^2*X_2^3-12436208418816*X_0^12*X_1^3*X_2^3+16307540361216*X_0^11*X_1^4*X_2^3-11946379640832*X_0^10*X_1^5*X_2^3+4686417395712*X_0^9*X_1^6*X_2^3-1063887372288*X_0^8*X_1^7*X_2^3+274070863872*X_0^7*X_1^8*X_2^3-75806834688*X_0^6*X_1^9*X_2^3+18141806592*X_0^5*X_1^10*X_2^3-1469486333952*X_0^13*X_1*X_2^4+6956734906368*X_0^12*X_1^2*X_2^4-14269178806272*X_0^11*X_1^3*X_2^4+16986562265088*X_0^10*X_1^4*X_2^4-12986293911552*X_0^9*X_1^5*X_2^4+6268642099200*X_0^8*X_1^6*X_2^4-1661919068160*X_0^7*X_1^7*X_2^4+204095324160*X_0^6*X_1^8*X_2^4-29156474880*X_0^5*X_1^9*X_2^4-734743166976*X_0^12*X_1*X_2^5+3848654684160*X_0^11*X_1^2*X_2^5-8397064765440*X_0^10*X_1^3*X_2^5+9971514408960*X_0^9*X_1^4*X_2^5-7172492820480*X_0^8*X_1^5*X_2^5+3358825906176*X_0^7*X_1^6*X_2^5-1049633095680*X_0^6*X_1^7*X_2^5+174938849280*X_0^5*X_1^8*X_2^5-122457194496*X_0^11*X_1*X_2^6+734743166976*X_0^10*X_1^2*X_2^6-1836857917440*X_0^9*X_1^3*X_2^6+2449143889920*X_0^8*X_1^4*X_2^6-1836857917440*X_0^7*X_1^5*X_2^6+734743166976*X_0^6*X_1^6*X_2^6-122457194496*X_0^5*X_1^7*X_2^6+67383853056*X_0^16*X_1*X_3-559804317696*X_0^15*X_1^2*X_3+544254197760*X_0^14*X_1^3*X_3+907090329600*X_0^13*X_1^4*X_3-1399510794240*X_0^12*X_1^5*X_3+295452278784*X_0^11*X_1^6*X_3+145134452736*X_0^10*X_1^7*X_3+393936371712*X_0^15*X_1*X_2*X_3-2368801603584*X_0^14*X_1^2*X_2*X_3+2705720868864*X_0^13*X_1^3*X_2*X_3+2975256281088*X_0^12*X_1^4*X_2*X_3-7235989143552*X_0^11*X_1^5*X_2*X_3+4058581303296*X_0^10*X_1^6*X_2*X_3-383569625088*X_0^9*X_1^7*X_2*X_3-145134452736*X_0^8*X_1^8*X_2*X_3+506674741248*X_0^14*X_1*X_2^2*X_3-3291442053120*X_0^13*X_1^2*X_2^2*X_3+3542835658752*X_0^12*X_1^3*X_2^2*X_3+6190243577856*X_0^11*X_1^4*X_2^2*X_3-14856843755520*X_0^10*X_1^5*X_2^2*X_3+10695890829312*X_0^9*X_1^6*X_2^2*X_3-3058190254080*X_0^8*X_1^7*X_2^2*X_3+234547642368*X_0^7*X_1^8*X_2^2*X_3+36283613184*X_0^6*X_1^9*X_2^2*X_3-1574449643520*X_0^13*X_1*X_2^3*X_3+6935353491456*X_0^12*X_1^2*X_2^3*X_3-13155401465856*X_0^11*X_1^3*X_2^3*X_3+16700828811264*X_0^10*X_1^4*X_2^3*X_3-17439459508224*X_0^9*X_1^5*X_2^3*X_3+12758873407488*X_0^8*X_1^6*X_2^3*X_3-4874962599936*X_0^7*X_1^7*X_2^3*X_3+707530457088*X_0^6*X_1^8*X_2^3*X_3-58312949760*X_0^5*X_1^9*X_2^3*X_3-1778544967680*X_0^12*X_1*X_2^4*X_3+10222260092928*X_0^11*X_1^2*X_2^4*X_3-23564262998016*X_0^10*X_1^3*X_2^4*X_3+28205973798912*X_0^9*X_1^4*X_2^4*X_3-19377393205248*X_0^8*X_1^5*X_2^4*X_3+8554509729792*X_0^7*X_1^6*X_2^4*X_3-2863165833216*X_0^6*X_1^7*X_2^4*X_3+600623382528*X_0^5*X_1^8*X_2^4*X_3-454841008128*X_0^11*X_1*X_2^5*X_3+3096417632256*X_0^10*X_1^2*X_2^5*X_3-8659473039360*X_0^9*X_1^3*X_2^5*X_3+12770535997440*X_0^8*X_1^4*X_2^5*X_3-10496330956800*X_0^7*X_1^5*X_2^5*X_3+4565903966208*X_0^6*X_1^6*X_2^5*X_3-822212591616*X_0^5*X_1^7*X_2^5*X_3+233251799040*X_0^14*X_1^2*X_3^2-233251799040*X_0^13*X_1^3*X_3^2-699755397120*X_0^12*X_1^4*X_3^2+1166258995200*X_0^11*X_1^5*X_3^2-466503598080*X_0^10*X_1^6*X_3^2+606454677504*X_0^14*X_1*X_2*X_3^2-2682395688960*X_0^13*X_1^2*X_2*X_3^2+3732028784640*X_0^12*X_1^3*X_2*X_3^2-1632762593280*X_0^11*X_1^4*X_2*X_3^2+233251799040*X_0^10*X_1^5*X_2*X_3^2-723080577024*X_0^9*X_1^6*X_2*X_3^2+466503598080*X_0^8*X_1^7*X_2*X_3^2+606454677504*X_0^13*X_1*X_2^2*X_3^2+1994302881792*X_0^12*X_1^2*X_2^2*X_3^2-9435035271168*X_0^11*X_1^3*X_2^2*X_3^2+11009484914688*X_0^10*X_1^4*X_2^2*X_3^2-5341466198016*X_0^9*X_1^5*X_2^2*X_3^2+1667750363136*X_0^8*X_1^6*X_2^2*X_3^2-384865468416*X_0^7*X_1^7*X_2^2*X_3^2-116625899520*X_0^6*X_1^8*X_2^2*X_3^2+151613669376*X_0^12*X_1*X_2^3*X_3^2+4425952886784*X_0^11*X_1^2*X_2^3*X_3^2-19173297881088*X_0^10*X_1^3*X_2^3*X_3^2+29628809773056*X_0^9*X_1^4*X_2^3*X_3^2-20607796445184*X_0^8*X_1^5*X_2^3*X_3^2+6210329149440*X_0^7*X_1^6*X_2^3*X_3^2-863031656448*X_0^6*X_1^7*X_2^3*X_3^2+227420504064*X_0^5*X_1^8*X_2^3*X_3^2+1364523024384*X_0^10*X_1^2*X_2^4*X_3^2-6822615121920*X_0^9*X_1^3*X_2^4*X_3^2+13645230243840*X_0^8*X_1^4*X_2^4*X_3^2-13645230243840*X_0^7*X_1^5*X_2^4*X_3^2+6822615121920*X_0^6*X_1^6*X_2^4*X_3^2-1364523024384*X_0^5*X_1^7*X_2^4*X_3^2-1819364032512*X_0^13*X_1^2*X_3^3+5458092097536*X_0^12*X_1^3*X_3^3-5458092097536*X_0^11*X_1^4*X_3^3+1819364032512*X_0^10*X_1^5*X_3^3-1819364032512*X_0^12*X_1^2*X_2*X_3^3+7277456130048*X_0^11*X_1^3*X_2*X_3^3-10916184195072*X_0^10*X_1^4*X_2*X_3^3+7277456130048*X_0^9*X_1^5*X_2*X_3^3-1819364032512*X_0^8*X_1^6*X_2*X_3^3-454841008128*X_0^11*X_1^2*X_2^2*X_3^3+2274205040640*X_0^10*X_1^3*X_2^2*X_3^3-4548410081280*X_0^9*X_1^4*X_2^2*X_3^3+4548410081280*X_0^8*X_1^5*X_2^2*X_3^3-2274205040640*X_0^7*X_1^6*X_2^2*X_3^3+454841008128*X_0^6*X_1^7*X_2^2*X_3^3"

assert(toString EQ == strEQ)
///


-- Test 18
-- Test isGoodDegree
TEST ///
 S = QQ[s,t];
f0 = t^2+3*s*t^2
f1 = t^2+2*t
f2 = t^3-s*t
f3 = s+s*t^3
l = {f0,f1,f2,f3};

myG := teToricRationalMap(l);

assert(not isGoodDegree(myG,0))
assert(isGoodDegree(myG,1))
assert(isGoodDegree(myG,2))
assert(isGoodDegree(myG,3))
///

-- Test 19
-- Test degreeImplicitEq
TEST ///
 S = QQ[s,t];
f0 = t^2+3*s*t^2
f1 = t^2+2*t
f2 = t^3-s*t
f3 = s+s*t^3
l = {f0,f1,f2,f3};

myG := teToricRationalMap(l);

assert(degreeImplicitEq(myG, 0) == 0)
assert(degreeImplicitEq(myG, 1) == 5)
assert(degreeImplicitEq(myG, 2) == 5)
assert(degreeImplicitEq(myG, 5) == 5)
assert(degreeImplicitEq(myG, 10) == 5)
///

end


--*******************************************************
-- EXAMPLES FOR PACKAGE
--*******************************************************
-----------------------------------------------
-- Example 0th example (2,2) 
-----------------------------------------------

	 S = QQ[s,t];
	f0 = s*t+t;
	f1 = s;
	f2 = t^2*s+t^2;
	f3 = s^2*t+s*t;


-----------------------------------------------
-- Example first example (2,3) 
-----------------------------------------------
	 S = QQ[s,t];
	f0 = s^2*t^2
	f1 = t^2+s^2*t
	f2 = t^3+s
	f3 = s^2+s^2*t^3


-----------------------------------------------
-- Example second example (1,3) 
-----------------------------------------------
	 S = QQ[s,t];
	f0 = s*t^2
	f1 = t^2+s*t
	f2 = t^3+s
	f3 = s+s*t^3


-----------------------------------------------
-- Example third example (1,3) 
-----------------------------------------------
	 S = QQ[s,t];
	f0 = t^2+3*s*t^2
	f1 = t^2+2*t
	f2 = t^3-s*t
	f3 = s+s*t^3


-----------------------------------------------
-- Example fourth example (4,6)   --- doesn't work, gcd is not 1
-----------------------------------------------
	 S = QQ[s,t];
	f0 = s^4*t^2
	f1 = s^4*t^6-s^3*t^2
	f2 = t^3+2*s^2*t^4
	f3 = s^5*t+s^3*t^5


-----------------------------------------------
-- Example fifth example (4,5)  without basepoints
-----------------------------------------------
	 S = QQ[s,t];
	f0 = s^4^5-2+t^5
	f1 = s^4*t^5+1
	f2 = t^3+2*s^2*t^3
	f3 = s*t+s^3*t^5


-----------------------------------------------
-- Example sixth example (2,4)  
-----------------------------------------------
	 S = QQ[s,t];
	f0 = s^2+t^4
	f1 = s^2*t^4+1
	f2 = t^2+2*s*t^3^1
	f3 = s*t+s^2*t^4


-----------------------------------------------
-- Example seventh example (4,6)  
-----------------------------------------------
	 S = QQ[s,t];
	f0 = s^4+t^6
	f1 = s^4*t^6+1
	f2 = t^3+2*s^3*t^5^1
	f3 = s^5*t+s^4*t^6



-----------------------------------------------
-- Example eighth example (3,6)  
-----------------------------------------------
	 S = QQ[s,t];
	f0 = s^3^6+t^6
	f1 = s^3*t^6+1
	f2 = t^3+2*s^2*t^5
	f3 = s^5*t+s^3*t^6

-----------------------------------------------
-- Example eighth (bis) example (3,6)  
-----------------------------------------------
	 S = QQ[s,t];
	f0 = s^2+s^3*t;
	f1 = s^3*t^6+1;
	f2 = s*t^2+2*s^3*t^5;
	f3 = s^2+s^3*t^6;

	-- Proposed examples
	load "Implicitization.m2";
	L={f0,f1,f2,f3};
	-----
	latticePoints polynomialsToPolytope L
	tEmb=newToricEmbedding L;
	-----
	g = teToricRationalMap L;
	rM=representationMatrix (teToricRationalMap L,2);
	R=QQ[s,t,X_0..X_3];
	Eq=sub(sub(rM,R), {X_0=>(sub(f0,R)), X_1=>(sub(f1,R)), X_2=>(sub(f2,R)), X_3=>(sub(f3,R))});
	rank rM
	rank Eq
	-----
	Q1=convexHull matrix{{0,0,1,1},{0,1,0,2}}
	gQ1 = teToricRationalMap (L,Q1)
	rMQ1=representationMatrix (gQ1,5);
	R=QQ[s,t,X_0..X_3];
	EqQ1=sub(sub(rMQ1,R), {X_0=>(sub(f0,R)), X_1=>(sub(f1,R)), X_2=>(sub(f2,R)), X_3=>(sub(f3,R))});
	rank rMQ1
	rank EqQ1
	-----
	Q2=convexHull matrix{{0,0,1},{0,1,0}}
	gQ2 = teToricRationalMap (L,Q2)
	rMQ2=representationMatrix (gQ2,16);
	R=QQ[s,t,X_0..X_3];
	EqQ2=sub(sub(rMQ2,R), {X_0=>(sub(f0,R)), X_1=>(sub(f1,R)), X_2=>(sub(f2,R)), X_3=>(sub(f3,R))});
	rank rMQ2
	rank EqQ2
	-----
	Q3=convexHull matrix{{0,0,1,1},{0,1,0,1}}
	gQ3 = teToricRationalMap (L,Q3)
	rMQ3=representationMatrix (gQ3,5);
	R=QQ[s,t,X_0..X_3];
	EqQ3=sub(sub(rMQ3,R), {X_0=>(sub(f0,R)), X_1=>(sub(f1,R)), X_2=>(sub(f2,R)), X_3=>(sub(f3,R))});
	rank rMQ3
	rank EqQ3

-----------------------------------------------
-- Example ninth example (2,3)  
-----------------------------------------------
	 S = QQ[s,t];
	f0 = 3*s^2*t^2-2*s*t^3-s^2*t+s*t^2-3*s*t-t^2+4*t-1
	f1 = 3*s^2*t^2-s^2*t-3*s*t^2-s*t+t^2+1
	f2 = 2*s^2*t^3-3*s^2*t^2-s^2*t+s*t^2+3*s*t-3*t^2+2*t-1
	f3 = 2*s^2*t^3-3*s^2*t^2-2*s*t^3+s^2*t+5*s*t^2-3*s*t-3*t^2+4*t-1

nu=1
-----------------------------------------------
-- Example ninth-bis example (2,3)
-----------------------------------------------
	 S = QQ[s,t];
	f0 = 3*s^2*t^2-2*s*t^3-s^2*t+s*t^2-3*s*t-t^2+4*t-1;
	f1 = 3*s^2*t^2-s^2*t-3*s*t^2-s*t+t^2+1+t^3+s^2;
	f2 = 2*s^2*t^3-3*s^2*t^2-s^2*t+s*t^2+3*s*t-3*t^2+2*t-1;
	f3 = 2*s^2*t^3-3*s^2*t^2-2*s*t^3+s^2*t+5*s*t^2-3*s*t-3*t^2+4*t-1;

	nu=2;
	L={f0,f1,f2,f3};
	rM=representationMatrix (teToricRationalMap ((newToricEmbedding L),L),nu);
	R=QQ[s,t,X_0..X_3];
	Eq=sub(sub(rM,R), {X_0=>(sub(f0,R)), X_1=>(sub(f1,R)), X_2=>(sub(f2,R)), X_3=>(sub(f3,R))});
	rank rM
	rank Eq

-----------------------------------------------
-- Example tenth example (3,4)   -- from D'Andrea/Khetan Example 11
-----------------------------------------------
	 S = QQ[s,t];
	f0 = (t+t^2)*(s-1)^2+(-1-s*t-s^2*t)*(t-1)^2+(t+s*t+s*t^2)*(s-1)*(t-1)
	f1 = (t+t^2)*(s-1)^2+(1+s*t-s^2*t)*(t-1)^2+(t+s*t+s*t^2)*(s-1)*(t-1)
	f2 = (-t-t^2)*(s-1)^2+(-1+s*t+s^2*t)*(t-1)^2+(t+s*t+s*t^2)*(s-1)*(t-1)
	f3 = (t-t^2)*(s-1)^2+(-1-s*t+s^2*t)*(t-1)^2+(t+s*t+s*t^2)*(s-1)*(t-1)

-----------------------------------------------
-- Example eleventh example (3,4)  -- from D'Andrea/Khetan Example 10 and I S SAC paper
-----------------------------------------------
	 S = QQ[s,t];
	f0 = (t+t^2)*(s-1)^2+(-1-s*t-s^2*t)*(t-1)^2
	f1 = (t+t^2)*(s-1)^2+(1+s*t-s^2*t)*(t-1)^2
	f2 = (-t-t^2)*(s-1)^2+(-1+s*t+s^2*t)*(t-1)^2
	f3 = (t-t^2)*(s-1)^2+(-1-s*t+s^2*t)*(t-1)^2


-----------------------------------------------
-- Example twelfth example (6,9) 
-----------------------------------------------
	 S = QQ[s,t];
	f0 = (s^2+t^2)*t^6*s^4+(-1-s^3*t^4-s^4*t^4)*(t-1)^5*(s^2-1)
	f1 = (s^2+t^2)*t^6*s^4+(1+s^3*t^4-s^4*t^4)*(t-1)^5*(s^2-1)
	f2 = (-s^2-t^2)*t^6*s^4+(-1+s^3*t^4+s^4*t^4)*(t-1)^5*(s^2-1)
	f3 = (s^2-t^2)*t^6*s^4+(-1-s^3*t^4+s^4*t^4)*(t-1)^5*(s^2-1)


-----------------------------------------------
-- Example 13th example (4,4)   
-----------------------------------------------
	 S = QQ[s,t];
	f0 = (s^2+t^2)*t*s^2+(-1-s*t^2-s^2*t^2)*(t-1)^2*(s^2-1)
	f1 = (s^2+t^2)*t*s^2+(1+s*t^2-s^2*t^2)*(t-1)^2*(s^2-1)
	f2 = (-s^2-t^2)*t*s^2+(-1+s*t^2+s^2*t^2)*(t-1)^2*(s^2-1)
	f3 = (s^2-t^2)*t*s^2+(-1-s*t^2+s^2*t^2)*(t-1)^2*(s^2-1)


-----------------------------------------------
-- Example 14th example (6,6)  
-----------------------------------------------
	 S = QQ[s,t];
	f0 = (-s^3*t+s^6)*(t-1)^4+(1+s^2*t^5-s^2*t^6)*(s-1)^2
	f1 = (s^3*t-s^6)*(t-1)^4+(1+s^2*t^5+s^2*t^6)*(s-1)^2
	f2 = (s^3*t+s^6)*(t-1)^4+(-1+s^2*t^5+s^2*t^6)*(s-1)^2
	f3 = (-s^3*t-s^6)*(t-1)^4+(1-s^2*t^5+s^2*t^6)*(s-1)^2


-----------------------------------------------
-- Example 15th example (6,6) 
-----------------------------------------------
	 S = QQ[s,t];
	f0 = (-s^3*t^5+s^6^6)*(t-1)^4+(1+s^2*t^7-s^2*t^10)*(s-1)^2
	f1 = (s^3*t^5-s^6^6)*(t-1)^4+(1+s^2*t^7+s^2*t^10)*(s-1)^2
	f2 = (s^3*t^5+s^6^6)*(t-1)^4+(-1+s^2*t^7+s^2*t^10)*(s-1)^2
	f3 = (-s^3*t^5-s^6^6)*(t-1)^4+(1-s^2*t^7+s^2*t^10)*(s-1)^2


-----------------------------------------------
-- Example 16th example (4,8)
-----------------------------------------------
	 S = QQ[s,t];
	f0 = (-2*s^3*t^5+s*t^3)*(s+1)*(t-1)^3+(-s^2*t^3+3)*(t+1)^4
	f1 = (2*s^3*t^5-s*t^3)*(s+1)*(t-1)^3+(s^2*t^3+3)*(t+1)^4
	f2 = (-2*s^3*t^5-s*t^3)*(s+1)*(t-1)^3+(s^2*t^3-3)*(t+1)^4
	f3 = (2*s^3*t^5+s*t^3)*(s+1)*(t-1)^3+(-s^2*t^3-3)*(t+1)^4


-----------------------------------------------
-- Example 17th example (4,8) 
-----------------------------------------------
	 S = QQ[s,t];
	f0 = (-s^4*t^6+s^2*t^2)*(t-1)^2+(-s^2*t^3+2*s)*(t+1)^4
	f1 = (s^4*t^6-s^2*t^2)*(t-1)^2+(s^2*t^3+2*s)*(t+1)^4
	f2 = (-s^4*t^6-s^2*t^2)*(t-1)^2+(s^2*t^3-2*s)*(t+1)^4
	f3 = (s^4*t^6+s^2*t^2)*(t-1)^2+(-s^2*t^3-2*s)*(t+1)^4


-----------------------------------------------
-- Example 18th example (4,6)  -- very sparse
-----------------------------------------------
	 S = QQ[s,t];
	f0 = 2*s^4*t^4+s^2*t^3
	f1 = s^2*t^6+2^6
	f2 = s*t^3-3*s^4*t^4
	f3 = s^5*t+5*s^3*t^5


-----------------------------------------------
-- Example 19th example (4,6)  -- very sparse
-----------------------------------------------
	 S = QQ[s,t];
	f0 = 2*s*t^3+s*t
	f1 = s*t^3+2
	f2 = s*t^2-3*s^2*t^2
	f3 = s*t+5*s^2*t^2



-----------------------------------------------
-- Example 20th example (2,6)  -- very sparse
-----------------------------------------------
	 S  =  QQ[s,t];
	f0  =  2^6+s^2*t^6
	f1  =  s*t^6+2^6
	f2  =  s*t^5-3*s*t^3
	f3  =  s*t^4+5*s^2*t^6
-----------------------------------------------
-- Example 21st example (4,2)  -- very sparse
-----------------------------------------------

	 S = QQ[s,t];
	f0 = 2*s+s^2*t
	f1 = s^4*t^2+2*s
	f2 = s^2*t-3*t
	f3 = s*t+5*s^4*t^2


-----------------------------------------------
-- Example 22nd example (3,3)  -- very sparse
-----------------------------------------------

	 S = QQ[s,t];
	f0 = 2*s+s^2*t^2
	f1 = s^3*t^3+2*s
	f2 = s^2*t^2-3*s*t
	f3 = t+5*s^3*t^3

-----------------------------------------------
-- Example 23st (?) example (2,2)  -- I S SAC example (almost dense)
-----------------------------------------------
	 S  =  QQ[s,t];
	f0  =  s^2*t
	f1  =  t+s^2*t
	f2  =  t^2+s
	f3  =  s^2+s^2*t^2


-----------------------------------------------
-- Ex2 from PerezDiaz (8,4)  --  non-proper
-----------------------------------------------
	 S = QQ[s,t];
	f0 = (t^4+2*t^2+5)*(s^4+1)
	f1 = (s^4*t^4+2*s^4*t^2+5*s^4+2*t^4+4*t^2+11)*(s^4+1)
	f2 = (s^4*t^4+2*s^4*t^2+5*s^4+t^4+2*t^2+6)
	f3 = -(s^4*t^4+2*s^4*t^2+5*s^4+t^4+2*t^2+3)*(s^4+1)


-----------------------------------------------
-- Ex2 from PerezDiaz (2,2)  --  proper 
-----------------------------------------------
	 S = QQ[s,t];
	f0 = (t-5)*(s-1)
	f1 = -(11+s*t-5*s-2*t)*(s-1)
	f2 = (6-t-5*s+s*t)
	f3 = (-t+s*t-5*s+3)*(s-1)