---------------------------------------------------------------------------
-- 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,
	degreeImplicitEq,
	ToricEmbedding,
	newToricEmbedding,
	teToricRing,
	getPolytope,
	teToricRationalMap,
	representationMatrix,
	implicitEq,
	MaxMinor,
	polynomialsToPolytope,
	relationNandP
}


---------------------------------------------------------------------------
needsPackage "Polyhedra"
---------------------------------------------------------------------------

ToricEmbedding = new Type of MutableHashTable
globalAssignment ToricEmbedding

newToricEmbedding = method()

newToricEmbedding (Polyhedron) := (aPolyhedron) -> (
	lat := latticePoints aPolyhedron;
	latLen := length(lat) - 1; -- for indexing
	new ToricEmbedding from {
		tePolyhedron => aPolyhedron,
		teLattice => lat,
		teLatticeLen => latLen,
		teRing => QQ[s,t,u,v,w,T_0..T_latLen,MonomialOrder=>Eliminate 5]
	}
)

newToricEmbedding (List) := (polynomialList) -> (
	newToricEmbedding(polynomialsToPolytope (polynomialList))
)

polynomialsToPolytope = method();

polynomialsToPolytope (List) := (polynomialList) -> (
	aPolytope := newtonPolytope polynomialList_0;
	scan(polynomialList, i -> aPolytope = convexHull(aPolytope, newtonPolytope i));
	return aPolytope;
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------

teMaxDegree = method()

teMaxDegree (ToricEmbedding, ZZ) := (te, indexVar) -> (
 	max apply(0..te#teLatticeLen, i -> (te#teLattice_i_0_indexVar))
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------

teBigradedMonomials = method()

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 | {(s^(te#teLattice_i_0_0)) * (t^(te#teLattice_i_0_1)) * (u^(maxDegreeS - (te#teLattice_i_0_0))) * (v^(maxDegreeT - (te#teLattice_i_0_1)))};
	);
	bigradedMonomials
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------

teToricRing = method()

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 - w * 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;
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------

getPolytope = method();

getPolytope(ToricEmbedding) := (te) -> (
	return te#tePolyhedron;
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------
-- Polytope N based on polytope P (where some homothety of P contains N)

-- 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();

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;
); 

homothety = method();

homothety (Polyhedron, Polyhedron) := (N,P) -> (
	latticeP := latticePoints P;
	h := 1;
	while not contains(convexHull(h * latticeP), N) do (
		h = h + 1;
	);
	return h;
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------

teToricRationalMap = method();

teToricRationalMap(List) := (polynomialList) -> (
	te := newToricEmbedding(polynomialList);

	use te#teRing;
	torVar := teToricRing(te);

	-- substitute monomials by variables Ti interpreted in the ring te#teRing
	-- OPEN ISSUE: how to do this if the toric variety is not defined by the newtonPolytope of polynomialList ?
	(M,C) := coefficients( sub( matrix {polynomialList}, te#teRing), Variables=>{s, t}, 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);
);

---------------------------------------------------------------------------

teToricRationalMap(List, Polyhedron) := (polynomialList, P) -> (
	te := newToricEmbedding(P);
------------------
	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
	-- OPEN ISSUE: how to do this if the toric variety is not defined by the newtonPolytope of polynomialList ?
	(M,C) := coefficients( sub( matrix {polynomialList}, te#teRing), Variables=>{s, t}, 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);
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------

matrixToList = method();

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;
)

---------------------------------------------------------------------------
---------------------------------------------------------------------------

representationMatrix = method();

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;
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------
implicitEq = method();

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)
);

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)
);

---------------------------------------------------------------------------
---------------------------------------------------------------------------
isGoodDegree = method();

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)
);

isGoodDegree(List, ZZ) := (polynomialList, nu) -> (
	return isGoodDegree (matrix {polynomialList},nu)
);
---------------------------------------------------------------------------
---------------------------------------------------------------------------
-- We compute here the degree of the MacRae's invariant = degree of implicit equation

degreeImplicitEq = method();

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;

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

--degreeImplicitEq(List, ZZ) := (polynomialList, nu) -> (
--	return degreeImplicitEq (matrix {polynomialList},nu)
--);

---------------------------------------------------------------------------
---------------------------------------------------------------------------

MaxCol = method();

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
};

---------------------------------------------------------------------------
---------------------------------------------------------------------------

MaxMinor = method();

MaxMinor (Matrix) := (M) -> {
	--Apply two times MaxCol to obtain a maximal minor
	--This is not very efficient...
	return (transpose MaxCol(transpose MaxCol(M)) )
};

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). 

///


document {
	Key => "getPolytope",
	Headline => "returns the Polytope associated to the ToricEmbedding",
	PARA {" INPUT: ToricEmbedding"},
	PARA {" OUTPUT: Polyhedron"},
	PARA {" "},

	EXAMPLE {" 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; ",
		" l = {f0,f1,f2,f3}; ",
		" TOR = newToricEmbedding(l); ",
		" getPolytope TOR "}
}

document {
	Key => "isGoodDegree",
	Headline => "verifies if the Z-complex is acyclic in the given degree",
	PARA {" INPUT: polynomialList '{f0,...,fn}', degreeZZ 'nu'"},
	PARA {" OUTPUT: A boolean 'True' or 'False'"},
	PARA {" The integer 'nu' needs to be a 'good degree' for the first parameter '{f0,...,fn}' that can be verified by doing isGoodDegree(polinomialList,nu) "},
	PARA { "isGoodDegree({f0,...,fn}, nu) verifies if the approximation complex Z associated to the polynomials given is acyclic in degree 'nu' "},
	PARA { "Given a list of polynomials {f_0,...,f_n}, 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+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; ",
		"isGoodDegree ({f1,f2,f3,f4},2)" }
}

document {
	Key => "degreeImplicitEq",
	Headline => "computes the degree of det((Z)_nu)",
	PARA {" INPUT: polynomialList '{f0,...,fn}', degreeZZ 'nu'"},
	PARA {" OUTPUT: an integer number in ZZ"},
	PARA {" The integer 'nu' needs to be a 'good degree' for the first parameter '{f0,...,fn}' that can be verified by doing isGoodDegree(polinomialList,nu) "},
	PARA { "degreeImplicitEq({f0,...,fn}, nu) computes the degree of the gcd of the maximal minors of the matrix representation of {f0,...,fn} in degree 'nu', equivalently, the degree of determinant of the Z-complex associated to {f0,...,fn} in degree 'nu'. This is:"},
	PARA { "degreeImplicitEq({f0,...,fn}, nu) = degree(representationMatrix({f0,...,fn},nu) that computes 'deg(det((Z(f0,...,fn))_nu))'"},
	PARA { "Given a list of polynomials {f_0,...,f_n}, it computes the degree of the gcd of the maximal minors of the right-most map of the strand of degree 'nu'."},
	EXAMPLE {" 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; ",
		"degreeImplicitEq ({f1,f2,f3,f4},2)" },
	SeeAlso => { representationMatrix }
}

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

document {
	Key => "newToricEmbedding",
	Headline => "constructs a ToricEmbedding",
	PARA {" INPUT: a polynomial list '(f0,...,f3)'"},
	PARA {" INPUT: a Polyhedron"},
	PARA {" OUTPUT: a ToricEmbedding instance"},
	PARA {" Given a polynedron P 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 {" Given a list of four polynomials {f_0..f_3} in a polynomial ring in two variables, namely QQ[s,t] we can compute its Newton polytope, namely, the convexHull of the union of the Newton polytope of all the polynomial in the list, with the method 'polynomialsToPolytope'. Take N= polynomialsToPolytope ({f0,...,f3}). Then teToricEmbedding ({f0,...,f3}) can be defined as ToricEmbedding (N)."},
	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), 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+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; ",
		" l = {f0,f1,f2,f3};",
		" teFromPolyList = newToricEmbedding(l);"  },
	PARA{},
	EXAMPLE {" P = convexHull(matrix{{1,0,0},{0,1,0}});",
		" teFromPolyhedron = newToricEmbedding(P);"},
	SeeAlso => {teToricRing, teToricRationalMap, polynomialsToPolytope}
}

document {
	Key => "teToricRing",
	Headline => "Returns the coordinate ring of the toric variety",
	PARA {" INPUT: ToricEmbedding"},
	PARA {" OUTPUT: Ring (the coordinate ring of the toric variety)"},
	PARA {" Given a polynedron P we construc 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 teToricRing gives coordinate ring QQ[T_0..T_p]/J."},
	EXAMPLE {" 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; ",
		" l = {f0,f1,f2,f3};",
		" TE = newToricEmbedding(l);",
		" V = teToricRing(TE);"},
	SeeAlso => {newToricEmbedding, teToricRationalMap }
}

document {
	Key => "teToricRationalMap",
	Headline => "",
	PARA {" INPUT: a polynomial list '(f0,...,f3)'"},
	PARA {" INPUT: a polynomial list '(f0,...,f3)' and a Polyhedron P" },
	PARA {" OUTPUT: a 1-row matrix with 4 entries"},
	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"},
	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 { "" },
	SeeAlso => {newToricEmbedding, teToricRing, polynomialsToPolytope}
}

document {
	Key => "representationMatrix",
	Headline => "computes the right-most map of the Z-complex in degree nu)",
	PARA {" INPUT: polynomialList '{f0,...,fn}', degreeZZ 'nu'"},
	PARA {" INPUT: matrix of homogeneous polynomials 'matrix{{f0,...,fn}}', degreeZZ 'nu'"},
	PARA {" OUTPUT: a matrix of linear functions in QQ[X_0,...,X_n]"},
	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,...,fn}' that can be verified by doing isGoodDegree(polinomialList,nu) "},
	PARA { "representationMatrix(teToricRationalMap{f0,...,fn},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"},
	EXAMPLE {" 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; ",
		"representationMatrix (teToricRationalMap{f1,f2,f3,f4},2)" },
	SeeAlso => { implicitEq }
}

document {
	Key => "implicitEq",
	Headline => "computes the gcd of the right-most map of the Z-complex in degree nu)",
	PARA {" INPUT: polynomialList '{f0,...,fn}', degreeList 'nu'"},
	PARA {" OUTPUT: An homogeneous polynomial in QQ[X_0,...,X_n]"},
	PARA {" The integer 'nu' needs to be a 'good degree' for the first parameter '{f0,...,fn}' that can be verified by doing isGoodDegree(polinomialList,nu) "},
	PARA { "implicitEq({f0,...,fn},nu) computes the determinant of the maximal minors of representationMatrix({f0,...,fn},nu). Equivalently, it computes the determinant of the Z-complex in degree 'nu'. This is:"},
	PARA { "implicitEq({f0,...,fn}, nu)=det((Z(f0,...,fn))_nu)"},
	PARA { "We also have that: degree(implicitEq({f0,...,fn}, nu))=degreeImplicitEq({f0,...,fn}, nu)"},
	EXAMPLE {" 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; ",
		"implicitEq ({f1,f2,f3,f4},2)" },
	PARA { "For the programmer: it computes the determinant of a maximal minor of 'representationMatrix (teToricRationalMap{f1,f2,f3,f4},2)' using the function 'MaxMinor'"},
	SeeAlso => { representationMatrix, degreeImplicitEq, MaxMinor}
}

document {
	Key => "MaxMinor",
	Headline => "Returns a maximal minor of the matrix of full rank.",
	PARA {" INPUT: m x n Matrix of rank r"},
	PARA {" OUTPUT: r x r full rank Matrix"},
	EXAMPLE { " M = matrix {{1,2,3},{1,2,3},{4,5,6},{4,5,6}}",
		    " MaxMinor M;"},
	SeeAlso => {implicitEq }
}

document {
	Key => "polynomialsToPolytope",
	Headline => "Return the convexHull of the union of the Newton polytope of all the polynomial in the list",
	PARA {" INPUT: List"},
	PARA {" OUTPUT: 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+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; ",
		" l = {f0,f1,f2,f3};",
		" polynomialsToPolytope(l);" },
	SeeAlso => {teToricRing, teToricRationalMap, representationMatrix}
}

document {
	Key => "relationNandP",
	Headline => "Returns the list of coordinates of N based on the smallest homothety of P containing N.",
	PARA {" INPUT: (Polyhedron N, Polyhedron P)"},
	PARA {" OUTPUT: 1-column 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+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; ",
		" 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}
}
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^1
	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;

needsPackage "Polyhedra";
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,1,1},{0,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)  AGREGAR TERMINOS AL ANTERIOR
-----------------------------------------------
	 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;
load "./Desktop/Implicitization.m2"

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)  -- Interesting: saturation index not verified
-- \\!***!// does not stops after 5 minutes \\!***!//
-----------------------------------------------
	 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)  -- 
-- \\!***!// *** buffer overflow detected ***: M2 terminated
-----------------------------------------------
	 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)  -- 
-- \\!***!// does not stops after 5 minutes \\!***!//
-----------------------------------------------
	 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)  -- strange: Euler char 0 in nu = 1,2 but not imp. eq.
-- \\!***!// does not stops after 5 minutes \\!***!//
-----------------------------------------------
	 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)  -- strange: Euler char 0 in nu = 1,2 but not imp. eq.
-----------------------------------------------
	 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+s^2*t^6
	f1  =  s*t^6+2
	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)