version="$Id: polymake.lib,v 1.13 2009-04-06 12:39:02 seelisch Exp $"; category="Tropical Geometry"; info=" LIBRARY: polymake.lib Computations with polytopes and fans, interface to polymake and TOPCOM AUTHOR: Thomas Markwig, email: keilen@mathematik.uni-kl.de WARNING: Most procedures will not work unless polymake or topcom is installed and if so, they will only work with the operating system LINUX! For more detailed information see IMPORTANT NOTE respectively consult the help string of the procedures. IMPORTANT NOTE: Even though this is a Singular library for computing polytopes and fans such as the Newton polytope or the Groebner fan of a polynomial, most of the hard computations are NOT done by Singular but by the program @* - polymake by Ewgenij Gawrilow, TU Berlin and Michael Joswig, TU Darmstadt @* (see http://www.math.tu-berlin.de/polymake/), @* respectively (only in the procedure triangularions) by the program @* - topcom by Joerg Rambau, Universitaet Bayreuth (see @* http://www.uni-bayreuth.de/departments/wirtschaftsmathematik/rambau/TOPCOM); @* this library should rather be seen as an interface which allows to use a (very limited) number of options which polymake respectively topcom offers to compute with polytopes and fans and to make the results available in Singular for further computations; moreover, the user familiar with Singular does not have to learn the syntax of polymake or topcom, if the options offered here are sufficient for his purposes. @* Note, though, that the procedures concerned with planar polygons are independent of both, polymake and topcom. PROCEDURES USING POLYMAKE: polymakePolytope computes the vertices of a polytope using polymake newtonPolytope computes the Newton polytope of a polynomial newtonPolytopeLP computes the lattice points of the Newton polytope normalFan computes the normal fan of a polytope groebnerFan computes the Groebner fan of a polynomial intmatToPolymake transforms an integer matrix into polymake format polymakeToIntmat transforms polymake output into an integer matrix PROCEDURES USING TOPCOM: triangulations computes all triangulations of a marked polytope secondaryPolytope computes the secondary polytope of a marked polytope PROCEDURES USING POLYMAKE AND TOPCOM: secondaryFan computes the secondary fan of a marked polytope PROCEDURES CONERNED WITH PLANAR POLYGONS: cycleLength computes the cycleLength of cycle splitPolygon splits a marked polygon into vertices, facets, interior points eta computes the eta-vector of a triangulation findOrientedBoundary computes the boundary of a convex hull cyclePoints computes lattice points connected to some lattice point latticeArea computes the lattice area of a polygon picksFormula computes the ingrediants of Pick's formula for a polygon ellipticNF computes the normal form of an elliptic polygon ellipticNFDB displays the 16 normal forms of elliptic polygons AUXILARY PROCEDURES: polymakeKeepTmpFiles determines if the files created in /tmp should be kept KEYWORDS: polytope; fan; secondary fan; secondary polytope; polymake; Newton polytope; Groebner fan "; //////////////////////////////////////////////////////////////////////////////// /// Auxilary Static Procedures in this Library //////////////////////////////////////////////////////////////////////////////// /// - scalarproduct /// - intmatcoldelete /// - intmatconcat /// - sortlist /// - minInList /// - stringdelete /// - abs /// - commondenominator /// - maxPosInIntvec /// - maxPosInIntmat /// - sortintvec /// - matrixtointmat //////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// LIB "poly.lib"; LIB "linalg.lib"; LIB "random.lib"; //////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////// /// PROCEDURES USING POLYMAKE ///////////////////////////////////////////////////////////////////////////// proc polymakePolytope (intmat polytope,list #) "USAGE: polymakePolytope(polytope[,#]); polytope list, # string ASSUME: each row of polytope gives the coordinates of a lattice point of a polytope with their affine coordinates as given by the output of secondaryPolytope PURPOSE: the procedure calls polymake to compute the vertices of the polytope as well as its dimension and information on its facets RETURN: list L with four entries @* L[1] : an integer matrix whose rows are the coordinates of vertices of the polytope @* L[2] : the dimension of the polytope @* L[3] : a list whose i-th entry explains to which vertices the ith vertex of the Newton polytope is connected -- i.e. L[3][i] is an integer vector and an entry k in there means that the vertex L[1][i] is connected to the vertex L[1][k] @* L[4] : an integer matrix whose rows mulitplied by (1,var(1),...,var(nvar)) give a linear system of equations describing the affine hull of the polytope, i.e. the smallest affine space containing the polytope NOTE: - for its computations the procedure calls the program polymake by Ewgenij Gawrilow, TU Berlin and Michael Joswig, TU Darmstadt; it therefore is necessary that this program is installed in order to use this procedure; see http://www.math.tu-berlin.de/polymake/ @* - note that in the vertex edge graph we have changed the polymake convention which starts indexing its vertices by zero while we start with one ! @* - the procedure creates the file /tmp/polytope.polymake which contains the polytope in polymake format; if you wish to use this for further computations with polymake, you have to use the procedure polymakeKeepTmpFiles in before @* - moreover, the procedure creates the file /tmp/polytope.output which it deletes again before ending @* - it is possible to provide an optional second argument as string which then will be used instead of 'polytope' in the name of the polymake output file EXAMPLE: example polymakePolytope; shows an example" { // the header for the file secendarypolytope.polymake string sp="_application polytope _version 2.2 _type RationalPolytope POINTS "; int i,j; // set the name for the polymake output file if (size(#)>0) { if (typeof(#[1])=="string") { string dateiname=#[1]; } else { string dateiname="polytope"; } } else { string dateiname="polytope"; } // create the lattice point list for polymake sp=sp+intmatToPolymake(polytope,"points"); // initialise dateiname.polymake and compute the vertices write(":w /tmp/"+dateiname+".polymake",sp); system("sh","cd /tmp; polymake "+dateiname+".polymake VERTICES > "+dateiname+".output"); string vertices=read("/tmp/"+dateiname+".output"); system("sh","/bin/rm /tmp/"+dateiname+".output"); intmat np=polymakeToIntmat(vertices,"affine"); // compute the dimension system("sh","cd /tmp; polymake "+dateiname+".polymake DIM > "+dateiname+".output"); string pdim=read("/tmp/"+dateiname+".output"); system("sh","/bin/rm /tmp/"+dateiname+".output"); pdim=pdim[5,size(pdim)-6]; execute("int nd="+pdim+";"); // compute the vertex-edge graph system("sh","cd /tmp; polymake "+dateiname+".polymake GRAPH > "+dateiname+".output"); string vertexedgegraph=read("/tmp/"+dateiname+".output"); system("sh","/bin/rm /tmp/"+dateiname+".output"); vertexedgegraph=vertexedgegraph[7,size(vertexedgegraph)-8]; string newveg; for (i=1;i<=size(vertexedgegraph);i++) { if (vertexedgegraph[i]=="{") { newveg=newveg+"intvec("; } else { if (vertexedgegraph[i]=="}") { newveg=newveg+"),"; } else { if (vertexedgegraph[i]==" ") { newveg=newveg+","; } else { newveg=newveg+vertexedgegraph[i]; } } } } newveg=newveg[1,size(newveg)-1]; execute("list nveg="+newveg+";"); // raise each entry in nveg by one for (i=1;i<=size(nveg);i++) { for (j=1;j<=size(nveg[i]);j++) { nveg[i][j]=nveg[i][j]+1; } } // compute the affine hull system("sh","cd /tmp; polymake "+dateiname+".polymake AFFINE_HULL > "+dateiname+".output"); string equations=read("/tmp/"+dateiname+".output"); system("sh","/bin/rm /tmp/"+dateiname+".output"); if (size(equations)>14) { intmat neq=polymakeToIntmat(equations,"cleardenom"); } else { intmat neq[1][ncols(polytope)+1]; } // delete the tmp-files, if polymakekeeptmpfiles is not set if (defined(polymakekeeptmpfiles)==0) { system("sh","/bin/rm /tmp/"+dateiname+".polymake"); } // return the files return(list(np,nd,nveg,neq)); } example { "EXAMPLE:"; echo=2; // the lattice points of the unit square in the plane list points=intvec(0,0),intvec(0,1),intvec(1,0),intvec(1,1); // the secondary polytope of this lattice point configuration is computed intmat secpoly=secondaryPolytope(points)[1]; list np=polymakePolytope(secpoly); // the vertices of the secondary polytope are: np[1]; // its dimension is np[2]; // np[3] contains information how the vertices are connected to each other, // e.g. the first vertex (number 0) is connected to the second one np[3][1]; // the affine hull has the equation ring r=0,x(1..4),dp; matrix M[5][1]=1,x(1),x(2),x(3),x(4); np[4]*M; } ///////////////////////////////////////////////////////////////////////////// proc newtonPolytope (poly f,list #) "USAGE: newtonPolytope(f[,#]); f poly, # string RETURN: list, L with four entries @* L[1] : an integer matrix whose rows are the coordinates of vertices of the Newton polytope of f @* L[2] : the dimension of the Newton polytope of f @* L[3] : a list whose ith entry explains to which vertices the ith vertex of the Newton polytope is connected -- i.e. L[3][i] is an integer vector and an entry k in there means that the vertex L[1][i] is connected to the vertex L[1][k] @* L[4] : an integer matrix whose rows mulitplied by (1,var(1),...,var(nvar)) give a linear system of equations describing the affine hull of the Newton polytope, i.e. the smallest affine space containing the Newton polytope NOTE: - if we replace the first column of L[4] by zeros, i.e. if we move the affine hull to the origin, then we get the equations for the orthogonal comploment of the linearity space of the normal fan dual to the Newton polytope, i.e. we get the EQUATIONS that we need as input for polymake when computing the normal fan @* - the procedure calls for its computation polymake by Ewgenij Gawrilow, TU Berlin and Michael Joswig, so it only works if polymake is installed; see http://www.math.tu-berlin.de/polymake/ @* - the procedure creates the file /tmp/newtonPolytope.polymake which contains the polytope in polymake format and which can be used for further computations with polymake @* - moreover, the procedure creates the file /tmp/newtonPolytope.output which it deletes again before ending @* - it is possible to give as an optional second argument as string which then will be used instead of 'newtonPolytope' in the name of the polymake output file EXAMPLE: example newtonPolytope; shows an example" { int i,j; // compute the list of exponent vectors of the polynomial, // which are the lattice points // whose convex hull is the Newton polytope of f intmat exponents[size(f)][nvars(basering)]; while (f!=0) { i++; exponents[i,1..nvars(basering)]=leadexp(f); f=f-lead(f); } if (size(#)==0) { #[1]="newtonPolytope"; } // call polymakePolytope with exponents return(polymakePolytope(exponents,#)); } example { "EXAMPLE:"; echo=2; ring r=0,(x,y,z),dp; matrix M[4][1]=1,x,y,z; poly f=y3+x2+xy+2xz+yz+z2+1; // the Newton polytope of f is list np=newtonPolytope(f); // the vertices of the Newton polytope are: np[1]; // its dimension is np[2]; // np[3] contains information how the vertices are connected to each other, // e.g. the first vertex (number 0) is connected to the second, third and // fourth vertex np[3][1]; ////////////////////////// f=x2-y3; // the Newton polytope of f is np=newtonPolytope(f); // the vertices of the Newton polytope are: np[1]; // its dimension is np[2]; // the Newton polytope is contained in the affine space given // by the equations np[4]*M; } ///////////////////////////////////////////////////////////////////////////// proc newtonPolytopeLP (poly f) "USAGE: newtonPolytopeLP(f); f poly RETURN: list, the exponent vectors of the monomials occuring in f, i.e. the lattice points of the Newton polytope of f EXAMPLE: example normalFan; shows an example" { list np; int i=1; while (f!=0) { np[i]=leadexp(f); f=f-lead(f); i++; } return(np); } example { "EXAMPLE:"; echo=2; ring r=0,(x,y,z),dp; poly f=y3+x2+xy+2xz+yz+z2+1; // the lattice points of the Newton polytope of f are newtonPolytopeLP(f); } ///////////////////////////////////////////////////////////////////////////// proc normalFan (intmat vertices,intmat affinehull,list graph,int er,list #) "USAGE: normalFan (vert,aff,graph,rays,[,#]); vert,aff intmat, graph list, rays int, # string ASSUME: - vert is an integer matrix whose rows are the coordinate of the vertices of a convex lattice polytope; @* - aff describes the affine hull of this polytope, i.e. the smallest affine space containing it, in the following sense: denote by n the number of columns of vert, then multiply aff by (1,x(1),...,x(n)) and set the resulting terms to zero in order to get the equations for the affine hull; @* - the ith entry of graph is an integer vector describing to which vertices the ith vertex is connected, i.e. a k as entry means that the vertex vert[i] is connected to vert[k]; @* - the integer rays is either one (if the extreme rays should be computed) or zero (otherwise) RETURN: list, the ith entry of L[1] contains information about the cone in the normal fan dual to the ith vertex of the polytope @* L[1][i][1] = integer matrix representing the inequalities which describe the cone dual to the ith vertex @* L[1][i][2] = a list which contains the inequalities represented by L[i][1] as a list of strings, where we use the variables x(1),...,x(n) @* L[1][i][3] = only present if 'er' is set to 1; in that case it is an interger matrix whose rows are the extreme rays of the cone @* L[2] = is an integer matrix whose rows span the linearity space of the fan, i.e. the linear space which is contained in each cone NOTE: - the procedure calls for its computation polymake by Ewgenij Gawrilow, TU Berlin and Michael Joswig, so it only works if polymake is installed; see http://www.math.tu-berlin.de/polymake/ @* - in the optional argument # it is possible to hand over other names for the variables to be used -- be carful, the format must be correct and that is not tested, e.g. if you want the variable names to be u00,u10,u01,u11 then you must hand over the string u11,u10,u01,u11 EXAMPLE: example normalFan; shows an example" { list ineq; // stores the inequalities of the cones int i,j,k; // we work over the following ring execute("ring ineqring=0,x(1.."+string(ncols(vertices))+"),lp;"); string greatersign=">"; // create the variable names if (size(#)>0) { if (typeof(#[1])=="string") { kill ineqring; execute("ring ineqring=0,("+#[1]+"),lp;"); } if (size(#)>1) { greatersign="<"; } } ////////////////////////////////////////////////////////////////// // Compute first the inequalities of the cones ////////////////////////////////////////////////////////////////// matrix VAR[1][ncols(vertices)]=maxideal(1); matrix EXP[ncols(vertices)][1]; poly p,pl,pr; // consider all vertices of the polytope for (i=1;i<=nrows(vertices);i++) { // first we produce for each vertex in the polytope // the inequalities describing the dual cone in the normal fan list pp; // contain strings representing the inequalities // describing the normal cone intmat ie[size(graph[i])][ncols(vertices)]; // contains the inequalities // as rows // consider all the vertices to which the ith vertex in the // polytope is connected by an edge for (j=1;j<=size(graph[i]);j++) { // produce the vector ie_j pointing from the jth vertex to the ith vertex; // this will be the jth inequality for the cone in the normal fan dual to // the ith vertex ie[j,1..ncols(vertices)]=vertices[i,1..ncols(vertices)]-vertices[graph[i][j],1..ncols(vertices)]; EXP=ie[j,1..ncols(vertices)]; // build a linear polynomial with the entries of ie_j as coefficients p=(VAR*EXP)[1,1]; pl,pr=0,0; // separate the terms with positive coefficients in p from // those with negative coefficients for (k=1;k<=size(p);k++) { if (leadcoef(p[k])<0) { pr=pr-p[k]; } else { pl=pl+p[k]; } } // build the string which represents the jth inequality // for the cone dual to the ith vertex // as polynomial inequality of type string, and store this // in the list pp as jth entry pp[j]=string(pl)+" "+greatersign+" "+string(pr); } // all inequalities for the ith vertex are stored in the list ineq ineq[i]=list(ie,pp); kill ie,pp; // kill certain lists } // remove the first column of affine hull to compute the linearity space intmat linearity=intmatcoldelete(affinehull,1); ////////////////////////////////////////////////////////////////// // Compute next the extreme rays of the cones ////////////////////////////////////////////////////////////////// if (er==1) { list extremerays; // keeps the result string polymake; // keeps polymake output // the header for ineq.polymake string head="_application polytope _version 2.2 _type RationalPolytope INEQUALITIES "; // the tail for both polymake files string tail=" EQUATIONS "; tail=tail+intmatToPolymake(linearity,"rays"); string ungleichungen; // keeps the inequalities for the polymake code intmat M; // the matrix keeping the inequalities // create the file ineq.output write(":w /tmp/ineq.output",""); int dimension; // keeps the dimension of the intersection the // bad cones with the u11tobeseencone for (i=1;i<=size(ineq);i++) { i,". Cone of ",nrows(vertices); // indicate how many // vertices have been dealt with ungleichungen=intmatToPolymake(ineq[i][1],"rays"); // write the inequalities to ineq.polymake and call polymake write(":w /tmp/ineq.polymake",head+ungleichungen+tail); ungleichungen=""; // clear ungleichungen system("sh","cd /tmp; /bin/rm ineq.output; polymake ineq.polymake VERTICES > ineq.output"); // read the result of polymake polymake=read("/tmp/ineq.output"); intmat VERT=polymakeToIntmat(polymake,"affine"); extremerays[i]=VERT; kill VERT; } for (i=1;i<=size(ineq);i++) { ineq[i]=ineq[i]+list(extremerays[i]); } } // delete the tmp-files, if polymakekeeptmpfiles is not set if (defined(polymakekeeptmpfiles)==0) { system("sh","/bin/rm /tmp/ineq.polymake"); system("sh","/bin/rm /tmp/ineq.output"); } // get the linearity space return(list(ineq,linearity)); } example { "EXAMPLE:"; echo=2; ring r=0,(x,y,z),dp; matrix M[4][1]=1,x,y,z; poly f=y3+x2+xy+2xz+yz+z2+1; // the Newton polytope of f is list np=newtonPolytope(f); // the Groebner fan of f, i.e. the normal fan of the Newton polytope list gf=normalFan(np[1],np[4],np[3],1,"x,y,z"); // the number of cones in the Groebner fan of f is: size(gf[1]); // the inequalities of the first cone as matrix are: print(gf[1][1][1]); // the inequalities of the first cone as string are: print(gf[1][1][2]); // the rows of the following matrix are the extreme rays of the first cone: print(gf[1][1][3]); // each cone contains the linearity space spanned by: print(gf[2]); } ///////////////////////////////////////////////////////////////////////////// proc groebnerFan (poly f,list #) "USAGE: groebnerFan(f[,#]); f poly, # string RETURN: list, the ith entry of L[1] contains information about the ith cone in the Groebner fan dual to the ith vertex in the Newton polytope of the f @* L[1][i][1] = integer matrix representing the inequalities which describe the cone @* L[1][i][2] = a list which contains the inequalities represented by L[1][i][1] as a list of strings @* L[1][i][3] = an interger matrix whose rows are the extreme rays of the cone @* L[2] = is an integer matrix whose rows span the linearity space of the fan, i.e. the linear space which is contained in each cone @* L[3] = the Newton polytope of f in the format of the procedure newtonPolytope @* L[4] = integer matrix where each row represents the exponet vector of one monomial occuring in the input polynomial NOTE: - if you have alread computed the Newton polytope of f then you might want to use the procedure normalFan instead in order to avoid doing costly computation twice @* - the procedure calls for its computation polymake by Ewgenij Gawrilow, TU Berlin and Michael Joswig, so it only works if polymake is installed; see http://www.math.tu-berlin.de/polymake/ @* - the procedure creates the file /tmp/newtonPolytope.polymake which contains the Newton polytope of f in polymake format and which can be used for further computations with polymake @* - it is possible to give as an optional second argument as string which then will be used instead of 'newtonPolytope' in the name of the polymake output file EXAMPLE: example groebnerFan; shows an example" { int i,j; // compute the list of exponent vectors of the polynomial, which are // the lattice points whose convex hull is the Newton polytope of f intmat exponents[size(f)][nvars(basering)]; while (f!=0) { i++; exponents[i,1..nvars(basering)]=leadexp(f); f=f-lead(f); } if (size(#)==0) { #[1]="newtonPolytope"; } // call polymakePolytope with exponents list newtonp=polymakePolytope(exponents,"newtonPolytope"); // get the variables as string string variablen; for (i=1;i<=nvars(basering);i++) { variablen=variablen+string(var(i))+","; } variablen=variablen[1,size(variablen)-1]; // call normalFan in order to compute the Groebner fan list gf=normalFan(newtonp[1],newtonp[4],newtonp[3],1,variablen); // append newtonp to gf gf[3]=newtonp; // append the exponent vectors to gf gf[4]=exponents; return(gf); } example { "EXAMPLE:"; echo=2; ring r=0,(x,y,z),dp; matrix M[4][1]=1,x,y,z; poly f=y3+x2+xy+2xz+yz+z2+1; // the Newton polytope of f is list gf=groebnerFan(f); // the exponent vectors of f are ordered as follows gf[4]; // the first cone of the groebner fan has the inequalities gf[1][1][1]; // as a string they look like gf[1][1][2]; // and it has the extreme rays print(gf[1][1][3]); // the linearity space is spanned by print(gf[2]); // the vertices of the Newton polytope are: gf[3][1]; // its dimension is gf[3][2]; // np[3] contains information how the vertices are connected to each other, // e.g. the 1st vertex is connected to the 2nd, 3rd and 4th vertex gf[3][3][1]; } ///////////////////////////////////////////////////////////////////////////// proc intmatToPolymake (intmat M,string art) "USAGE: intmatToPolymake(M,art); M intmat, art string ASSUME: - M is an integer matrix which should be transformed into polymake format; @* - art is one of the following strings: @* + 'rays' : indicating that a first column of 0's should be added @* + 'points' : indicating that a first column of 1's should be added RETURN: string, the matrix is transformed in a string and a first column has been added EXAMPLE: example intmatToPolymake; shows an example" { if (art=="rays") { string anf="0 "; } else { string anf="1 "; } string sp; int i,j; // create the lattice point list for polymake for (i=1;i<=nrows(M);i++) { sp=sp+anf; for (j=1;j<=ncols(M);j++) { sp=sp+string(M[i,j])+" "; if (j==ncols(M)) { sp=sp+" "; } } } return(sp); } example { "EXAMPLE:"; echo=2; intmat M[3][4]=1,2,3,4,5,6,7,8,9,10,11,12; intmatToPolymake(M,"rays"); intmatToPolymake(M,"points"); } ///////////////////////////////////////////////////////////////////////////// proc polymakeToIntmat (string pm,string art) "USAGE: polymakeToIntmat(pm,art); pm, art string ASSUME: pm is the result of calling polymake with one 'argument' like VERTICES, AFFINE_HULL, etc., so that the first row of the string is the name of the corresponding 'argument' and the further rows contain the result which consists of vectors either over the integers or over the rationals RETURN: intmat, the rows of the matrix are basically the vectors in pm from the second row on where each row has been multiplied with the lowest common multiple of the denominators of its entries so as to be an integer matrix; moreover, if art=='affine', then the first column is omitted since we only want affine coordinates EXAMPLE: example polymakeToIntmat; shows an example" { // we need a line break string zeilenumbruch=" "; // remove the 'argment' name, i.e. the first row of pm while (pm[1]!=zeilenumbruch) { pm=stringdelete(pm,1); } pm=stringdelete(pm,1); // find out how many entries each vector has, namely one more // than 'spaces' in a row int i=1; int s=1; int z=1; while (pm[i]!=zeilenumbruch) { if (pm[i]==" ") { s++; } i++; } // if we want to have affine coordinates if (art=="affine") { s--; // then there is one column less // and the entry of the first column (in the first row) has to be removed while (pm[1]!=" ") { pm=stringdelete(pm,1); } pm=stringdelete(pm,1); } // we add two line breaks at the end in order to have this as // a stopping criterion pm=pm+zeilenumbruch+zeilenumbruch; // we now have to work through each row for (i=1;i<=size(pm);i++) { // if there are two consecutive line breaks we are done if ((pm[i]==zeilenumbruch) and (pm[i+1]==zeilenumbruch)) { i=size(pm)+1; } else { // a line break has to be replaced by a comma if (pm[i]==zeilenumbruch) { z++; pm[i]=","; // if we want to have affine coordinates, // then we have to delete the first entry in each row if (art=="affine") { while (pm[i+1]!=" ") { pm=stringdelete(pm,i+1); } pm=stringdelete(pm,i+1); } } // a space has to be replaced by a comma if (pm[i]==" ") { pm[i]=","; } } } // if we have introduced superflous commata at the end, they should be removed while (pm[size(pm)]==",") { pm=stringdelete(pm,size(pm)); } // since the matrix could be over the rationals, // we need a ring with rational coefficients ring zwischering=0,x,lp; // create the matrix with the elements of pm as entries execute("matrix mm["+string(z)+"]["+string(s)+"]="+pm+";"); // transform this into an integer matrix matrix M[1][ncols(mm)]; // takes a row of mm int cm; // takes a lowest common multiple // multiply each row by an integer such that its entries are integers for (int j=1;j<=nrows(mm);j++) { M=mm[j,1..ncols(mm)]; cm=commondenominator(M); for (i=1;i<=ncols(mm);i++) { mm[j,i]=cm*mm[j,i]; } } // transform the matrix mm into an integer matrix execute("intmat im["+string(z)+"]["+string(s)+"]="+string(mm)+";"); return(im); } example { "EXAMPLE:"; echo=2; // this is the usual output of some polymake computation string pm="VERTICES 0 1 3 5/3 1/3 -1 -23/3 -1/3 5/3 1/3 1 0 1 3 -23/3 5/3 1 5/3 1/3 1/3 -1/3 -1 0 1 1 1/3 -1/3 -1 5/3 1/3 -23/3 5/3 3 0 1 1 5/3 -23/3 3 1/3 5/3 -1/3 1/3 -1 0 1 -1 1/3 5/3 3 -1/3 -23/3 1/3 5/3 1 0 1 -1 -1/3 1/3 1 1/3 5/3 5/3 -23/3 3 0 1 -1 1 3 -5 -1 3 -1 1 -1 0 1 -1 -1 -1 -1 1 1 3 3 -5 0 1 -5 3 1 -1 3 -1 1 -1 -1 "; intmat PM=polymakeToIntmat(pm,"affine"); // note that the first column has been removed, since we asked for // affine coordinates, and the denominators have been cleared print(PM); } /////////////////////////////////////////////////////////////////////////////// /// PROCEDURES USING TOPCOM /////////////////////////////////////////////////////////////////////////////// proc triangulations (list polygon) "USAGE: triangulations(polygon); list polygon ASSUME: polygon is a list of integer vectors of the same size representing the affine coordinates of the lattice points PURPOSE: the procedure considers the marked polytope given as the convex hull of the lattice points and with these lattice points as markings; it then computes all possible triangulations of this marked polytope RETURN: list, each entry corresponds to one triangulation and the ith entry is itself a list of integer vectors of size three, where each integer vector defines one triangle in the triangulation by telling which points of the input are the vertices of the triangle NOTE:- the procedure calls for its computations the program points2triangs from the program topcom by Joerg Rambau, Universitaet Bayreuth; it therefore is necessary that this program is installed in order to use this procedure; see @* http://www.uni-bayreuth.de/departments/wirtschaftsmathematik/rambau/TOPCOM @* - the procedure creates the files /tmp/triangulationsinput and /tmp/triangulationsoutput; the former is used as input for points2triangs and the latter is its output containing the triangulations of corresponding to points in the format of points2triangs; if you wish to use this for further computations with topcom, you have to use the procedure polymakeKeepTmpFiles in before @* - note that an integer i in an integer vector representing a triangle refers to the ith lattice point, i.e. polygon[i]; this convention is different from TOPCOM's convention, where i would refer to the i-1st lattice point EXAMPLE: example triangulations; shows an example" { int i,j; // prepare the input for points2triangs by writing the input polygon in the // necessary format string spi="["; for (i=1;i<=size(polygon);i++) { polygon[i][size(polygon[i])+1]=1; spi=spi+"["+string(polygon[i])+"]"; if (i triangulationsoutput"); string p2t=read("/tmp/triangulationsoutput"); // takes result of points2triangs // delete the tmp-files, if polymakekeeptmpfiles is not set if (defined(polymakekeeptmpfiles)==0) { system("sh","cd /tmp; rm -f triangulationsinput; rm -f triangulationsoutput"); } // preprocessing of p2t if points2triangs is version >= 0.15 // brings p2t to the format of version 0.14 string np2t; // takes the triangulations in Singular format for (i=1;i<=size(p2t)-2;i++) { if ((p2t[i]==":") and (p2t[i+1]=="=") and (p2t[i+2]=="[")) { np2t=np2t+p2t[i]+p2t[i+1]; i=i+3; while (p2t[i]!=":") { i=i+1; } } else { if ((p2t[i]=="]") and (p2t[i+1]==";")) { np2t=np2t+p2t[i+1]; i=i+1; } else { np2t=np2t+p2t[i]; } } } if (p2t[size(p2t)-1]=="]") { np2t=np2t+p2t[size(p2t)]; } else { if (np2t[size(np2t)]!=";") { np2t=np2t+p2t[size(p2t)-1]+p2t[size(p2t)]; } } p2t=np2t; np2t=""; // transform the points2triangs output of version 0.14 into Singular format for (i=1;i<=size(p2t);i++) { if (p2t[i]=="=") { np2t=np2t+p2t[i]+"list("; i++; } else { if (p2t[i]!=":") { if ((p2t[i]=="}") and (p2t[i+1]=="}")) { np2t=np2t+"))"; i++; } else { if (p2t[i]=="{") { np2t=np2t+"intvec("; } else { if (p2t[i]=="}") { np2t=np2t+")"; } else { np2t=np2t+p2t[i]; } } } } } } list T; execute(np2t); // raise each index by one for (i=1;i<=size(T);i++) { for (j=1;j<=size(T[i]);j++) { T[i][j]=T[i][j]+1; } } return(T); } example { "EXAMPLE:"; echo=2; // the lattice points of the unit square in the plane list polygon=intvec(0,0),intvec(0,1),intvec(1,0),intvec(1,1); // the triangulations of this lattice point configuration are computed list triang=triangulations(polygon); triang; } ///////////////////////////////////////////////////////////////////////////// proc secondaryPolytope (list polygon,list #) "USAGE: secondaryPolytope(polygon[,#]); list polygon, list # ASSUME: - polygon is a list of integer vectors of the same size representing the affine coordinates of lattice points @* - if the triangulations of the corresponding polygon have already been computed with the procedure triangulations then these can be given as a second (optional) argument in order to avoid doing this computation again PURPOSE: the procedure considers the marked polytope given as the convex hull of the lattice points and with these lattice points as markings; it then computes the lattice points of the secondary polytope given by this marked polytope which correspond to the triangulations computed by the procedure triangulations RETURN: list, say L, such that: @* L[1] = intmat, each row gives the affine coordinates of a lattice point in the secondary polytope given by the marked polytope corresponding to polygon @* L[2] = the list of corresponding triangulations NOTE: if the triangluations are not handed over as optional argument the procedure calls for its computation of these triangulations the program points2triangs from the program topcom by Joerg Rambau, Universitaet Bayreuth; it therefore is necessary that this program is installed in order to use this procedure; see @* http://www.uni-bayreuth.de/departments/wirtschaftsmathematik/rambau/TOPCOM EXAMPLE: example secondaryPolytope; shows an example" { // compute the triangulations of the point configuration with points2triangs if (size(#)==0) { list triangs=triangulations(polygon); } else { list triangs=#; } int i,j,k,l; intmat N[2][2]; // is used to compute areas of triangles intvec vertex; // stores a point in the secondary polytope as // intermediate result int eintrag; int halt; intmat secpoly[size(triangs)][size(polygon)]; // stores all lattice points // of the secondary polytope // consider each triangulation and compute the corresponding point // in the secondary polytope for (i=1;i<=size(triangs);i++) { // for each triangulation we have to compute the coordinates // corresponding to each marked point for (j=1;j<=size(polygon);j++) { eintrag=0; // for each marked point we have to consider all triangles in the // triangulation which involve this particular point for (k=1;k<=size(triangs[i]);k++) { halt=0; for (l=1;(l<=3) and (halt==0);l++) { if (triangs[i][k][l]==j) { halt=1; N[1,1]=polygon[triangs[i][k][3]][1]-polygon[triangs[i][k][1]][1]; N[1,2]=polygon[triangs[i][k][2]][1]-polygon[triangs[i][k][1]][1]; N[2,1]=polygon[triangs[i][k][3]][2]-polygon[triangs[i][k][1]][2]; N[2,2]=polygon[triangs[i][k][2]][2]-polygon[triangs[i][k][1]][2]; eintrag=eintrag+abs(det(N)); } } } vertex[j]=eintrag; } secpoly[i,1..size(polygon)]=vertex; } return(list(secpoly,triangs)); } example { "EXAMPLE:"; echo=2; // the lattice points of the unit square in the plane list polygon=intvec(0,0),intvec(0,1),intvec(1,0),intvec(1,1); // the secondary polytope of this lattice point configuration is computed list secpoly=secondaryPolytope(polygon); // the points in the secondary polytope print(secpoly[1]); // the corresponding triangulations secpoly[2]; } /////////////////////////////////////////////////////////////////////////////// /// PROCEDURES USING POLYMAKE AND TOPCOM /////////////////////////////////////////////////////////////////////////////// proc secondaryFan (list polygon,list #) "USAGE: secondaryFan(polygon[,#]); list polygon, list # ASSUME: - polygon is a list of integer vectors of the same size representing the affine coordinates of lattice points @* - if the triangulations of the corresponding polygon have already been computed with the procedure triangulations then these can be given as a second (optional) argument in order to avoid doing this computation again PURPOSE: the procedure considers the marked polytope given as the convex hull of the lattice points and with these lattice points as markings; it then computes the lattice points of the secondary polytope given by this marked polytope which correspond to the triangulations computed by the procedure triangulations RETURN: list, the ith entry of L[1] contains information about the ith cone in the secondary fan of the polygon, i.e. the cone dual to the ith vertex of the secondary polytope @* L[1][i][1] = integer matrix representing the inequalities which describe the cone dual to the ith vertex @* L[1][i][2] = a list which contains the inequalities represented by L[1][i][1] as a list of strings, where we use the variables x(1),...,x(n) @* L[1][i][3] = only present if 'er' is set to 1; in that case it is an interger matrix whose rows are the extreme rays of the cone @* L[2] = is an integer matrix whose rows span the linearity space of the fan, i.e. the linear space which is contained in each cone @* L[3] = the secondary polytope in the format of the procedure polymakePolytope @* L[4] = the list of triangulations corresponding to the vertices of the secondary polytope NOTE:- the procedure calls for its computation polymake by Ewgenij Gawrilow, TU Berlin and Michael Joswig, so it only works if polymake is installed; see http://www.math.tu-berlin.de/polymake/ @* - in the optional argument # it is possible to hand over other names for the variables to be used -- be carful, the format must be correct and that is not tested, e.g. if you want the variable names to be u00,u10,u01,u11 then you must hand over the string u11,u10,u01,u11 @* - if the triangluations are not handed over as optional argument the procedure calls for its computation of these triangulations the program points2triangs from the program topcom by Joerg Rambau, Universitaet Bayreuth; it therefore is necessary that this program is installed in order to use this procedure; see @* http://www.uni-bayreuth.de/departments/wirtschaftsmathematik/rambau/TOPCOM EXAMPLE: example secondaryFan; shows an example" { if (size(#)==0) { list triang=triangulations(polygon); } else { list triang=#[1]; } list sp=secondaryPolytope(polygon,triang); list spp=polymakePolytope(sp[1]); list sf=normalFan(spp[1],spp[4],spp[3],1); return(list(sf[1],sf[2],spp,triang)); } example { "EXAMPLE:"; echo=2; // the lattice points of the unit square in the plane list polygon=intvec(0,0),intvec(0,1),intvec(1,0),intvec(1,1); // the secondary polytope of this lattice point configuration is computed list secfan=secondaryFan(polygon); // the number of cones in the secondary fan of the polygon size(secfan[1]); // the inequalities of the first cone as matrix are: print(secfan[1][1][1]); // the inequalities of the first cone as string are: print(secfan[1][1][2]); // the rows of the following matrix are the extreme rays of the first cone: print(secfan[1][1][3]); // each cone contains the linearity space spanned by: print(secfan[2]); // the points in the secondary polytope print(secfan[3][1]); // the corresponding triangulations secfan[4]; } //////////////////////////////////////////////////////////////////////////////// /// PROCEDURES CONCERNED WITH PLANAR POLYGONS //////////////////////////////////////////////////////////////////////////////// proc cycleLength (list boundary,intvec interior) "USAGE: cycleLength(boundary,interior); list boundary, intvec interior ASSUME: boundary is a list of integer vectors describing a cycle in some convex lattice polygon around the lattice point interior ordered clock wise RETURN: string, the cycle length of the corresponding cycle in the dual tropical curve EXAMPLE: example cycleLength; shows an example" { int j; // create a ring whose variables are indexed by the points in // boundary resp. by interior string rst="ring cyclering=0,(u"+string(interior[1])+string(interior[2]); for (j=1;j<=size(boundary);j++) { rst=rst+",u"+string(boundary[j][1])+string(boundary[j][2]); } rst=rst+"),lp;"; execute(rst); // add the first and second point at the end of boundary boundary[size(boundary)+1]=boundary[1]; boundary[size(boundary)+1]=boundary[2]; poly cl,summand; // takes the cycle length matrix N1[2][2]; // used to compute the area of a triangle matrix N2[2][2]; // used to compute the area of a triangle matrix N3[2][2]; // used to compute the area of a triangle // for each original point in boundary compute its contribution to the cycle for (j=2;j<=size(boundary)-1;j++) { N1=boundary[j-1]-interior,boundary[j]-interior; N2=boundary[j]-interior,boundary[j+1]-interior; N3=boundary[j+1]-interior,boundary[j-1]-interior; execute("summand=-u"+string(boundary[j][1])+string(boundary[j][2])+"+u"+string(interior[1])+string(interior[2])+";"); summand=summand*(det(N1)+det(N2)+det(N3))/(det(N1)*det(N2)); cl=cl+summand; } return(string(cl)); } example { "EXAMPLE:"; echo=2; // the integer vectors in boundary are lattice points on the boundary // of a convex lattice polygon in the plane list boundary=intvec(0,0),intvec(0,1),intvec(0,2),intvec(2,2), intvec(2,1),intvec(2,0); // interior is a lattice point in the interior of this lattice polygon intvec interior=1,1; // compute the general cycle length of a cycle of the corresponding cycle // in the dual tropical curve, note that (0,1) and (2,1) do not contribute cycleLength(boundary,interior); } ///////////////////////////////////////////////////////////////////////////// proc splitPolygon (list markings) "USAGE: splitPolygon (markings); markings list ASSUME: markings is a list of integer vectors representing lattice points in the plane which we consider as the marked points of the convex lattice polytope spanned by them PURPOSE: split the marked points in the vertices, the points on the facets which are not vertices, and the interior points RETURN: list, L consisting of three lists @* L[1] : represents the vertices the polygon ordered clockwise @* L[1][i][1] = intvec, the coordinates of the ith vertex @* L[1][i][2] = int, the position of L[1][i][1] in markings @* L[2][i] : represents the lattice points on the facet of the polygon with endpoints L[1][i] and L[1][i+1] (i considered modulo size(L[1])) @* L[2][i][j][1] = intvec, the coordinates of the jth lattice point on that facet @* L[2][i][j][2] = int, the position of L[2][i][j][1] in markings @* L[3] : represents the interior lattice points of the polygon @* L[3][i][1] = intvec, coordinates of ith interior point @* L[3][i][2] = int, the position of L[3][i][1] in markings EXAMPLE: example splitPolygon; shows an example" { list vert; // stores the result // compute the boundary of the polygon in an oriented way list pb=findOrientedBoundary(markings); // the vertices are just the second entry of pb vert[1]=pb[2]; int i,j,k; // indices list boundary; // stores the points on the facets of the // polygon which are not vertices // append to the boundary points as well as to the vertices // the first vertex a second time pb[1]=pb[1]+list(pb[1][1]); pb[2]=pb[2]+list(pb[2][1]); // for each vertex find all points on the facet of the polygon with this vertex // and the next vertex as endpoints int z=2; for (i=1;i<=size(vert[1]);i++) { j=1; list facet; // stores the points on this facet which are not vertices // while the next vertex is not reached, store the boundary lattice point while (pb[1][z]!=pb[2][i+1]) { facet[j]=pb[1][z]; j++; z++; } // store the points on the ith facet as boundary[i] boundary[i]=facet; kill facet; z++; } // store the information on the boundary in vert[2] vert[2]=boundary; // find the remaining points in the input which are not on // the boundary by checking // for each point in markings if it is contained in pb[1] list interior=markings; for (i=size(interior);i>=1;i--) { for (j=1;j<=size(pb[1])-1;j++) { if (interior[i]==pb[1][j]) { interior=delete(interior,i); j=size(pb[1]); } } } // store the interior points in vert[3] vert[3]=interior; // add to each point in vert the index which it gets from // its position in the input markings; // do so for ver[1] for (i=1;i<=size(vert[1]);i++) { j=1; while (markings[j]!=vert[1][i]) { j++; } vert[1][i]=list(vert[1][i],j); } // do so for ver[2] for (i=1;i<=size(vert[2]);i++) { for (k=1;k<=size(vert[2][i]);k++) { j=1; while (markings[j]!=vert[2][i][k]) { j++; } vert[2][i][k]=list(vert[2][i][k],j); } } // do so for ver[3] for (i=1;i<=size(vert[3]);i++) { j=1; while (markings[j]!=vert[3][i]) { j++; } vert[3][i]=list(vert[3][i],j); } return(vert); } example { "EXAMPLE:"; echo=2; // the lattice polygon spanned by the points (0,0), (3,0) and (0,3) // with all integer points as markings list polygon=intvec(1,1),intvec(3,0),intvec(2,0),intvec(1,0), intvec(0,0),intvec(2,1),intvec(0,1),intvec(1,2), intvec(0,2),intvec(0,3); // split the polygon in its vertices, its facets and its interior points list sp=splitPolygon(polygon); // the vertices sp[1]; // the points on facets which are not vertices sp[2]; // the interior points sp[3]; } ///////////////////////////////////////////////////////////////////////////// proc eta (list triang,list polygon) "USAGE: eta(triang,polygon); triang, polygon list ASSUME: polygon has the format of the output of splitPolygon, i.e. it is a list with three entries describing a convex lattice polygon in the following way: @* polygon[1] : is a list of lists; for each i the entry polygon[1][i][1] is a lattice point which is a vertex of the lattice polygon, and polygon[1][i][2] is an integer assigned to this lattice point as identifying index @* polygon[2] : is a list of lists; for each vertex of the polygon, i.e. for each entry in polygon[1], it contains a list polygon[2][i], which contains the lattice points on the facet with endpoints polygon[1][i] and polygon[1][i+1] - i considered mod size(polygon[1]); each such lattice point contributes an entry polygon[2][i][j][1] which is an integer vector giving the coordinate of the lattice point and an entry polygon[2][i][j][2] which is the identifying index @* polygon[3] : is a list of lists, where each entry corresponds to a lattice point in the interior of the polygon, with polygon[3][j][1] being the coordinates of the point and polygon[3][j][2] being the identifying index; @* triang is a list of integer vectors all of size three describing a triangulation of the polygon described by polygon; if an entry of triang is the vector (i,j,k) then the triangle is build by the vertices with indices i, j and k RETURN: intvec, the integer vector eta describing that vertex of the Newton polytope discriminant of the polygone whose dual cone in the Groebner fan contains the cone of the secondary fan of the polygon corresponding to the given triangulation NOTE: for a better description of eta see either Gelfand, Kapranov, Zelevinski: Discriminants, Resultants and multidimensional Determinants. Chapter 10. EXAMPLE: example eta; shows an example" { int i,j,k,l,m,n; // index variables list ordpolygon; // stores the lattice points in the order // used in the triangulation list triangarea; // stores the areas of the triangulations intmat N[2][2]; // used to compute triangle areas // 1) store the lattice points in the order used in the triangulation // go first through all vertices of the polytope for (j=1;j<=size(polygon[1]);j++) { ordpolygon[polygon[1][j][2]]=polygon[1][j][1]; } // then consider all inner points for (j=1;j<=size(polygon[3]);j++) { ordpolygon[polygon[3][j][2]]=polygon[3][j][1]; } // finally consider all lattice points on the boundary which are not vertices for (j=1;j<=size(polygon[2]);j++) { for (i=1;i<=size(polygon[2][j]);i++) { ordpolygon[polygon[2][j][i][2]]=polygon[2][j][i][1]; } } // 2) compute for each triangle in the triangulation the area of the triangle for (i=1;i<=size(triang);i++) { // Note that the ith lattice point in orderedpolygon has the // number i-1 in the triangulation! N=ordpolygon[triang[i][1]]-ordpolygon[triang[i][2]],ordpolygon[triang[i][1]]-ordpolygon[triang[i][3]]; triangarea[i]=abs(det(N)); } intvec ETA; // stores the eta_ij int etaij; // stores the part of eta_ij during computations // which comes from triangle areas int seitenlaenge; // stores the part of eta_ij during computations // which comes from boundary facets list seiten; // stores the lattice points on facets of the polygon intvec v; // used to compute a facet length // 3) store first in seiten[i] all lattice points on the facet // connecting the ith vertex, // i.e. polygon[1][i], with the i+1st vertex, i.e. polygon[1][i+1], // where we replace i+1 // 1 if i=size(polygon[1]); // then append the last entry of seiten once more at the very // beginning of seiten, so // that the index is shifted by one for (i=1;i<=size(polygon[1]);i++) { if (imp[1]) { abstand[i]=laenge[mp[1],i]; } } polygon=sortlistbyintvec(polygon,abstand); return(list(polygon,endpoints)); } /////////////////////////////////////////////////////////////// list orderedvertices; // stores the vertices in an ordered way list minimisedorderedvertices; // stores the vertices in an ordered way; // redundant ones removed list comparevertices; // stores vertices which should be compared to // the testvertex orderedvertices[1]=polygon[1]; // set the starting vertex minimisedorderedvertices[1]=polygon[1]; // set the starting vertex intvec testvertex=polygon[1]; //vertex to which the others have to be compared intvec startvertex=polygon[1]; // keep the starting vertex to test, // when the end is reached int endtest; // is set to one, when the end is reached int startvertexfound;// is 1, once for some testvertex a candidate // for the next vertex has been found polygon=delete(polygon,1); // delete the testvertex intvec v,w; int l=1; // counts the vertices // the basic idea is that a vertex can be // the next one on the boundary if all other vertices // ly to the right of the vector v pointing // from the testvertex to this one; this can be tested // by checking if the determinant of the 2x2-matrix // with first column v and second column the vector w, // pointing from the testvertex to the new vertex, // is non-positive; if this is the case for all // new vertices, then the one in consideration is // a possible choice for the next vertex on the boundary // and it is stored in naechste; we can then order // the candidates according to their distance from // the testvertex; then they occur on the boundary in that order! while (endtest==0) { list naechste; // stores the possible choices for the next vertex k=1; for (i=1;i<=size(polygon);i++) { d=0; // stores the value of the determinant of (v,w) v=polygon[i]-testvertex; // points from the testvertex to the ith vertex comparevertices=delete(polygon,i); // we needn't compare v to itself // we should compare v to the startvertex-testvertex; // in the first calling of the loop // this is irrelevant since the difference will be zero; // however, later on it will // be vital, since we delete the vertices // which we have already tested from the list // of all vertices, and when all vertices // on the boundary have been found we would // therefore find a vertex in the interior // as candidate; but always testing against // the starting vertex, this can not happen comparevertices[size(comparevertices)+1]=startvertex; for (j=1;(j<=size(comparevertices)) and (d<=0);j++) { w=comparevertices[j]-testvertex; // points form the testvertex // to the jth vertex D=v,w; d=det(D); } if (d<=0) // if all determinants are non-positive, { // then the ith vertex is a candidate naechste[k]=list(polygon[i],i,scalarproduct(v,v));// we store the vertex, //its position, and its k++; // distance from the testvertex } } if (size(naechste)>0) // then a candidate for the next vertex has been found { startvertexfound=1; // at least once a candidate has been found naechste=sortlist(naechste,3); // we order the candidates according // to their distance from testvertex; for (j=1;j<=size(naechste);j++) // then we store them in this { // order in orderedvertices l++; orderedvertices[l]=naechste[j][1]; } testvertex=naechste[size(naechste)][1]; // we store the last one as // next testvertex; // store the next corner of NSD minimisedorderedvertices[size(minimisedorderedvertices)+1]=testvertex; naechste=sortlist(naechste,2); // then we reorder the vertices // according to their position for (j=size(naechste);j>=1;j--) // and we delete them from the vertices { polygon=delete(polygon,naechste[j][2]); } } else // that means either that the vertex was inside the polygon, { // or that we have reached the last vertex on the boundary // of the polytope if (startvertexfound==0) // the vertex was in the interior; { // we delete it and start all over again orderedvertices[1]=polygon[1]; minimisedorderedvertices[1]=polygon[1]; testvertex=polygon[1]; startvertex=polygon[1]; polygon=delete(polygon,1); } else // we have reached the last vertex on the boundary of { // the polytope and can stop endtest=1; } } kill naechste; } // test if the first vertex in minimisedorderedvertices // is on the same line with the second and // the last, i.e. if we started our search in the // middle of a face; if so, delete it v=minimisedorderedvertices[2]-minimisedorderedvertices[1]; w=minimisedorderedvertices[size(minimisedorderedvertices)]-minimisedorderedvertices[1]; D=v,w; if (det(D)==0) { minimisedorderedvertices=delete(minimisedorderedvertices,1); } // test if the first vertex in minimisedorderedvertices // is on the same line with the two // last ones, i.e. if we started our search at the end of a face; // if so, delete it v=minimisedorderedvertices[size(minimisedorderedvertices)-1]-minimisedorderedvertices[1]; w=minimisedorderedvertices[size(minimisedorderedvertices)]-minimisedorderedvertices[1]; D=v,w; if (det(D)==0) { minimisedorderedvertices=delete(minimisedorderedvertices,size(minimisedorderedvertices)); } return(list(orderedvertices,minimisedorderedvertices)); } example { "EXAMPLE:"; echo=2; // the following lattice points in the plane define a polygon list polygon=intvec(0,0),intvec(3,1),intvec(1,0),intvec(2,0), intvec(1,1),intvec(3,2),intvec(1,2),intvec(2,3), intvec(2,4); // we compute its boundary list boundarypolygon=findOrientedBoundary(polygon); // the points on the boundary ordered clockwise are boundarypolygon[1] boundarypolygon[1]; // the vertices of the boundary are boundarypolygon[2] boundarypolygon[2]; } ///////////////////////////////////////////////////////////////////////////// proc cyclePoints (list triang,list points,int pt) "USAGE: cyclePoints(triang,points,pt) triang,points list, pt int ASSUME: - points is a list of integer vectors describing the lattice points of a marked polygon; @* - triang is a list of integer vectors describing a triangulation of the marked polygon in the sense that an integer vector of the form (i,j,k) describes the triangle formed by polygon[i], polygon[j] and polygon[k]; @* - pt is an integer between 1 and size(points), singling out a lattice point among the marked points PURPOSE: consider the convex lattice polygon, say P, spanned by all lattice points in points which in the triangulation triang are connected to the point points[pt]; the procedure computes all marked points in points which ly on the boundary of that polygon, ordered clockwise RETURN: list, of integer vectors which are the coordinates of the lattice points on the boundary of the above mentioned polygon P, if this polygon is not the empty set (that would be the case if points[pt] is not a vertex of any triangle in the triangulation); otherwise return the empty list EXAMPLE: example cyclePoints; shows an example" { int i,j; // indices list v; // saves the indices of lattice points connected to the // interior point in the triangulation // save all points in triangulations containing pt in v for (i=1;i<=size(triang);i++) { if ((triang[i][1]==pt) or (triang[i][2]==pt) or (triang[i][3]==pt)) { j++; v[3*j-2]=triang[i][1]; v[3*j-1]=triang[i][2]; v[3*j]=triang[i][3]; } } if (size(v)==0) { return(list()); } // remove pt itself and redundancies in v for (i=size(v);i>=1;i--) { j=1; while ((j16)) { ERROR("n is not between 1 and 16."); } if (size(#)>0) { if ((defined(x)==0) or (defined(y)==0)) { ERROR("The variables x and y are not defined."); } } if ((defined(x)==0) or (defined(y)==0)) { ring nfring=0,(x,y),dp; } // store the normal forms as polynomials list nf=x2+y+xy2,x2+x+1+xy2,x3+x2+x+1+y2+y,x4+x3+x2+x+1+y2+y+x2y,x3+x2+x+1+x2y+y+xy2+y2+y3, x2+x+x2y+y2,x2y+x+y+xy2,x2+x+1+y+xy2,x2+x+1+y+xy2+y2,x3+x2+x+1+x2y+y+y2,x3+x2+x+1+x2y+y+xy2+y2, x2+x+1+x2y+y+x2y2+xy2+y2,x+1+x2y+y+xy2,x2+x+1+x2y+y+xy2,x2+x+1+x2y+y+xy2+y2,x2+x+x2y+y+xy2+y2; list pg=newtonPolytopeLP(nf[n]); if (size(#)==0) { return(list(findOrientedBoundary(pg)[2],pg)); } else { return(list(findOrientedBoundary(pg)[2],pg,nf[n])); } } example { "EXAMPLE:"; echo=2; list nf=ellipticNFDB(5); // the vertices of the 5th normal form are nf[1]; // its lattice points are nf[2]; } ///////////////////////////////////////////////////////////////////////////////// /// AUXILARY PROCEDURES ///////////////////////////////////////////////////////////////////////////////// proc polymakeKeepTmpFiles (int i) "USAGE: polymakeKeepTmpFiles(int i); i int PURPOSE: some procedures create files in the directory /tmp which are used for computations with polymake respectively topcom; these will be removed when the corresponding procedure is left; however, it might be desireable to keep them for further computations with either polymake or topcom; this can be achieved by this procedure; call the procedure as: @* - polymakeKeepTmpFiles(1); - then the files will be kept @* - polymakeKeepTmpFiles(0); - then files will be removed in the future RETURN: none" { if ((i==1) and (defined(polymakekeeptmpfiles)==0)) { int polymakekeeptmpfiles; export(polymakekeeptmpfiles); } if (i!=1) { if (defined(polymakekeeptmpfiles)) { kill polymakekeeptmpfiles; } } } ///////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////// /// AUXILARY PROCEDURES, WHICH ARE DECLARED STATIC ///////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////// /// - scalarproduct /// - intmatcoldelete /// - intmatconcat /// - sortlist /// - minInList /// - stringdelete /// - abs /// - commondenominator /// - maxPosInIntvec /// - maxPosInIntmat /// - sortintvec /// - matrixtointmat ///////////////////////////////////////////////////////////////////////////////// static proc scalarproduct (intvec w,intvec v) "USAGE: scalarproduct(w,v); w,v intvec ASSUME: w and v are integer vectors of the same length RETURN: int, the scalarproduct of v and w NOTE: the procedure is called by findOrientedBoundary" { int sp; for (int i=1;i<=size(w);i++) { sp=sp+v[i]*w[i]; } return(sp); } static proc intmatcoldelete (intmat w,int i) "USAGE: intmatcoldelete(w,i); w intmat, i int RETURN: intmat, the integer matrix w with the ith comlumn deleted NOTE: the procedure is called by intmatsort and normalFan" { if ((i<1) or (i>ncols(w)) or (ncols(w)==1)) { return(w); } if (i==1) { intmat M[nrows(w)][ncols(w)-1]=w[1..nrows(w),2..ncols(w)]; return(M); } if (i==ncols(w)) { intmat M[nrows(w)][ncols(w)-1]=w[1..nrows(w),1..ncols(w)-1]; return(M); } else { intmat M[nrows(w)][i-1]=w[1..nrows(w),1..i-1]; intmat N[nrows(w)][ncols(w)-i]=w[1..nrows(w),i+1..ncols(w)]; return(intmatconcat(M,N)); } } static proc intmatconcat (intmat M,intmat N) "USAGE: intmatconcat(M,N); M,N intmat RETURN: intmat, M and N concatenated NOTE: the procedure is called by intmatcoldelete and sortintmat" { if (nrows(M)>=nrows(N)) { int m=nrows(M); } else { int m=nrows(N); } intmat P[m][ncols(M)+ncols(N)]; P[1..nrows(M),1..ncols(M)]=M[1..nrows(M),1..ncols(M)]; P[1..nrows(N),ncols(M)+1..ncols(M)+ncols(N)]=N[1..nrows(N),1..ncols(N)]; return(P); } static proc sortlist (list v,int pos) "USAGE: sortlist(v,pos); v list, pos int RETURN: list, the list L ordered in an ascending way according to the pos-th entries NOTE: called by tropicalCurve" { if(size(v)==1) { return(v); } list w=minInList(v,pos); v=delete(v,w[2]); v=sortlist(v,pos); v=list(w[1])+v; return(v); } static proc minInList (list v,int pos) "USAGE: minInList(v,pos); v list, pos int RETURN: list, (v[i],i) such that v[i][pos] is minimal NOTE: called by sortlist" { int min=v[1][pos]; int minpos=1; for (int i=2;i<=size(v);i++) { if (v[i][pos]size(w)) or (i<=0)) { return(w); } if ((size(w)==1) and (i==1)) { return(""); } if (i==1) { return(w[2..size(w)]); } if (i==size(w)) { return(w[1..size(w)-1]); } else { string erg=w[1..i-1],w[i+1..size(w)]; return(erg); } } static proc abs (def n) "USAGE: abs(n); n poly or int RETURN: poly or int, the absolute value of n" { if (n>=0) { return(n); } else { return(-n); } } static proc commondenominator (matrix M) "USAGE: commondenominator(M); M matrix ASSUME: the base ring has characteristic zero RETURN: int, the lowest common multiple of the denominators of the leading coefficients of the entries in M NOTE: the procedure is called from polymakeToIntmat" { int i,j; int kgV=1; // successively build the lowest common multiple of the denominators of the leading coefficients // of the entries in M for (i=1;i<=nrows(M);i++) { for (j=1;j<=ncols(M);j++) { kgV=lcm(kgV,int(denominator(leadcoef(M[i,j])))); } } return(kgV); } static proc maxPosInIntvec (intvec v) "USAGE: maxPosInIntvec(v); v intvec RETURN: int, the first position of a maximal entry in v NOTE: called by sortintmat" { int max=v[1]; int maxpos=1; for (int i=2;i<=size(v);i++) { if (v[i]>max) { max=v[i]; maxpos=i; } } return(maxpos); } static proc maxPosInIntmat (intmat v) "USAGE: maxPosInIntmat(v); v intmat ASSUME: v has a unique maximal entry RETURN: intvec, the position (i,j) of the maximal entry in v NOTE: called by findOrientedBoundary" { int max=v[1,1]; intvec maxpos=1,1; int i,j; for (i=1;i<=nrows(v);i++) { for (j=1;j<=ncols(v);j++) { if (v[i,j]>max) { max=v[i,j]; maxpos=i,j; } } } return(maxpos); } static proc sortintvec (intvec w) "USAGE: sortintvec(v); v intvec RETURN: intvec, the entries of v are ordered in an ascending way NOTE: called from ellipticNF" { int j,k,stop; intvec v=w[1]; for (j=2;j<=size(w);j++) { k=1; stop=0; while ((k<=size(v)) and (stop==0)) { if (v[k]