source: git/Singular/LIB/polymake.lib @ e84f07d

fieker-DuValspielwiese
Last change on this file since e84f07d was e84f07d, checked in by Thomas Markwig <keilen@…>, 15 years ago
einige Hilfstexte verbessert git-svn-id: file:///usr/local/Singular/svn/trunk@11561 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 88.4 KB
Line 
1version="$Id: polymake.lib,v 1.11 2009-03-16 11:34:10 keilen Exp $";
2category="Tropical Geometry";
3info="
4LIBRARY:        polymake.lib    Computations with polytopes and fans, interface to polymake and TOPCOM
5AUTHOR:         Thomas Markwig,  email: keilen@mathematik.uni-kl.de
6
7WARNING:
8     Most procedures will not work unless polymake or topcom is installed and
9     if so, they will only work with the operating system LINUX!
10     For more detailed information see IMPORTANT NOTE respectively consult the
11     help string of the procedures.
12
13IMPORTANT NOTE:
14     Even though this is a Singular library for computing polytopes and fans such
15     as the Newton polytope or the Groebner fan of a polynomial, most of the hard
16     computations are NOT done by Singular but by the program
17@*   -  polymake by Ewgenij Gawrilow, TU Berlin and Michael Joswig, TU Darmstadt
18@*      (see http://www.math.tu-berlin.de/polymake/),
19@*   respectively (only in the procedure triangularions) by the program
20@*   -  topcom by Joerg Rambau, Universitaet Bayreuth
21@*      (see http://www.uni-bayreuth.de/departments/wirtschaftsmathematik/rambau/TOPCOM);
22@*   this library should rather be seen as an interface which allows to use a (very limited)
23     number of options which polymake respectively topcom offers to compute with polytopes
24     and fans and to make the results available in Singular for further computations;
25     moreover, the user familiar with Singular does not have to learn the syntax of polymake
26     or topcom, if the options offered here are sufficient for his purposes.
27@*   Note, though, that the procedures concerned with planar polygons are independent of
28     both, polymake and topcom.
29
30PROCEDURES USING POLYMAKE:
31     polymakePolytope(list)            computes the vertices of a polytope using polymake
32     newtonPolytope(poly)              computes the Newton polytope of the polynomial
33     newtonPolytopeLP(poly)            computes the lattice points of the Newton polytope of the polynomial
34     normalFan(intmat,intmat,list)     computes the normal fan of a polytope
35     groebnerFan(poly)                 computes the Groebner fan of a polynomial
36     intmatToPolymake(intmat,string)   transforms an integer matrix into polymake format
37     polymakeToIntmat(string,string)   transforms polymake output into an integer matrix
38
39PROCEDURES USING TOPCOM:
40     triangulations(list)          computes all triangulations of a marked polytope
41     secondaryPolytope(list)       computes the secondary polytope of a marked polytope
42
43PROCEDURES USING POLYMAKE AND TOPCOM:
44     secondaryFan(list)            computes the secondary fan of a marked polytope
45
46PROCEDURES CONERNED WITH PLANAR POLYGONS:
47     cycleLength(list,intvec)      computes the cycleLength of cycle dual to list with interior point intvec
48     splitPolygon(list)            splits a marked polygon into vertices, facets and interior points
49     eta(list,list)                computes the eta-vector of a triangulation
50     findOrientedBoundary(list)    computes the boundary of the convex hull of a list of lattice points
51     cyclePoints(list,list,int)    computes lattice points connected to a lattice point in a triangulation
52     latticeArea(list)             computes the lattice area of a polygon
53     picksFormula(list)            computes the ingrediants of Pick's formula for a polygon
54     ellipticNF(list)              computes the normal form of an elliptic polygon
55     ellipticNFDB(int)             displays the 16 normal forms of elliptic polygons
56
57AUXILARY PROCEDURES:
58     polymakeKeepTmpFiles(int)     determines whether the files created in /tmp should be kept
59
60KEYWORDS:    polytope; fan; secondary fan; secondary polytope; polymake;
61             Newton polytope; Groebner fan
62
63";
64
65////////////////////////////////////////////////////////////////////////////////////
66/// Auxilary Static Procedures in this Library
67////////////////////////////////////////////////////////////////////////////////////
68/// - scalarproduct
69/// - intmatcoldelete
70/// - intmatconcat
71/// - sortlist
72/// - minInList
73/// - stringdelete
74/// - abs
75/// - commondenominator
76/// - maxPosInIntvec
77/// - maxPosInIntmat
78/// - sortintvec
79/// - matrixtointmat
80/////////////////////////////////////////////////////////////////////////////////////
81
82////////////////////////////////////////////////////////////////////////////////////
83LIB "poly.lib";
84LIB "linalg.lib";
85LIB "random.lib";
86////////////////////////////////////////////////////////////////////////////////////
87
88
89////////////////////////////////////////////////////////////////////////////////////////
90/// PROCEDURES USING POLYMAKE
91////////////////////////////////////////////////////////////////////////////////////////
92
93proc polymakePolytope (intmat polytope,list #)
94"USAGE:  polymakePolytope(polytope[,#]);   polytope list, # string
95ASSUME:  each row of polytope gives the coordinates of a lattice point of a polytope
96         with their affine coordinates as given by the output of secondaryPolytope
97PURPOSE: the procedure calls polymake to compute the vertices of the polytope as well
98         as its dimension and information on its facets
99RETURN:  list, L with four entries
100@*            L[1] : an integer matrix whose rows are the coordinates of vertices
101                     of the polytope
102@*            L[2] : the dimension of the polytope
103@*            L[3] : a list whose ith entry explains to which vertices the ith vertex
104                     of the Newton polytope is connected -- i.e. L[3][i] is an integer
105                     vector and an entry k in there means that the vertex L[1][i] is
106                     connected to the vertex L[1][k]
107@*            L[4] : an integer matrix whose rows mulitplied by (1,var(1),...,var(nvar)) give
108                     a linear system of equations describing the affine hull of the polytope,
109                     i.e. the smallest affine space containing the polytope
110NOTE: -  for its computations the procedure calls the program polymake by Ewgenij Gawrilow,
111         TU Berlin and Michael Joswig, TU Darmstadt; it therefore is necessary that
112         this program is installed in order to use this procedure;
113         see http://www.math.tu-berlin.de/polymake/
114@*    -  note that in the vertex edge graph we have changed the polymake convention which
115         starts indexing its vertices by zero while we start with one !
116@*    -  the procedure creates the file  /tmp/polytope.polymake which contains the polytope
117         in polymake format; if you wish to use this for further computations with polymake,
118         you have to use the procedure polymakeKeepTmpFiles in before
119@*    -  moreover, the procedure creates the file /tmp/polytope.output which it deletes
120         again before ending
121@*    -  it is possible to give as an optional second argument as string which then will be
122         used instead of 'polytope' in the name of the polymake output file
123EXAMPLE: example polymakePolytope;   shows an example"
124{
125  // the header for the file secendarypolytope.polymake
126  string sp="_application polytope
127_version 2.2
128_type RationalPolytope
129
130POINTS
131";
132  int i,j;
133  // set the name for the polymake output file
134  if (size(#)>0)
135  {
136    if (typeof(#[1])=="string")
137    {
138      string dateiname=#[1];
139    }
140    else
141    {
142      string dateiname="polytope";
143    }
144  }
145  else
146  {
147    string dateiname="polytope";
148  }
149  // create the lattice point list for polymake
150  sp=sp+intmatToPolymake(polytope,"points");
151  // initialise dateiname.polymake and compute the vertices
152  write(":w /tmp/"+dateiname+".polymake",sp);
153  system("sh","cd /tmp; polymake "+dateiname+".polymake VERTICES > "+dateiname+".output");
154  string vertices=read("/tmp/"+dateiname+".output");
155  system("sh","/bin/rm /tmp/"+dateiname+".output");
156  intmat np=polymakeToIntmat(vertices,"affine");
157  // compute the dimension
158  system("sh","cd /tmp; polymake "+dateiname+".polymake DIM > "+dateiname+".output");
159  string pdim=read("/tmp/"+dateiname+".output");
160  system("sh","/bin/rm /tmp/"+dateiname+".output");
161  pdim=pdim[5,size(pdim)-6];
162  execute("int nd="+pdim+";");
163  // compute the vertex-edge graph
164  system("sh","cd /tmp; polymake "+dateiname+".polymake GRAPH > "+dateiname+".output");
165  string vertexedgegraph=read("/tmp/"+dateiname+".output");
166  system("sh","/bin/rm /tmp/"+dateiname+".output");
167  vertexedgegraph=vertexedgegraph[7,size(vertexedgegraph)-8];
168  string newveg;
169  for (i=1;i<=size(vertexedgegraph);i++)
170  {
171    if (vertexedgegraph[i]=="{")
172    {
173      newveg=newveg+"intvec(";
174    }
175    else
176    {
177      if (vertexedgegraph[i]=="}")
178      {
179        newveg=newveg+"),";
180      }
181      else
182      {
183        if (vertexedgegraph[i]==" ")
184        {
185          newveg=newveg+",";
186        }
187        else
188        {
189          newveg=newveg+vertexedgegraph[i];
190        }
191      }
192    }
193  }
194  newveg=newveg[1,size(newveg)-1];
195  execute("list nveg="+newveg+";");
196  // raise each entry in nveg by one
197  for (i=1;i<=size(nveg);i++)
198  {
199    for (j=1;j<=size(nveg[i]);j++)
200    {
201      nveg[i][j]=nveg[i][j]+1;
202    }
203  }
204  // compute the affine hull
205  system("sh","cd /tmp; polymake "+dateiname+".polymake AFFINE_HULL > "+dateiname+".output");
206  string equations=read("/tmp/"+dateiname+".output");
207  system("sh","/bin/rm /tmp/"+dateiname+".output");
208  if (size(equations)>14)
209  {
210    intmat neq=polymakeToIntmat(equations,"cleardenom");
211  }
212  else
213  {
214    intmat neq[1][ncols(polytope)+1];
215  }
216  // delete the tmp-files, if polymakekeeptmpfiles is not set
217  if (defined(polymakekeeptmpfiles)==0)
218  {
219    system("sh","/bin/rm /tmp/"+dateiname+".polymake");
220  }
221  // return the files
222  return(list(np,nd,nveg,neq));
223}
224example
225{
226   "EXAMPLE:";
227   echo=2;
228   // the lattice points of the unit square in the plane
229   list points=intvec(0,0),intvec(0,1),intvec(1,0),intvec(1,1);
230   // the secondary polytope of this lattice point configuration is computed
231   intmat secpoly=secondaryPolytope(points)[1];
232   list np=polymakePolytope(secpoly);
233   // the vertices of the secondary polytope are:
234   np[1];
235   // its dimension is
236   np[2];
237   // np[3] contains information how the vertices are connected to each other,
238   // e.g. the first vertex (number 0) is connected to the second one
239   np[3][1];
240   // the affine hull has the equation
241   ring r=0,x(1..4),dp;
242   matrix M[5][1]=1,x(1),x(2),x(3),x(4);
243   np[4]*M;
244}
245
246proc newtonPolytope (poly f,list #)
247"USAGE: newtonPolytope(f[,#]);  f poly, # string
248RETURN: list, L with four entries
249@*            L[1] : an integer matrix whose rows are the coordinates of vertices
250                     of the Newton polytope of f
251@*            L[2] : the dimension of the Newton polytope of f
252@*            L[3] : a list whose ith entry explains to which vertices the ith vertex
253                     of the Newton polytope is connected -- i.e. L[3][i] is an integer
254                     vector and an entry k in there means that the vertex L[1][i] is
255                     connected to the vertex L[1][k]
256@*            L[4] : an integer matrix whose rows mulitplied by (1,var(1),...,var(nvar)) give
257                     a linear system of equations describing the affine hull of the Newton
258                     polytope, i.e. the smallest affine space containing the Newton polytope
259NOTE: -  if we replace the first column of L[4] by zeros, i.e. if we move the affine hull to
260         the origin, then we get the equations for the orthogonal comploment of the linearity
261         space of the normal fan dual to the Newton polytope, i.e. we get the EQUATIONS that
262         we need as input for polymake when computing the normal fan
263@*    -  the procedure calls for its computation polymake by Ewgenij Gawrilow,
264         TU Berlin and Michael Joswig, so it only works if polymake is installed;
265         see http://www.math.tu-berlin.de/polymake/
266@*    -  the procedure creates the file  /tmp/newtonPolytope.polymake which contains the polytope
267         in polymake format and which can be used for further computations with polymake
268@*    -  moreover, the procedure creates the file /tmp/newtonPolytope.output which it deletes
269         again before ending
270@*    -  it is possible to give as an optional second argument as string which then will be
271         used instead of 'newtonPolytope' in the name of the polymake output file
272EXAMPLE: example newtonPolytope;   shows an example"
273{
274  int i,j;
275  // compute the list of exponent vectors of the polynomial, which are the lattice points
276  // whose convex hull is the Newton polytope of f
277  intmat exponents[size(f)][nvars(basering)];
278  while (f!=0)
279  {
280    i++;
281    exponents[i,1..nvars(basering)]=leadexp(f);
282    f=f-lead(f);
283  }
284  if (size(#)==0)
285  {
286    #[1]="newtonPolytope";
287  }
288  // call polymakePolytope with exponents
289  return(polymakePolytope(exponents,#));
290}
291example
292{
293   "EXAMPLE:";
294   echo=2;
295   ring r=0,(x,y,z),dp;
296   matrix M[4][1]=1,x,y,z;
297   poly f=y3+x2+xy+2xz+yz+z2+1;
298   // the Newton polytope of f is
299   list np=newtonPolytope(f);
300   // the vertices of the Newton polytope are:
301   np[1];
302   // its dimension is
303   np[2];
304   // np[3] contains information how the vertices are connected to each other,
305   // e.g. the first vertex (number 0) is connected to the second, third and
306   //      fourth vertex
307   np[3][1];
308   //////////////////////////
309   f=x2-y3;
310   // the Newton polytope of f is
311   np=newtonPolytope(f);
312   // the vertices of the Newton polytope are:
313   np[1];
314   // its dimension is
315   np[2];
316   // the Newton polytope is contained in the affine space given
317   //     by the equations
318   np[4]*M;
319}
320
321proc newtonPolytopeLP (poly f)
322"USAGE:  newtonPolytopeLP(f);  f poly
323RETURN: list, the exponent vectors of the monomials occuring in f, i.e. the lattice
324              points of the Newton polytope of f
325EXAMPLE: example normalFan;   shows an example"
326{
327  list np;
328  int i=1;
329  while (f!=0)
330  {
331    np[i]=leadexp(f);
332    f=f-lead(f);
333    i++;
334  }
335  return(np);
336}
337example
338{
339   "EXAMPLE:";
340   echo=2;
341   ring r=0,(x,y,z),dp;
342   poly f=y3+x2+xy+2xz+yz+z2+1;
343   // the lattice points of the Newton polytope of f are
344   newtonPolytopeLP(f);
345}
346
347proc normalFan (intmat vertices,intmat affinehull,list graph,int er,list #)
348"USAGE:  normalFan (vert,aff,graph,rays,[,#]);   vert,aff intmat,  graph list, rays int, # string
349ASSUME:  - vert is an integer matrix whose rows are the coordinate of the vertices of
350           a convex lattice polytope;
351@*       - aff describes the affine hull of this polytope, i.e.
352           the smallest affine space containing it, in the following sense:
353           denote by n the number of columns of vert, then multiply aff by (1,x(1),...,x(n))
354           and set the resulting terms to zero in order to get the equations for the affine hull;
355@*       - the ith entry of graph is an integer vector describing to which vertices
356           the ith vertex is connected, i.e. a k as entry means that the vertex vert[i] is
357           connected to vert[k];
358@*       - the integer rays is either one (if the extreme rays should be computed) or zero (otherwise)
359RETURN:  list, the ith entry of L[1] contains information about the cone in the normal fan
360               dual to the ith vertex of the polytope
361@*             L[1][i][1] = integer matrix representing the inequalities which describe the
362                            cone dual to the ith vertex
363@*             L[1][i][2] = a list which contains the inequalities represented by L[i][1]
364                            as a list of strings, where we use the variables x(1),...,x(n)
365@*             L[1][i][3] = only present if 'er' is set to 1; in that case it is an interger matrix
366                            whose rows are the extreme rays of the cone
367@*             L[2] = is an integer matrix whose rows span the linearity space of the fan,
368                      i.e. the linear space which is contained in each cone
369NOTE:    - the procedure calls for its computation polymake by Ewgenij Gawrilow,
370           TU Berlin and Michael Joswig, so it only works if polymake is installed;
371           see http://www.math.tu-berlin.de/polymake/
372@*       - in the optional argument # it is possible to hand over other names for the
373           variables to be used -- be carful, the format must be correct and that is
374           not tested, e.g. if you want the variable names to be u00,u10,u01,u11
375           then you must hand over the string u11,u10,u01,u11
376EXAMPLE: example normalFan;   shows an example"
377{
378  list ineq; // stores the inequalities of the cones
379  int i,j,k;
380  // we work over the following ring
381  execute("ring ineqring=0,x(1.."+string(ncols(vertices))+"),lp;");
382  string greatersign=">";
383  // create the variable names
384  if (size(#)>0)
385  {
386    if (typeof(#[1])=="string")
387    {
388      kill ineqring;
389      execute("ring ineqring=0,("+#[1]+"),lp;");
390    }
391    if (size(#)>1)
392    {
393      greatersign="<";
394    }
395  }
396  //////////////////////////////////////////////////////////////////
397  // Compute first the inequalities of the cones
398  //////////////////////////////////////////////////////////////////
399  matrix VAR[1][ncols(vertices)]=maxideal(1);
400  matrix EXP[ncols(vertices)][1];
401  poly p,pl,pr;
402  // consider all vertices of the polytope
403  for (i=1;i<=nrows(vertices);i++)
404  {
405    // first we produce for each vertex in the polytope
406    // the inequalities describing the dual cone in the normal fan
407    list pp;  // contain strings representing the inequalities describing the normal cone
408    intmat ie[size(graph[i])][ncols(vertices)]; // contains the inequalities as rows
409    // consider all the vertices to which the ith vertex in the polytope is connected by an edge
410    for (j=1;j<=size(graph[i]);j++)
411    {
412      // produce the vector ie_j pointing from the jth vertex to the ith vertex;
413      // this will be the jth inequality for the cone in the normal fan dual to the ith vertex
414      ie[j,1..ncols(vertices)]=vertices[i,1..ncols(vertices)]-vertices[graph[i][j],1..ncols(vertices)];
415      EXP=ie[j,1..ncols(vertices)];
416      // build a linear polynomial with the entries of ie_j as coefficients
417      p=(VAR*EXP)[1,1];
418      pl,pr=0,0;
419      // separate the terms with positive coefficients in p from those with negative coefficients
420      for (k=1;k<=size(p);k++)
421      {
422        if (leadcoef(p[k])<0)
423        {
424          pr=pr-p[k];
425        }
426        else
427        {
428          pl=pl+p[k];
429        }
430      }
431      // build the string which represents the jth inequality for the cone dual to the ith vertex
432      // as polynomial inequality of type string, and store this in the list pp as jth entry
433      pp[j]=string(pl)+" "+greatersign+" "+string(pr);
434    }
435    // all inequalities for the ith vertex are stored in the list ineq
436    ineq[i]=list(ie,pp);
437    kill ie,pp; // kill certain lists
438  }
439  // remove the first column of affine hull to compute the linearity space
440  intmat linearity=intmatcoldelete(affinehull,1);
441  //////////////////////////////////////////////////////////////////
442  // Compute next the extreme rays of the cones
443  //////////////////////////////////////////////////////////////////
444  if (er==1)
445  {
446    list extremerays; // keeps the result
447    string polymake; // keeps polymake output
448    // the header for ineq.polymake
449    string head="_application polytope
450_version 2.2
451_type RationalPolytope
452
453INEQUALITIES
454";
455    // the tail for both polymake files
456    string tail="
457EQUATIONS
458";
459    tail=tail+intmatToPolymake(linearity,"rays");
460    string ungleichungen; // keeps the inequalities for the polymake code
461    intmat M; // the matrix keeping the inequalities
462    // create the file ineq.output
463    write(":w /tmp/ineq.output","");
464    int dimension; // keeps the dimension of the intersection the bad cones with the u11tobeseencone
465    for (i=1;i<=size(ineq);i++)
466    {
467      i,". Cone of ",nrows(vertices); // indicate how many vertices have been dealt with
468      ungleichungen=intmatToPolymake(ineq[i][1],"rays");
469      // write the inequalities to ineq.polymake and call polymake
470      write(":w /tmp/ineq.polymake",head+ungleichungen+tail);
471      ungleichungen=""; // clear ungleichungen
472      system("sh","cd /tmp; /bin/rm ineq.output; polymake ineq.polymake VERTICES > ineq.output");
473      // read the result of polymake
474      polymake=read("/tmp/ineq.output");
475      intmat VERT=polymakeToIntmat(polymake,"affine");
476      extremerays[i]=VERT;
477      kill VERT;
478    }
479    for (i=1;i<=size(ineq);i++)
480    {
481      ineq[i]=ineq[i]+list(extremerays[i]);
482    }
483  }
484  // delete the tmp-files, if polymakekeeptmpfiles is not set
485  if (defined(polymakekeeptmpfiles)==0)
486  {
487    system("sh","/bin/rm /tmp/ineq.polymake");
488    system("sh","/bin/rm /tmp/ineq.output");
489  }
490  // get the linearity space
491  return(list(ineq,linearity));
492}
493example
494{
495   "EXAMPLE:";
496   echo=2;
497   ring r=0,(x,y,z),dp;
498   matrix M[4][1]=1,x,y,z;
499   poly f=y3+x2+xy+2xz+yz+z2+1;
500   // the Newton polytope of f is
501   list np=newtonPolytope(f);
502   // the Groebner fan of f, i.e. the normal fan of the Newton polytope
503   list gf=normalFan(np[1],np[4],np[3],1,"x,y,z");
504   // the number of cones in the Groebner fan of f is:
505   size(gf[1]);
506   // the inequalities of the first cone as matrix are:
507   print(gf[1][1][1]);
508   // the inequalities of the first cone as string are:
509   print(gf[1][1][2]);
510   // the rows of the following matrix are the extreme rays of the first cone:
511   print(gf[1][1][3]);
512   // each cone contains the linearity space spanned by:
513   print(gf[2]);
514}
515
516proc groebnerFan (poly f,list #)
517"USAGE:  groebnerFan(f[,#]);  f poly, # string
518RETURN:  list, the ith entry of L[1] contains information about the ith cone in the Groebner fan
519               dual to the ith vertex in the Newton polytope of the f
520@*             L[1][i][1] = integer matrix representing the inequalities which describe the cone
521@*             L[1][i][2] = a list which contains the inequalities represented by L[1][i][1]
522                            as a list of strings
523@*             L[1][i][3] = an interger matrix whose rows are the extreme rays of the cone
524@*             L[2] = is an integer matrix whose rows span the linearity space of the fan,
525                      i.e. the linear space which is contained in each cone
526@*             L[3] = the Newton polytope of f in the format of the procedure newtonPolytope
527@*             L[4] = integer matrix where each row represents the exponet vector of one monomial
528                      occuring in the input polynomial
529NOTE:    - if you have alread computed the Newton polytope of f then you might want
530           to use the procedure normalFan instead in order to avoid doing costly computation
531           twice
532@*       - the procedure calls for its computation polymake by Ewgenij Gawrilow,
533           TU Berlin and Michael Joswig, so it only works if polymake is installed;
534           see http://www.math.tu-berlin.de/polymake/
535@*       - the procedure creates the file  /tmp/newtonPolytope.polymake which contains the
536           Newton polytope of f in polymake format and which can be used for further
537           computations with polymake
538@*       - it is possible to give as an optional second argument as string which then will be
539           used instead of 'newtonPolytope' in the name of the polymake output file
540EXAMPLE: example groebnerFan;   shows an example"
541{
542  int i,j;
543  // compute the list of exponent vectors of the polynomial, which are the lattice points
544  // whose convex hull is the Newton polytope of f
545  intmat exponents[size(f)][nvars(basering)];
546  while (f!=0)
547  {
548    i++;
549    exponents[i,1..nvars(basering)]=leadexp(f);
550    f=f-lead(f);
551  }
552  if (size(#)==0)
553  {
554    #[1]="newtonPolytope";
555  }
556  // call polymakePolytope with exponents
557  list newtonp=polymakePolytope(exponents,"newtonPolytope");
558  // get the variables as string
559  string variablen;
560  for (i=1;i<=nvars(basering);i++)
561  {
562    variablen=variablen+string(var(i))+",";
563  }
564  variablen=variablen[1,size(variablen)-1];
565  // call normalFan in order to compute the Groebner fan
566  list gf=normalFan(newtonp[1],newtonp[4],newtonp[3],1,variablen);
567  // append newtonp to gf
568  gf[3]=newtonp;
569  // append the exponent vectors to gf
570  gf[4]=exponents;
571  return(gf);
572}
573example
574{
575   "EXAMPLE:";
576   echo=2;
577   ring r=0,(x,y,z),dp;
578   matrix M[4][1]=1,x,y,z;
579   poly f=y3+x2+xy+2xz+yz+z2+1;
580   // the Newton polytope of f is
581   list gf=groebnerFan(f);
582   // the exponent vectors of f are ordered as follows
583   gf[4];
584   // the first cone of the groebner fan has the inequalities
585   gf[1][1][1];
586   // as a string they look like
587   gf[1][1][2];
588   // and it has the extreme rays
589   print(gf[1][1][3]);
590   // the linearity space is spanned by
591   print(gf[2]);
592   // the vertices of the Newton polytope are:
593   gf[3][1];
594   // its dimension is
595   gf[3][2];
596   // np[3] contains information how the vertices are connected to each other,
597   // e.g. the 1st vertex is connected to the 2nd, 3rd and 4th vertex
598   gf[3][3][1];
599}
600
601
602
603proc intmatToPolymake (intmat M,string art)
604"USAGE:  intmatToPolymake(M,art);  M intmat, art string
605ASSUME:  - M is an integer matrix which should be transformed into polymake format;
606@*       - art is one of the following strings:
607@*         + 'rays'   : indicating that a first column of 0's should be added
608@*         + 'points' : indicating that a first column of 1's should be added
609RETURN:  string, the matrix is transformed in a string and a first column has been added
610EXAMPLE: example intmatToPolymake;   shows an example"
611{
612  if (art=="rays")
613  {
614    string anf="0 ";
615  }
616  else
617  {
618    string anf="1 ";
619  }
620  string sp;
621  int i,j;
622  // create the lattice point list for polymake
623  for (i=1;i<=nrows(M);i++)
624  {
625    sp=sp+anf;
626    for (j=1;j<=ncols(M);j++)
627    {
628      sp=sp+string(M[i,j])+" ";
629      if (j==ncols(M))
630      {
631        sp=sp+"
632";
633      }
634    }
635  }
636  return(sp);
637}
638example
639{
640   "EXAMPLE:";
641   echo=2;
642   intmat M[3][4]=1,2,3,4,5,6,7,8,9,10,11,12;
643   intmatToPolymake(M,"rays");
644   intmatToPolymake(M,"points");
645}
646
647proc polymakeToIntmat (string pm,string art)
648"USAGE:  polymakeToIntmat(pm,art);  pm, art string
649ASSUME:  pm is the result of calling polymake with one 'argument' like VERTICES, AFFINE_HULL, etc.,
650         so that the first row of the string is the name of the corresponding 'argument' and
651         the further rows contain the result which consists of vectors either over the integers
652         or over the rationals
653RETURN:  intmat, the rows of the matrix are basically the vectors in pm from the second row on
654                 where each row has been multiplied with the lowest common multiple of the
655                 denominators of its entries so as to be an integer matrix; moreover,
656                 if art=='affine', then the first column is omitted since we only want affine
657                 coordinates
658EXAMPLE: example polymakeToIntmat;   shows an example"
659{
660  // we need a line break
661  string zeilenumbruch="
662";
663  // remove the 'argment' name, i.e. the first row of pm
664  while (pm[1]!=zeilenumbruch)
665  {
666    pm=stringdelete(pm,1);
667  }
668  pm=stringdelete(pm,1);
669  // find out how many entries each vector has, namely one more than 'spaces' in a row
670  int i=1;
671  int s=1;
672  int z=1;
673  while (pm[i]!=zeilenumbruch)
674  {
675    if (pm[i]==" ")
676    {
677      s++;
678    }
679    i++;
680  }
681  // if we want to have affine coordinates
682  if (art=="affine")
683  {
684    s--; // then there is one column less
685    // and the entry of the first column (in the first row) has to be removed
686    while (pm[1]!=" ")
687    {
688      pm=stringdelete(pm,1);
689    }
690    pm=stringdelete(pm,1);
691  }
692  // we add two line breaks at the end in order to have this as a stopping criterion
693  pm=pm+zeilenumbruch+zeilenumbruch;
694  // we now have to work through each row
695  for (i=1;i<=size(pm);i++)
696  {
697    // if there are two consecutive line breaks we are done
698    if ((pm[i]==zeilenumbruch) and (pm[i+1]==zeilenumbruch))
699    {
700      i=size(pm)+1;
701    }
702    else
703    {
704      // a line break has to be replaced by a comma
705      if (pm[i]==zeilenumbruch)
706      {
707        z++;
708        pm[i]=",";
709        // if we want to have affine coordinates, then we have to delete the first entry in each row
710        if (art=="affine")
711        {
712         while (pm[i+1]!=" ")
713         {
714           pm=stringdelete(pm,i+1);
715         }
716         pm=stringdelete(pm,i+1);
717        }
718      }
719      // a space has to be replaced by a comma
720      if (pm[i]==" ")
721      {
722        pm[i]=",";
723      }
724    }
725  }
726  // if we have introduced superflous commata at the end, they should be removed
727  while (pm[size(pm)]==",")
728  {
729    pm=stringdelete(pm,size(pm));
730  }
731  // since the matrix could be over the rationals, we need a ring with rational coefficients
732  ring zwischering=0,x,lp;
733  // create the matrix with the elements of pm as entries
734  execute("matrix mm["+string(z)+"]["+string(s)+"]="+pm+";");
735  // transform this into an integer matrix
736  matrix M[1][ncols(mm)]; // takes a row of mm
737  int cm;  // takes a lowest common multiple
738  // multiply each row by an integer such that its entries are integers
739  for (int j=1;j<=nrows(mm);j++)
740  {
741    M=mm[j,1..ncols(mm)];
742    cm=commondenominator(M);
743    for (i=1;i<=ncols(mm);i++)
744    {
745      mm[j,i]=cm*mm[j,i];
746    }
747  }
748  // transform the matrix mm into an integer matrix
749  execute("intmat im["+string(z)+"]["+string(s)+"]="+string(mm)+";");
750  return(im);
751}
752example
753{
754   "EXAMPLE:";
755   echo=2;
756   // this is the usual output of some polymake computation
757   string pm="VERTICES
7580 1 3 5/3 1/3 -1 -23/3 -1/3 5/3 1/3 1
7590 1 3 -23/3 5/3 1 5/3 1/3 1/3 -1/3 -1
7600 1 1 1/3 -1/3 -1 5/3 1/3 -23/3 5/3 3
7610 1 1 5/3 -23/3 3 1/3 5/3 -1/3 1/3 -1
7620 1 -1 1/3 5/3 3 -1/3 -23/3 1/3 5/3 1
7630 1 -1 -1/3 1/3 1 1/3 5/3 5/3 -23/3 3
7640 1 -1 1 3 -5 -1 3 -1 1 -1
7650 1 -1 -1 -1 -1 1 1 3 3 -5
7660 1 -5 3 1 -1 3 -1 1 -1 -1
767
768";
769   intmat PM=polymakeToIntmat(pm,"affine");
770   // note that the first column has been removed, since we asked for
771   // affine coordinates, and the denominators have been cleared
772   print(PM);
773}
774
775////////////////////////////////////////////////////////////////////////////////////////
776/// PROCEDURES USING TOPCOM
777////////////////////////////////////////////////////////////////////////////////////////
778
779proc triangulations (list polygon)
780"USAGE:  triangulations(polygon); list polygon
781ASSUME:  polygon is a list of integer vectors of the same size representing the affine
782         coordinates of the lattice points
783PURPOSE: the procedure considers the marked polytope given as the convex hull of
784         the lattice points and with these lattice points as markings; it then
785         computes all possible triangulations of this marked polytope
786RETURN:  list, each entry corresponds to one triangulation and the ith entry is
787               itself a list of integer vectors of size three, where each integer
788               vector defines one triangle in the triangulation by telling which
789               points of the input are the vertices of the triangle
790NOTE:  - the procedure calls for its computations the program points2triangs
791         from the program topcom by Joerg Rambau, Universitaet Bayreuth; it
792         therefore is necessary that this program is installed in order to use this
793         procedure; see
794@*       http://www.uni-bayreuth.de/departments/wirtschaftsmathematik/rambau/TOPCOM
795@*     - the procedure creates the files /tmp/triangulationsinput and /tmp/triangulationsoutput;
796         the former is used as input for points2triangs and the latter is its output
797         containing the triangulations of corresponding to points in the format
798         of points2triangs; if you wish to use this for further computations with topcom,
799         you have to use the procedure polymakeKeepTmpFiles in before
800@*     - note that an integer i in an integer vector representing a triangle refers to
801         the ith lattice point, i.e. polygon[i]; this convention is different from
802         TOPCOM's convention, where i would refer to the i-1st lattice point
803EXAMPLE: example triangulations;   shows an example"
804{
805  int i,j;
806  // prepare the input for points2triangs by writing the input polygon in the
807  // necessary format
808  string spi="[";
809  for (i=1;i<=size(polygon);i++)
810  {
811    polygon[i][size(polygon[i])+1]=1;
812    spi=spi+"["+string(polygon[i])+"]";
813    if (i<size(polygon))
814    {
815      spi=spi+",";
816    }
817  }
818  spi=spi+"]";
819  write(":w /tmp/triangulationsinput",spi);
820  // call points2triangs
821  system("sh","cd /tmp; points2triangs < triangulationsinput > triangulationsoutput");
822  string p2t=read("/tmp/triangulationsoutput"); // takes the result of points2triangs
823  // delete the tmp-files, if polymakekeeptmpfiles is not set
824  if (defined(polymakekeeptmpfiles)==0)
825  {
826    system("sh","cd /tmp; rm -f triangulationsinput; rm -f triangulationsoutput");
827  }
828  // preprocessing of p2t if points2triangs is version >= 0.15 brings p2t to the format of version 0.14
829  string np2t; // takes the triangulations in Singular format
830  for (i=1;i<=size(p2t)-2;i++)
831  {
832    if ((p2t[i]==":") and (p2t[i+1]=="=") and (p2t[i+2]=="["))
833    {
834      np2t=np2t+p2t[i]+p2t[i+1];
835      i=i+3;
836      while (p2t[i]!=":")
837      {
838        i=i+1;
839      }
840    }
841    else
842    {
843      if ((p2t[i]=="]") and (p2t[i+1]==";"))
844      {
845        np2t=np2t+p2t[i+1];
846        i=i+1;
847      }
848      else
849      {
850        np2t=np2t+p2t[i];
851      }
852    }
853  }
854  if (p2t[size(p2t)-1]=="]")
855  {
856    np2t=np2t+p2t[size(p2t)];
857  }
858  else
859  {
860    if (np2t[size(np2t)]!=";")
861    {
862      np2t=np2t+p2t[size(p2t)-1]+p2t[size(p2t)];
863    }
864  }
865  p2t=np2t;
866  np2t="";
867  // transform the points2triangs output of version 0.14 into Singular format
868  for (i=1;i<=size(p2t);i++)
869  {
870    if (p2t[i]=="=")
871    {
872      np2t=np2t+p2t[i]+"list(";
873      i++;
874    }
875    else
876    {
877      if (p2t[i]!=":")
878      {
879        if ((p2t[i]=="}") and (p2t[i+1]=="}"))
880        {
881          np2t=np2t+"))";
882          i++;
883        }
884        else
885        {
886          if (p2t[i]=="{")
887          {
888            np2t=np2t+"intvec(";
889          }
890          else
891          {
892            if (p2t[i]=="}")
893            {
894              np2t=np2t+")";
895            }
896            else
897            {
898              np2t=np2t+p2t[i];
899            }
900          }
901        }
902      }
903    }
904  }
905  list T;
906  execute(np2t);
907  // raise each index by one
908  for (i=1;i<=size(T);i++)
909  {
910    for (j=1;j<=size(T[i]);j++)
911    {
912      T[i][j]=T[i][j]+1;
913    }
914  }
915  return(T);
916}
917example
918{
919   "EXAMPLE:";
920   echo=2;
921   // the lattice points of the unit square in the plane
922   list polygon=intvec(0,0),intvec(0,1),intvec(1,0),intvec(1,1);
923   // the triangulations of this lattice point configuration are computed
924   list triang=triangulations(polygon);
925   triang;
926}
927
928proc secondaryPolytope (list polygon,list #)
929"USAGE:  secondaryPolytope(polygon[,#]); list polygon, list #
930ASSUME:  - polygon is a list of integer vectors of the same size representing the affine
931           coordinates of lattice points
932@*       - if the triangulations of the corresponding polygon have already been computed
933           with the procedure triangulations then these can be given as a second (optional)
934           argument in order to avoid doing this computation again
935PURPOSE: the procedure considers the marked polytope given as the convex hull of
936         the lattice points and with these lattice points as markings; it then
937         computes the lattice points of the secondary polytope given by this
938         marked polytope which correspond to the triangulations computed by
939         the procedure triangulations
940RETURN:  list, say L, such that:
941@*             L[1] = intmat, each row gives the affine coordinates of a lattice point
942                      in the secondary polytope given by the marked
943                      polytope corresponding to polygon
944@*             L[2] = the list of corresponding triangulations
945NOTE:    if the triangluations are not handed over as optional argument the procedure calls
946         for its computation of these triangulations the program points2triangs
947         from the program topcom by Joerg Rambau, Universitaet Bayreuth; it
948         therefore is necessary that this program is installed in order to use this
949         procedure; see
950@*       http://www.uni-bayreuth.de/departments/wirtschaftsmathematik/rambau/TOPCOM
951EXAMPLE: example secondaryPolytope;   shows an example"
952{
953  // compute the triangulations of the point configuration with points2triangs
954  if (size(#)==0)
955  {
956    list triangs=triangulations(polygon);
957  }
958  else
959  {
960    list triangs=#;
961  }
962  int i,j,k,l;
963  intmat N[2][2]; // is used to compute areas of triangles
964  intvec vertex;  // stores a point in the secondary polytope as intermediate result
965  int eintrag;
966  int halt;
967  intmat secpoly[size(triangs)][size(polygon)];   // stores all lattice points of the secondary polytope
968  // consider each triangulation and compute the corresponding point in the secondary polytope
969  for (i=1;i<=size(triangs);i++)
970  {
971    // for each triangulation we have to compute the coordinates corresponding to each marked point
972    for (j=1;j<=size(polygon);j++)
973    {
974      eintrag=0;
975      // for each marked point we have to consider all triangles in the triangulation
976      // which involve this particular point
977      for (k=1;k<=size(triangs[i]);k++)
978      {
979        halt=0;
980        for (l=1;(l<=3) and (halt==0);l++)
981        {
982          if (triangs[i][k][l]==j)
983          {
984            halt=1;
985            N[1,1]=polygon[triangs[i][k][3]][1]-polygon[triangs[i][k][1]][1];
986            N[1,2]=polygon[triangs[i][k][2]][1]-polygon[triangs[i][k][1]][1];
987            N[2,1]=polygon[triangs[i][k][3]][2]-polygon[triangs[i][k][1]][2];
988            N[2,2]=polygon[triangs[i][k][2]][2]-polygon[triangs[i][k][1]][2];
989            eintrag=eintrag+abs(det(N));
990          }
991        }
992      }
993      vertex[j]=eintrag;
994    }
995    secpoly[i,1..size(polygon)]=vertex;
996  }
997  return(list(secpoly,triangs));
998}
999example
1000{
1001   "EXAMPLE:";
1002   echo=2;
1003   // the lattice points of the unit square in the plane
1004   list polygon=intvec(0,0),intvec(0,1),intvec(1,0),intvec(1,1);
1005   // the secondary polytope of this lattice point configuration is computed
1006   list secpoly=secondaryPolytope(polygon);
1007   // the points in the secondary polytope
1008   print(secpoly[1]);
1009   // the corresponding triangulations
1010   secpoly[2];
1011}
1012
1013////////////////////////////////////////////////////////////////////////////////////////
1014/// PROCEDURES USING POLYMAKE AND TOPCOM
1015////////////////////////////////////////////////////////////////////////////////////////
1016
1017proc secondaryFan (list polygon,list #)
1018"USAGE:  secondaryFan(polygon[,#]); list polygon, list #
1019ASSUME:  - polygon is a list of integer vectors of the same size representing the affine
1020           coordinates of lattice points
1021@*       - if the triangulations of the corresponding polygon have already been computed
1022           with the procedure triangulations then these can be given as a second (optional)
1023           argument in order to avoid doing this computation again
1024PURPOSE: the procedure considers the marked polytope given as the convex hull of
1025         the lattice points and with these lattice points as markings; it then
1026         computes the lattice points of the secondary polytope given by this
1027         marked polytope which correspond to the triangulations computed by
1028         the procedure triangulations
1029RETURN:  list, the ith entry of L[1] contains information about the ith cone in the
1030               secondary fan of the polygon, i.e. the cone dual to the ith vertex of the
1031               secondary polytope
1032@*             L[1][i][1] = integer matrix representing the inequalities which describe the
1033                            cone dual to the ith vertex
1034@*             L[1][i][2] = a list which contains the inequalities represented by L[1][i][1]
1035                            as a list of strings, where we use the variables x(1),...,x(n)
1036@*             L[1][i][3] = only present if 'er' is set to 1; in that case it is an interger matrix
1037                            whose rows are the extreme rays of the cone
1038@*             L[2] = is an integer matrix whose rows span the linearity space of the fan,
1039                      i.e. the linear space which is contained in each cone
1040@*             L[3] = the secondary polytope in the format of the procedure polymakePolytope
1041@*             L[4] = the list of triangulations corresponding to the vertices
1042                      of the secondary polytope
1043NOTE:    - the procedure calls for its computation polymake by Ewgenij Gawrilow,
1044           TU Berlin and Michael Joswig, so it only works if polymake is installed;
1045           see http://www.math.tu-berlin.de/polymake/
1046@*       - in the optional argument # it is possible to hand over other names for the
1047           variables to be used -- be carful, the format must be correct and that is
1048           not tested, e.g. if you want the variable names to be u00,u10,u01,u11
1049           then you must hand over the string u11,u10,u01,u11
1050@*       - if the triangluations are not handed over as optional argument the procedure calls
1051           for its computation of these triangulations the program points2triangs
1052           from the program topcom by Joerg Rambau, Universitaet Bayreuth; it
1053           therefore is necessary that this program is installed in order to use this
1054           procedure; see
1055@*         http://www.uni-bayreuth.de/departments/wirtschaftsmathematik/rambau/TOPCOM
1056EXAMPLE: example secondaryFan;   shows an example"
1057{
1058  if (size(#)==0)
1059  {
1060    list triang=triangulations(polygon);
1061  }
1062  else
1063  {
1064    list triang=#[1];
1065  }
1066  list sp=secondaryPolytope(polygon,triang);
1067  list spp=polymakePolytope(sp[1]);
1068  list sf=normalFan(spp[1],spp[4],spp[3],1);
1069  return(list(sf[1],sf[2],spp,triang));
1070}
1071example
1072{
1073   "EXAMPLE:";
1074   echo=2;
1075   // the lattice points of the unit square in the plane
1076   list polygon=intvec(0,0),intvec(0,1),intvec(1,0),intvec(1,1);
1077   // the secondary polytope of this lattice point configuration is computed
1078   list secfan=secondaryFan(polygon);
1079   // the number of cones in the secondary fan of the polygon
1080   size(secfan[1]);
1081   // the inequalities of the first cone as matrix are:
1082   print(secfan[1][1][1]);
1083   // the inequalities of the first cone as string are:
1084   print(secfan[1][1][2]);
1085   // the rows of the following matrix are the extreme rays of the first cone:
1086   print(secfan[1][1][3]);
1087   // each cone contains the linearity space spanned by:
1088   print(secfan[2]);
1089   // the points in the secondary polytope
1090   print(secfan[3][1]);
1091   // the corresponding triangulations
1092   secfan[4];
1093}
1094
1095
1096////////////////////////////////////////////////////////////////////////////////////////
1097/// PROCEDURES CONCERNED WITH PLANAR POLYGONS
1098////////////////////////////////////////////////////////////////////////////////////////
1099
1100proc cycleLength (list boundary,intvec interior)
1101"USAGE:  cycleLength(boundary,interior); list boundary, intvec interior
1102ASSUME:  boundary is a list of integer vectors describing a cycle in some convex lattice
1103         polygon around the lattice point interior ordered clock wise
1104RETURN:  string, the cycle length of the corresponding cycle in the dual tropical curve
1105EXAMPLE: example cycleLength;   shows an example"
1106{
1107  int j;
1108  // create a ring whose variables are indexed by the points in boundary resp. by interior
1109  string rst="ring cyclering=0,(u"+string(interior[1])+string(interior[2]);
1110  for (j=1;j<=size(boundary);j++)
1111  {
1112    rst=rst+",u"+string(boundary[j][1])+string(boundary[j][2]);
1113  }
1114  rst=rst+"),lp;";
1115  execute(rst);
1116  // add the first and second point at the end of boundary
1117  boundary[size(boundary)+1]=boundary[1];
1118  boundary[size(boundary)+1]=boundary[2];
1119  poly cl,summand; // takes the cycle length
1120  matrix N1[2][2]; // used to compute the area of a triangle
1121  matrix N2[2][2]; // used to compute the area of a triangle
1122  matrix N3[2][2]; // used to compute the area of a triangle
1123  // for each original point in the boundary compute its contribution to the cycle
1124  for (j=2;j<=size(boundary)-1;j++)
1125  {
1126    N1=boundary[j-1]-interior,boundary[j]-interior;
1127    N2=boundary[j]-interior,boundary[j+1]-interior;
1128    N3=boundary[j+1]-interior,boundary[j-1]-interior;
1129    execute("summand=-u"+string(boundary[j][1])+string(boundary[j][2])+"+u"+string(interior[1])+string(interior[2])+";");
1130    summand=summand*(det(N1)+det(N2)+det(N3))/(det(N1)*det(N2));
1131    cl=cl+summand;
1132  }
1133  return(string(cl));
1134}
1135example
1136{
1137   "EXAMPLE:";
1138   echo=2;
1139   // the integer vectors in boundary are lattice points on the boundary
1140   // of a convex lattice polygon in the plane
1141   list boundary=intvec(0,0),intvec(0,1),intvec(0,2),intvec(2,2),
1142                 intvec(2,1),intvec(2,0);
1143   // interior is a lattice point in the interior of this lattice polygon
1144   intvec interior=1,1;
1145   // compute the general cycle length of a cycle of the corresponding cycle
1146   // in the dual tropical curve, note that (0,1) and (2,1) do not contribute
1147   cycleLength(boundary,interior);
1148}
1149
1150proc splitPolygon (list markings)
1151"USAGE:  splitPolygon (markings);  markings list
1152ASSUME:  markings is a list of integer vectors representing lattice points in the plane
1153         which we consider as the marked points of the convex lattice polytope spanned by them
1154PURPOSE: split the marked points in the vertices, the points on the facets which are not vertices,
1155         and the interior points
1156RETURN:  list, L consisting of three lists
1157@*             L[1]    : represents the vertices the polygon ordered clockwise
1158@*                       L[1][i][1] = intvec, the coordinates of the ith vertex
1159@*                       L[1][i][2] = int, the position of L[1][i][1] in markings
1160@*             L[2][i] : represents the lattice points on the facet of the polygon with
1161                         endpoints L[1][i] and L[1][i+1] (i considered modulo size(L[1]))
1162@*                       L[2][i][j][1] = intvec, the coordinates of the jth lattice point on that facet
1163@*                       L[2][i][j][2] = int, the position of L[2][i][j][1] in markings
1164@*             L[3]    : represents the interior lattice points of the polygon
1165@*                       L[3][i][1] = intvec, the coordinates of the ith interior point
1166@*                       L[3][i][2] = int, the position of L[3][i][1] in markings
1167EXAMPLE: example splitPolygon;   shows an example"
1168{
1169  list vert; // stores the result
1170  // compute the boundary of the polygon in an oriented way
1171  list pb=findOrientedBoundary(markings);
1172  // the vertices are just the second entry of pb
1173  vert[1]=pb[2];
1174  int i,j,k;      // indices
1175  list boundary;  // stores the points on the facets of the polygon which are not vertices
1176  // append to the boundary points as well as to the vertices the first vertex a second time
1177  pb[1]=pb[1]+list(pb[1][1]);
1178  pb[2]=pb[2]+list(pb[2][1]);
1179  // for each vertex find all points on the facet of the polygon with this vertex
1180  // and the next vertex as endpoints
1181  int z=2;
1182  for (i=1;i<=size(vert[1]);i++)
1183  {
1184    j=1;
1185    list facet; // stores the points on this facet which are not vertices
1186    // while the next vertex is not reached, store the boundary lattice point
1187    while (pb[1][z]!=pb[2][i+1])
1188    {
1189      facet[j]=pb[1][z];
1190      j++;
1191      z++;
1192    }
1193    // store the points on the ith facet as boundary[i]
1194    boundary[i]=facet;
1195    kill facet;
1196    z++;
1197  }
1198  // store the information on the boundary in vert[2]
1199  vert[2]=boundary;
1200  // find the remaining points in the input which are not on the boundary by checking
1201  // for each point in markings if it is contained in pb[1]
1202  list interior=markings;
1203  for (i=size(interior);i>=1;i--)
1204  {
1205    for (j=1;j<=size(pb[1])-1;j++)
1206    {
1207      if (interior[i]==pb[1][j])
1208      {
1209        interior=delete(interior,i);
1210        j=size(pb[1]);
1211      }
1212    }
1213  }
1214  // store the interior points in vert[3]
1215  vert[3]=interior;
1216  // add to each point in vert the index which it gets from its position in the input markings;
1217  // do so for ver[1]
1218  for (i=1;i<=size(vert[1]);i++)
1219  {
1220    j=1;
1221    while (markings[j]!=vert[1][i])
1222    {
1223      j++;
1224    }
1225    vert[1][i]=list(vert[1][i],j);
1226  }
1227  // do so for ver[2]
1228  for (i=1;i<=size(vert[2]);i++)
1229  {
1230    for (k=1;k<=size(vert[2][i]);k++)
1231    {
1232      j=1;
1233      while (markings[j]!=vert[2][i][k])
1234      {
1235        j++;
1236      }
1237      vert[2][i][k]=list(vert[2][i][k],j);
1238    }
1239  }
1240  // do so for ver[3]
1241  for (i=1;i<=size(vert[3]);i++)
1242  {
1243    j=1;
1244    while (markings[j]!=vert[3][i])
1245    {
1246      j++;
1247    }
1248    vert[3][i]=list(vert[3][i],j);
1249  }
1250  return(vert);
1251}
1252example
1253{
1254   "EXAMPLE:";
1255   echo=2;
1256   // the lattice polygon spanned by the points (0,0), (3,0) and (0,3)
1257   // with all integer points as markings
1258   list polygon=intvec(1,1),intvec(3,0),intvec(2,0),intvec(1,0),
1259                intvec(0,0),intvec(2,1),intvec(0,1),intvec(1,2),
1260                intvec(0,2),intvec(0,3);
1261   // split the polygon in its vertices, its facets and its interior points
1262   list sp=splitPolygon(polygon);
1263   // the vertices
1264   sp[1];
1265   // the points on facets which are not vertices
1266   sp[2];
1267   // the interior points
1268   sp[3];
1269}
1270
1271
1272proc eta (list triang,list polygon)
1273"USAGE:  eta(triang,polygon);  triang, polygon list
1274ASSUME:  polygon has the format of the output of splitPolygon, i.e. it is a list with three
1275         entries describing a convex lattice polygon in the following way:
1276@*       polygon[1] : is a list of lists; for each i the entry polygon[1][i][1] is a lattice point which is
1277                      a vertex of the lattice polygon, and polygon[1][i][2] is an integer assigned to
1278                      this lattice point as identifying index
1279@*       polygon[2] : is a list of lists; for each vertex of the polygon, i.e. for each entry in polygon[1],
1280                      it contains a list polygon[2][i], which contains the lattice points on the facet
1281                      with endpoints polygon[1][i] and polygon[1][i+1] - i considered mod size(polygon[1]);
1282                      each such lattice point contributes an entry polygon[2][i][j][1] which is an integer
1283                      vector giving the coordinate of the lattice point and an entry polygon[2][i][j][2]
1284                      which is the identifying index
1285@*       polygon[3] : is a list of lists, where each entry corresponds to a lattice point in the
1286                      interior of the polygon, with polygon[3][j][1] being the coordinates of the point
1287                      and polygon[3][j][2] being the identifying index;
1288@*       triang is a list of integer vectors all of size three describing a triangulation of the
1289         polygon described by polygon; if an entry of triang is the vector (i,j,k) then the triangle
1290         is build by the vertices with indices i, j and k
1291RETURN:  intvec, the integer vector eta describing that vertex of the Newton polytope discriminant
1292                 of the polygone whose dual cone in the Groebner fan contains the cone of the
1293                 secondary fan of the polygon corresponding to the given triangulation
1294NOTE:    for a better description of eta see either Gelfand, Kapranov, Zelevinski: Discriminants,
1295         Resultants and multidimensional Determinants. Chapter 10.
1296EXAMPLE: example eta;   shows an example"
1297{
1298  int i,j,k,l,m,n; // index variables
1299  list ordpolygon;   // stores the lattice points in the order used in the triangulation
1300  list triangarea; // stores the areas of the triangulations
1301  intmat N[2][2];  // used to compute triangle areas
1302  // 1) store the lattice points in the order used in the triangulation
1303  // go first through all vertices of the polytope
1304  for (j=1;j<=size(polygon[1]);j++)
1305  {
1306    ordpolygon[polygon[1][j][2]]=polygon[1][j][1];
1307  }
1308  // then consider all inner points
1309  for (j=1;j<=size(polygon[3]);j++)
1310  {
1311    ordpolygon[polygon[3][j][2]]=polygon[3][j][1];
1312  }
1313  // finally consider all lattice points on the boundary which are not vertices
1314  for (j=1;j<=size(polygon[2]);j++)
1315  {
1316    for (i=1;i<=size(polygon[2][j]);i++)
1317    {
1318      ordpolygon[polygon[2][j][i][2]]=polygon[2][j][i][1];
1319    }
1320  }
1321  // 2) compute for each triangle in the triangulation the area of the triangle
1322  for (i=1;i<=size(triang);i++)
1323  {
1324    // Note that the ith lattice point in orderedpolygon has the number i-1 in the triangulation!
1325    N=ordpolygon[triang[i][1]]-ordpolygon[triang[i][2]],ordpolygon[triang[i][1]]-ordpolygon[triang[i][3]];
1326    triangarea[i]=abs(det(N));
1327  }
1328  intvec ETA;        // stores the eta_ij
1329  int etaij;         // stores the part of eta_ij during computations which comes from triangle areas
1330  int seitenlaenge;  // stores the part of eta_ij during computations which comes from boundary facets
1331  list seiten;       // stores the lattice points on facets of the polygon
1332  intvec v;          // used to compute a facet length
1333  // 3) store first in seiten[i] all lattice points on the facet connecting the ith vertex,
1334  //    i.e. polygon[1][i], with the i+1st vertex, i.e. polygon[1][i+1], where we replace i+1
1335  //    1 if i=size(polygon[1]);
1336  //    then append the last entry of seiten once more at the very beginning of seiten, so
1337  //    that the index is shifted by one
1338  for (i=1;i<=size(polygon[1]);i++)
1339  {
1340    if (i<size(polygon[1]))
1341    {
1342      seiten[i]=list(polygon[1][i])+polygon[2][i]+list(polygon[1][i+1]);
1343    }
1344    else
1345    {
1346      seiten[i]=list(polygon[1][i])+polygon[2][i]+list(polygon[1][1]);
1347    }
1348  }
1349  seiten=insert(seiten,seiten[size(seiten)],0);
1350  // 4) compute the eta_ij for all vertices of the polygon
1351  for (j=1;j<=size(polygon[1]);j++)
1352  {
1353    // the vertex itself contributes a 1
1354    etaij=1;
1355    // check for each triangle in the triangulation ...
1356    for (k=1;k<=size(triang);k++)
1357    {
1358      // ... if the vertex is actually a vertex of the triangle ...
1359      if ((polygon[1][j][2]==triang[k][1]) or (polygon[1][j][2]==triang[k][2]) or (polygon[1][j][2]==triang[k][3]))
1360      {
1361        // ... if so, add the area of the triangle to etaij
1362        etaij=etaij+triangarea[k];
1363        // then check if that triangle has a facet which is contained in one of the
1364        // two facets of the polygon which are adjecent to the given vertex ...
1365        // these two facets are seiten[j] and seiten[j+1]
1366        for (n=j;n<=j+1;n++)
1367        {
1368          // check for each lattice point in the facet of the polygon ...
1369          for (l=1;l<=size(seiten[n]);l++)
1370          {
1371            // ... and for each lattice point in the triangle ...
1372            for (m=1;m<=size(triang[k]);m++)
1373            {
1374              // ... if they coincide and are not the vertex itself ...
1375              if ((seiten[n][l][2]==triang[k][m]) and (seiten[n][l][2]!=polygon[1][j][2]))
1376              {
1377                // if so, then compute the vector pointing from this lattice point to the vertex
1378                v=polygon[1][j][1]-seiten[n][l][1];
1379                // and the lattice length of this vector has to be subtracted from etaij
1380                etaij=etaij-abs(gcd(v[1],v[2]));
1381              }
1382            }
1383          }
1384        }
1385      }
1386    }
1387    // store etaij in the list
1388    ETA[polygon[1][j][2]]=etaij;
1389  }
1390  // 5) compute the eta_ij for all lattice points on the facets of the polygon which are not vertices,
1391  //    these are the lattice points in polygon[2][1] to polygon[2][size(polygon[1])]
1392  for (i=1;i<=size(polygon[2]);i++)
1393  {
1394    for (j=1;j<=size(polygon[2][i]);j++)
1395    {
1396      // initialise etaij
1397      etaij=0;
1398      // initialise seitenlaenge
1399      seitenlaenge=0;
1400      // check for each triangle in the triangulation ...
1401      for (k=1;k<=size(triang);k++)
1402      {
1403        // ... if the vertex is actually a vertex of the triangle ...
1404        if ((polygon[2][i][j][2]==triang[k][1]) or (polygon[2][i][j][2]==triang[k][2]) or (polygon[2][i][j][2]==triang[k][3]))
1405        {
1406          // ... if so, add the area of the triangle to etaij
1407          etaij=etaij+triangarea[k];
1408          // then check if that triangle has a facet which is contained in the
1409          // facet of the polygon which contains the lattice point in question,
1410          // this is the facet seiten[i+1];
1411          // check for each lattice point in the facet of the polygon ...
1412          for (l=1;l<=size(seiten[i+1]);l++)
1413          {
1414            // ... and for each lattice point in the triangle ...
1415            for (m=1;m<=size(triang[k]);m++)
1416            {
1417              // ... if they coincide and are not the vertex itself ...
1418              if ((seiten[i+1][l][2]==triang[k][m]) and (seiten[i+1][l][2]!=polygon[2][i][j][2]))
1419              {
1420                // if so, then compute the vector pointing from this lattice point to the vertex
1421                v=polygon[2][i][j][1]-seiten[i+1][l][1];
1422                // and the lattice length of this vector contributes to seitenlaenge
1423                seitenlaenge=seitenlaenge+abs(gcd(v[1],v[2]));
1424              }
1425            }
1426          }
1427        }
1428      }
1429      // if the lattice point was a vertex of any triangle in the triangulation ...
1430      if (etaij!=0)
1431      {
1432        // then eta_ij is the sum of the triangle areas minus seitenlaenge
1433        ETA[polygon[2][i][j][2]]=etaij-seitenlaenge;
1434      }
1435      else
1436      {
1437        // otherwise it is just zero
1438        ETA[polygon[2][i][j][2]]=0;
1439      }
1440    }
1441  }
1442  // 4) compute the eta_ij for all inner lattice points of the polygon
1443  for (j=1;j<=size(polygon[3]);j++)
1444  {
1445    // initialise etaij
1446    etaij=0;
1447    // check for each triangle in the triangulation ...
1448    for (k=1;k<=size(triang);k++)
1449    {
1450      // ... if the vertex is actually a vertex of the triangle ...
1451      if ((polygon[3][j][2]==triang[k][1]) or (polygon[3][j][2]==triang[k][2]) or (polygon[3][j][2]==triang[k][3]))
1452      {
1453        // ... if so, add the area of the triangle to etaij
1454        etaij=etaij+triangarea[k];
1455      }
1456    }
1457    // store etaij in ETA
1458    ETA[polygon[3][j][2]]=etaij;
1459  }
1460  return(ETA);
1461}
1462example
1463{
1464   "EXAMPLE:";
1465   echo=2;
1466   // the lattice polygon spanned by the points (0,0), (3,0) and (0,3)
1467   // with all integer points as markings
1468   list polygon=intvec(1,1),intvec(3,0),intvec(2,0),intvec(1,0),
1469                intvec(0,0),intvec(2,1),intvec(0,1),intvec(1,2),
1470                intvec(0,2),intvec(0,3);
1471   // split the polygon in its vertices, its facets and its interior points
1472   list sp=splitPolygon(polygon);
1473   // define a triangulation by connecting the only interior point
1474   //        with the vertices
1475   list triang=intvec(1,2,5),intvec(1,5,10),intvec(1,5,10);
1476   // compute the eta-vector of this triangulation
1477   eta(triang,sp);
1478}
1479
1480proc findOrientedBoundary (list polygon)
1481"USAGE:      findOrientedBoundary(polygon); polygon list
1482ASSUME:      polygon is a list of integer vectors defining integer lattice points in the plane
1483RETURN:      list, l with the followin interpretation
1484@*                 l[1] = list of integer vectors such that the polygonal path defined by
1485                          these is the boundary of the convex hull of the lattice points in polygon
1486@*                 l[2] = list, the redundant points in l[1] have been removed
1487EXAMPLE:     example findOrientedBoundary;   shows an example"
1488{
1489  // Order the vertices such that passing from one to the next we travel along
1490  // the boundary of the convex hull of the vertices clock wise
1491  int d,k,i,j;
1492  intmat D[2][2];
1493  /////////////////////////////////////
1494  // Treat first the pathological cases that the polygon is not two-dimensional:
1495  /////////////////////////////////////
1496  // if the polygon is empty or only one point or a line segment of two points
1497  if (size(polygon)<=2)
1498  {
1499    return(list(polygon,polygon));
1500  }
1501  // check is the polygon is only a line segment given by more than two points;
1502  // for this first compute sum of the absolute values of the determinants of the matrices whose
1503  // rows are the vectors pointing from the first to the second point and from the
1504  // the first point to the ith point for i=3,...,size(polygon); if this sum is zero
1505  // then the polygon is a line segment and we have to find its end points
1506  d=0;
1507  for (i=3;i<=size(polygon);i++)
1508  {
1509    D=polygon[2]-polygon[1],polygon[i]-polygon[1];
1510    d=d+abs(det(D));
1511  }
1512  if (d==0) // then polygon is a line segment
1513  {
1514    intmat laenge[size(polygon)][size(polygon)];
1515    intvec mp;
1516    //   for this collect first all vectors pointing from one lattice point to the next,
1517    //   compute their pairwise angles and their lengths
1518    for (i=1;i<=size(polygon)-1;i++)
1519    {
1520      for (j=i+1;j<=size(polygon);j++)
1521      {
1522        mp=polygon[i]-polygon[j];
1523        laenge[i,j]=abs(gcd(mp[1],mp[2]));
1524      }
1525    }
1526    mp=maxPosInIntmat(laenge);
1527    list endpoints=polygon[mp[1]],polygon[mp[2]];
1528    intvec abstand;
1529    for (i=1;i<=size(polygon);i++)
1530    {
1531      abstand[i]=0;
1532      if (i<mp[1])
1533      {
1534        abstand[i]=laenge[i,mp[1]];
1535      }
1536      if (i>mp[1])
1537      {
1538        abstand[i]=laenge[mp[1],i];
1539      }
1540    }
1541    polygon=sortlistbyintvec(polygon,abstand);
1542    return(list(polygon,endpoints));
1543  }
1544  ///////////////////////////////////////////////////////////////
1545  list orderedvertices;  // stores the vertices in an ordered way
1546  list minimisedorderedvertices;  // stores the vertices in an ordered way; redundant ones removed
1547  list comparevertices; // stores vertices which should be compared to the testvertex
1548  orderedvertices[1]=polygon[1]; // set the starting vertex
1549  minimisedorderedvertices[1]=polygon[1]; // set the starting vertex
1550  intvec testvertex=polygon[1];  // the vertex to which the others have to be compared
1551  intvec startvertex=polygon[1]; // keep the starting vertex to test, when the end is reached
1552  int endtest;                    // is set to one, when the end is reached
1553  int startvertexfound;// is 1, once for some testvertex a candidate for the next vertex has been found
1554  polygon=delete(polygon,1);    // delete the testvertex
1555  intvec v,w;
1556  int l=1;  // counts the vertices
1557  // the basic idea is that a vertex can be the next one on the boundary if all other vertices
1558  // ly to the right of the vector v pointing from the testvertex to this one; this can be tested
1559  // by checking if the determinant of the 2x2-matrix with first column v and second column the vector w,
1560  // pointing from the testvertex to the new vertex, is non-positive; if this is the case for all
1561  // new vertices, then the one in consideration is a possible choice for the next vertex on the boundary
1562  // and it is stored in naechste; we can then order the candidates according to their distance from
1563  // the testvertex; then they occur on the boundary in that order!
1564  while (endtest==0)
1565  {
1566    list naechste;  // stores the possible choices for the next vertex
1567    k=1;
1568    for (i=1;i<=size(polygon);i++)
1569    {
1570      d=0;  // stores the value of the determinant of (v,w)
1571      v=polygon[i]-testvertex; // points from the testvertex to the ith vertex
1572      comparevertices=delete(polygon,i); // we needn't compare v to itself
1573      // we should compare v to the startvertex-testvertex; in the first calling of the loop
1574      // this is irrelevant since the difference will be zero; however, later on it will
1575      // be vital, since we delete the vertices which we have already tested from the list
1576      // of all vertices, and when all vertices on the boundary have been found we would
1577      // therefore find a vertex in the interior as candidate; but always testing against
1578      // the starting vertex, this can not happen
1579      comparevertices[size(comparevertices)+1]=startvertex;
1580      for (j=1;(j<=size(comparevertices)) and (d<=0);j++)
1581      {
1582        w=comparevertices[j]-testvertex; // points form the testvertex to the jth vertex
1583        D=v,w;
1584        d=det(D);
1585      }
1586      if (d<=0) // if all determinants are non-positive, then the ith vertex is a candidate
1587      {
1588        naechste[k]=list(polygon[i],i,scalarproduct(v,v)); // we store the vertex, its position, and its
1589        k++;                                                // distance from the testvertex
1590      }
1591    }
1592    if (size(naechste)>0) // then a candidate for the next vertex has been found
1593    {
1594      startvertexfound=1; // at least once a candidate has been found
1595      naechste=sortlist(naechste,3);  //we order the candidates according to their distance from testvertex;
1596      for (j=1;j<=size(naechste);j++) // then we store them in this order in orderedvertices
1597      {
1598        l++;
1599        orderedvertices[l]=naechste[j][1];
1600      }
1601      testvertex=naechste[size(naechste)][1];  // we store the last one as next testvertex;
1602      minimisedorderedvertices[size(minimisedorderedvertices)+1]=testvertex; // store the next corner of NSD
1603      naechste=sortlist(naechste,2);           // then we reorder the vertices according to their position
1604      for (j=size(naechste);j>=1;j--)          // and we delete them from the vertices
1605      {
1606        polygon=delete(polygon,naechste[j][2]);
1607      }
1608    }
1609    else // that means either that the vertex was inside the polygon,
1610    {    // or that we have reached the last vertex on the boundary of the polytope
1611      if (startvertexfound==0) // the vertex was in the interior; we delete it and start all over again
1612      {
1613        orderedvertices[1]=polygon[1];
1614        minimisedorderedvertices[1]=polygon[1];
1615        testvertex=polygon[1];
1616        startvertex=polygon[1];
1617        polygon=delete(polygon,1);
1618      }
1619      else // we have reached the last vertex on the boundary of the polytope and can stop
1620      {
1621        endtest=1;
1622      }
1623    }
1624    kill naechste;
1625  }
1626  // test if the first vertex in minimisedorderedvertices is on the same line with the second and
1627  // the last, i.e. if we started our search in the middle of a face; if so, delete it
1628  v=minimisedorderedvertices[2]-minimisedorderedvertices[1];
1629  w=minimisedorderedvertices[size(minimisedorderedvertices)]-minimisedorderedvertices[1];
1630  D=v,w;
1631  if (det(D)==0)
1632  {
1633    minimisedorderedvertices=delete(minimisedorderedvertices,1);
1634  }
1635  // test if the first vertex in minimisedorderedvertices is on the same line with the two
1636  // last ones, i.e. if we started our search at the end of a face; if so, delete it
1637  v=minimisedorderedvertices[size(minimisedorderedvertices)-1]-minimisedorderedvertices[1];
1638  w=minimisedorderedvertices[size(minimisedorderedvertices)]-minimisedorderedvertices[1];
1639  D=v,w;
1640  if (det(D)==0)
1641  {
1642    minimisedorderedvertices=delete(minimisedorderedvertices,size(minimisedorderedvertices));
1643  }
1644  return(list(orderedvertices,minimisedorderedvertices));
1645}
1646example
1647{
1648   "EXAMPLE:";
1649   echo=2;
1650// the following lattice points in the plane define a polygon
1651   list polygon=intvec(0,0),intvec(3,1),intvec(1,0),intvec(2,0),
1652                intvec(1,1),intvec(3,2),intvec(1,2),intvec(2,3),
1653                intvec(2,4);
1654// we compute its boundary
1655   list boundarypolygon=findOrientedBoundary(polygon);
1656// the points on the boundary ordered clockwise are boundarypolygon[1]
1657   boundarypolygon[1];
1658// the vertices of the boundary are boundarypolygon[2]
1659   boundarypolygon[2];
1660}
1661
1662
1663proc cyclePoints (list triang,list points,int pt)
1664"USAGE:      cyclePoints(triang,points,pt)  triang,points list, pt int
1665ASSUME:      - points is a list of integer vectors describing the lattice points of a marked polygon;
1666@*           - triang is a list of integer vectors describing a triangulation of the marked polygon
1667               in the sense that an integer vector of the form (i,j,k) describes the triangle formed
1668               by polygon[i], polygon[j] and polygon[k];
1669@*           - pt is an integer between 1 and size(points), singling out a lattice point among
1670               the marked points
1671PURPOSE:     consider the convex lattice polygon, say P, spanned by all lattice points in points which
1672             in the triangulation triang are connected to the point points[pt]; the procedure
1673             computes all marked points in points which ly on the boundary of that polygon, ordered
1674             clockwise
1675RETURN:      list, of integer vectors which are the coordinates of the lattice points on
1676                   the boundary of the above mentioned polygon P, if this polygon is not the
1677                   empty set (that would be the case if points[pt] is not a vertex of any
1678                   triangle in the triangulation); otherwise return the empty list
1679EXAMPLE:     example cyclePoints;   shows an example"
1680{
1681  int i,j; // indices
1682  list v;  // saves the indices of lattice points connected to the interior point in the triangulation
1683  // save all points in triangulations containing pt in v
1684  for (i=1;i<=size(triang);i++)
1685  {
1686    if ((triang[i][1]==pt) or (triang[i][2]==pt) or (triang[i][3]==pt))
1687    {
1688      j++;
1689      v[3*j-2]=triang[i][1];
1690      v[3*j-1]=triang[i][2];
1691      v[3*j]=triang[i][3];
1692    }
1693  }
1694  if (size(v)==0)
1695  {
1696    return(list());
1697  }
1698  // remove pt itself and redundancies in v
1699  for (i=size(v);i>=1;i--)
1700  {
1701    j=1;
1702    while ((j<i) and (v[i]!=v[j]))
1703    {
1704      j++;
1705    }
1706    if ((j<i) or (v[i]==pt))
1707    {
1708      v=delete(v,i);
1709    }
1710  }
1711  // save in pts the coordinates of the points with indices in v
1712  list pts;
1713  for (i=1;i<=size(v);i++)
1714  {
1715    pts[i]=points[v[i]];
1716  }
1717  // consider the convex polytope spanned by the points in pts, find the points on the
1718  // boundary and order them clockwise
1719  return(findOrientedBoundary(pts)[1]);
1720}
1721example
1722{
1723   "EXAMPLE:";
1724   echo=2;
1725   // the lattice polygon spanned by the points (0,0), (3,0) and (0,3)
1726   // with all integer points as markings
1727   list points=intvec(1,1),intvec(3,0),intvec(2,0),intvec(1,0),
1728               intvec(0,0),intvec(2,1),intvec(0,1),intvec(1,2),
1729               intvec(0,2),intvec(0,3);
1730   // define a triangulation
1731   list triang=intvec(1,2,5),intvec(1,5,7),intvec(1,7,9),intvec(8,9,10),
1732               intvec(1,8,9),intvec(1,2,8);
1733   // compute the points connected to (1,1) in triang
1734   cyclePoints(triang,points,1);
1735}
1736
1737proc latticeArea (list polygon)
1738"USAGE:  latticeArea(polygon);   polygon list
1739ASSUME:  polygon is a list of integer vectors in the plane
1740RETURN:  int, the lattice area of the convex hull of the lattice points in polygon,
1741              i.e. twice the Euclidean area
1742EXAMPLE: example polygonlatticeArea;   shows an example"
1743{
1744  list pg=findOrientedBoundary(polygon)[2];
1745  int area;
1746  intmat M[2][2];
1747  for (int i=2;i<=size(pg)-1;i++)
1748  {
1749    M[1,1..2]=pg[i]-pg[1];
1750    M[2,1..2]=pg[i+1]-pg[1];
1751    area=area+abs(det(M));
1752  }
1753  return(area);
1754}
1755example
1756{
1757   "EXAMPLE:";
1758   echo=2;
1759   // define a polygon with lattice area 5
1760   list polygon=intvec(1,2),intvec(1,0),intvec(2,0),intvec(1,1),
1761                intvec(2,1),intvec(0,0);
1762   latticeArea(polygon);
1763}
1764
1765proc picksFormula (list polygon)
1766"USAGE:  picksFormula(polygon);   polygon list
1767ASSUME:  polygon is a list of integer vectors in the plane and consider their convex hull C
1768RETURN:  list, L of three integersthe
1769@*             L[1] : the lattice area of C, i.e. twice the Euclidean area
1770@*             L[2] : the number of lattice points on the boundary of C
1771@*             L[3] : the number of interior lattice points of C
1772NOTE:    the integers in L are related by Pick's formula, namely: L[1]=L[2]+2*L[3]-2
1773EXAMPLE: example picksFormula;   shows an example"
1774{
1775  list pg=findOrientedBoundary(polygon)[2];
1776  int area,bdpts,i;
1777  intmat M[2][2];
1778  // compute the lattice area of the polygon, i.e. twice the Euclidean area
1779  for (i=2;i<=size(pg)-1;i++)
1780  {
1781    M[1,1..2]=pg[i]-pg[1];
1782    M[2,1..2]=pg[i+1]-pg[1];
1783    area=area+abs(det(M));
1784  }
1785  // compute the number of lattice points on the boundary
1786  intvec edge;
1787  pg[size(pg)+1]=pg[1];
1788  for (i=1;i<=size(pg)-1;i++)
1789  {
1790    edge=pg[i]-pg[i+1];
1791    bdpts=bdpts+abs(gcd(edge[1],edge[2]));
1792  }
1793  // Pick's formula says that the lattice area A, the number g of interior points and
1794  // the number b of boundary points are connected by the formula: A=b+2g-2
1795  return(list(area,bdpts,(area-bdpts+2)/2));
1796}
1797example
1798{
1799   "EXAMPLE:";
1800   echo=2;
1801   // define a polygon with lattice area 5
1802   list polygon=intvec(1,2),intvec(1,0),intvec(2,0),intvec(1,1),
1803                intvec(2,1),intvec(0,0);
1804   list pick=picksFormula(polygon);
1805   // the lattice area of the polygon is:
1806   pick[1];
1807   // the number of lattice points on the boundary is:
1808   pick[2];
1809   // the number of interior lattice points is:
1810   pick[3];
1811   // the number's are related by Pick's formula:
1812   pick[1]-pick[2]-2*pick[3]+2;
1813}
1814
1815proc ellipticNF (list polygon)
1816"USAGE:  ellipticNF(polygon);   polygon list
1817ASSUME:  polygon is a list of integer vectors in the plane such that their convex hull C
1818         has precisely one interior lattice point; i.e. C is the Newton polygon of an
1819         elliptic curve
1820PURPOSE: compute the normal form of the polygon with respect to the unimodular affine
1821         transformations T=A*x+v; there are sixteen different normal forms
1822         (see e.g. Bjorn Poonen, Fernando Rodriguez-Villegas: Lattice Polygons and
1823                   the number 12.  Amer. Math. Monthly  107  (2000),  no. 3, 238--250.)
1824RETURN:  list, L such that
1825@*             L[1] : list whose entries are the vertices of the normal form of the polygon
1826@*             L[2] : the matrix A of the unimodular transformation
1827@*             L[3] : the translation vector v of the unimodular transformation
1828@*             L[4] : list such that the ith entry is the image of polygon[i] under the
1829                      unimodular transformation T
1830EXAMPLE: example ellipticNF;   shows an example"
1831{
1832  int i;            // index
1833  intvec edge;      // stores the vector of an edge
1834  intvec boundary;  // stores the lattice lengths of the edges of the Newton cylce
1835  // find the vertices of the Newton cycle and order it clockwise
1836  list pg=findOrientedBoundary(polygon)[2];
1837  // check if there is precisely one interior point in the Newton polygon
1838  if (picksFormula(pg)[3]!=1)
1839  {
1840    ERROR("The polygon has not precisely one interior point!");
1841  }
1842  // insert the first vertex at the end once again
1843  pg[size(pg)+1]=pg[1];
1844  // compute the number of lattice points on each edge
1845  for (i=1;i<=size(pg)-1;i++)
1846  {
1847    edge=pg[i]-pg[i+1];
1848    boundary[i]=1+abs(gcd(edge[1],edge[2]));
1849  }
1850  // store the values of boundary once more adding the first two at the end
1851  intvec tboundary=boundary,boundary[1],boundary[2];
1852  // sort boundary in an asecending way
1853  intvec sbd=sortintvec(boundary);
1854  // find the first edge having the maximal number of lattice points
1855  int max=maxPosInIntvec(boundary);
1856  // some computations have to be done over the rationals
1857  ring transformationring=0,x,lp;
1858  intvec trans;    // stores the vector by which we have to translate the polygon
1859  intmat A[2][2];  // stores the matrix by which we have to transform the polygon
1860  matrix M[3][3];  // stores the projective coordinates of the points which are to be transformed
1861  matrix N[3][3];  // stores the projective coordinates of the points to which M is to be transformed
1862  intmat T[3][3];  // stores the unimodular affine transformation in projective form
1863  // add the second point of pg once again at the end
1864  pg=insert(pg,pg[2],size(pg));
1865  // if there is only one edge which has the maximal number of lattice points, then M should be:
1866  M=pg[max],1,pg[max+1],1,pg[max+2],1;
1867  // consider the 16 different cases which can occur:
1868  // Case 1:
1869  if (sbd==intvec(2,2,2))
1870  {
1871    N=0,1,1,1,2,1,2,0,1;
1872  }
1873  // Case 2:
1874  if (sbd==intvec(2,2,3))
1875  {
1876    N=2,0,1,0,0,1,1,2,1;
1877  }
1878  // Case 3:
1879  if (sbd==intvec(2,3,4))
1880  {
1881    // here the orientation of the Newton polygon is important !
1882    if (tboundary[max+1]==3)
1883    {
1884      N=3,0,1,0,0,1,0,2,1;
1885    }
1886    else
1887    {
1888      N=0,0,1,3,0,1,0,2,1;
1889    }
1890  }
1891  // Case 4:
1892  if (sbd==intvec(3,3,5))
1893  {
1894    N=4,0,1,0,0,1,0,2,1;
1895  }
1896  // Case 5:
1897  if (sbd==intvec(4,4,4))
1898  {
1899    N=3,0,1,0,0,1,0,3,1;
1900  }
1901  // Case 6+7:
1902  if (sbd==intvec(2,2,2,2))
1903  {
1904    // there are two different polygons which has four edges all of length 2,
1905    // but only one of them has two edges whose direction vectors form a matrix
1906    // of determinant 3
1907    A=pg[1]-pg[2],pg[3]-pg[2];
1908    while ((max<4) and (det(A)!=3))
1909    {
1910      max++;
1911      A=pg[max]-pg[max+1],pg[max+2]-pg[max+1];
1912    }
1913    // Case 6:
1914    if (det(A)==3)
1915    {
1916      M=pg[max],1,pg[max+1],1,pg[max+2],1;
1917      N=1,0,1,0,2,1,2,1,1;
1918    }
1919    // Case 7:
1920    else
1921    {
1922      N=2,1,1,1,0,1,0,1,1;
1923    }
1924  }
1925  // Case 8:
1926  if (sbd==intvec(2,2,2,3))
1927  {
1928    // the orientation of the polygon is important
1929    A=pg[max]-pg[max+1],pg[max+2]-pg[max+1];
1930    if (det(A)==2)
1931    {
1932      N=2,0,1,0,0,1,0,1,1;
1933    }
1934    else
1935    {
1936      N=0,0,1,2,0,1,1,2,1;
1937    }
1938  }
1939  // Case 9:
1940  if (sbd==intvec(2,2,3,3))
1941  {
1942    // if max==1, then the 5th entry in tboundary is the same as the first
1943    if (max==1)
1944    {
1945      max=5;
1946    }
1947    // if boundary=3,2,2,3 then set max=4
1948    if (tboundary[max+1]!=3)
1949    {
1950      max=4;
1951    }
1952    M=pg[max],1,pg[max+1],1,pg[max+2],1;
1953    // the orientation of the polygon matters
1954    A=pg[max-1]-pg[max],pg[max+1]-pg[max];
1955    if (det(A)==4)
1956    {
1957      N=2,0,1,0,0,1,0,2,1;
1958    }
1959    else
1960    {
1961      N=0,2,1,0,0,1,2,0,1;
1962    }
1963  }
1964  // Case 10:
1965  if (sbd==intvec(2,2,3,4))
1966  {
1967    // the orientation of the polygon matters
1968    if (tboundary[max+1]==3)
1969    {
1970      N=3,0,1,0,0,1,0,2,1;
1971    }
1972    else
1973    {
1974      N=0,0,1,3,0,1,2,1,1;
1975    }
1976  }
1977  // Case 11:
1978  if (sbd==intvec(2,3,3,4))
1979  {
1980    N=3,0,1,0,0,1,0,2,1;
1981  }
1982  // Case 12:
1983  if (sbd==intvec(3,3,3,3))
1984  {
1985    N=2,0,1,0,0,1,0,2,1;
1986  }
1987  // Case 13:
1988  if (sbd==intvec(2,2,2,2,2))
1989  {
1990    // compute the angles of the polygon vertices
1991    intvec dt;
1992    for (i=1;i<=5;i++)
1993    {
1994      A=pg[i]-pg[i+1],pg[i+2]-pg[i+1];
1995      dt[i]=det(A);
1996    }
1997    dt[6]=dt[1];
1998    // find the vertex to be mapped to (0,1)
1999    max=1;
2000    while ((dt[max]!=2) or (dt[max+1]!=2))
2001    {
2002      max++;
2003    }
2004    M=pg[max],1,pg[max+1],1,pg[max+2],1;
2005    N=0,1,1,1,2,1,2,1,1;
2006  }
2007  // Case 14:
2008  if (sbd==intvec(2,2,2,2,3))
2009  {
2010    N=2,0,1,0,0,1,0,1,1;
2011  }
2012  // Case 15:
2013  if (sbd==intvec(2,2,2,3,3))
2014  {
2015    // find the vertix to be mapped to (2,0)
2016    if (tboundary[max+1]!=3)
2017    {
2018      max=5;
2019      M=pg[max],1,pg[max+1],1,pg[max+2],1;
2020    }
2021    N=2,0,1,0,0,1,0,2,1;
2022  }
2023  // Case 16:
2024  if (sbd==intvec(2,2,2,2,2,2))
2025  {
2026    N=2,0,1,1,0,1,0,1,1;
2027  }
2028  // we have to transpose the matrices M and N
2029  M=transpose(M);
2030  N=transpose(N);
2031  // compute the unimodular affine transformation, which is of the form
2032  // A11 A12 | T1
2033  // A21 A22 | T2
2034  //  0   0  | 1
2035  T=matrixtointmat(N*inverse(M));
2036  // the upper-left 2x2-block is A
2037  A=T[1..2,1..2];
2038  // the upper-right 2x1-block is the translation vector
2039  trans=T[1,3],T[2,3];
2040  // transform now the lattice points of the polygon with respect to A and T
2041  list nf;
2042  for (i=1;i<=size(polygon);i++)
2043  {
2044    intmat V[2][1]=polygon[i];
2045    V=A*V;
2046    nf[i]=intvec(V[1,1]+trans[1],V[2,1]+trans[2]);
2047    kill V;
2048  }
2049  return(list(findOrientedBoundary(nf)[2],A,trans,nf));
2050}
2051example
2052{
2053   "EXAMPLE:";
2054   echo=2;
2055   ring r=0,(x,y),dp;
2056   // the Newton polygon of the following polynomial
2057   //     has precisely one interior point
2058   poly f=x22y11+x19y10+x17y9+x16y9+x12y7+x9y6+x7y5+x2y3;
2059   list polygon=newtonPolytopeLP(f);
2060   // its lattice points are
2061   polygon;
2062   // find its normal form
2063   list nf=ellipticNF(polygon);
2064   // the vertices of the normal form are
2065   nf[1];
2066   // it has been transformed by the unimodular affine transformation A*x+v
2067   // with matrix A
2068   nf[2];
2069   // and translation vector v
2070   nf[3];
2071   // the 3rd lattice point ...
2072   polygon[3];
2073   // ... has been transformed to
2074   nf[4][3];
2075}
2076
2077
2078proc ellipticNFDB (int n,list #)
2079"USAGE:  ellipticNFDB(n[,#]);   n int, # list
2080ASSUME:  n is an intger between 1 and 16
2081PURPOSE: this is a database storing the 16 normal forms of planar polygons with
2082         precisely one interior point up to unimodular affine transformations
2083@*       (see e.g. Bjorn Poonen, Fernando Rodriguez-Villegas: Lattice Polygons and
2084                   the number 12.  Amer. Math. Monthly  107  (2000),  no. 3, 238--250.)
2085RETURN:  list, L such that
2086@*             L[1] : list whose entries are the vertices of the nth normal form
2087@*             L[2] : list whose entries are all the lattice points of the nth normal form
2088@*             L[3] : only present if the optional parameter # is present, and then
2089                      it is a polynomial in the variables (x,y) whose Newton polygon
2090                      is the nth normal form
2091NOTE:    the optional parameter is only allowed if the basering has the variables x and y
2092EXAMPLE: example ellipticNFDB;   shows an example"
2093{
2094  if ((n<1) or (n>16))
2095  {
2096    ERROR("n is not between 1 and 16.");
2097  }
2098  if (size(#)>0)
2099  {
2100    if ((defined(x)==0) or (defined(y)==0))
2101    {
2102      ERROR("The variables x and y are not defined.");
2103    }
2104  }
2105  if ((defined(x)==0) or (defined(y)==0))
2106  {
2107    ring nfring=0,(x,y),dp;
2108  }
2109  // store the normal forms as polynomials
2110  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,
2111    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,
2112    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;
2113  list pg=newtonPolytopeLP(nf[n]);
2114  if (size(#)==0)
2115  {
2116    return(list(findOrientedBoundary(pg)[2],pg));
2117  }
2118  else
2119  {
2120    return(list(findOrientedBoundary(pg)[2],pg,nf[n]));
2121  }
2122}
2123example
2124{
2125   "EXAMPLE:";
2126   echo=2;
2127   list nf=ellipticNFDB(5);
2128   // the vertices of the 5th normal form are
2129   nf[1];
2130   // its lattice points are
2131   nf[2];
2132}
2133
2134
2135/////////////////////////////////////////////////////////////////////////////////
2136/// AUXILARY PROCEDURES
2137/////////////////////////////////////////////////////////////////////////////////
2138
2139proc polymakeKeepTmpFiles (int i)
2140"USAGE:  polymakeKeepTmpFiles(int i);   i int
2141PURPOSE: some procedures create files in the directory /tmp which are used for
2142         computations with polymake respectively topcom; these will be removed
2143         when the corresponding procedure is left; however, it might be desireable
2144         to keep them for further computations with either polymake or topcom; this
2145         can be achieved by this procedure; call the procedure as:
2146@*       - polymakeKeepTmpFiles(1);  - then the files will be kept
2147@*       - polymakeKeepTmpFiles(0);  - then the files will be removed in the future
2148RETURN:  none"
2149{
2150  if ((i==1) and (defined(polymakekeeptmpfiles)==0))
2151  {
2152    int polymakekeeptmpfiles;
2153    export(polymakekeeptmpfiles);
2154  }
2155  if (i!=1)
2156  {
2157    if (defined(polymakekeeptmpfiles))
2158    {
2159      kill polymakekeeptmpfiles;
2160    }
2161  }
2162}
2163
2164
2165/////////////////////////////////////////////////////////////////////////////////
2166/////////////////////////////////////////////////////////////////////////////////
2167/// AUXILARY PROCEDURES, WHICH ARE DECLARED STATIC
2168/////////////////////////////////////////////////////////////////////////////////
2169/////////////////////////////////////////////////////////////////////////////////
2170/// - scalarproduct
2171/// - intmatcoldelete
2172/// - intmatconcat
2173/// - sortlist
2174/// - minInList
2175/// - stringdelete
2176/// - abs
2177/// - commondenominator
2178/// - maxPosInIntvec
2179/// - maxPosInIntmat
2180/// - sortintvec
2181/// - matrixtointmat
2182/////////////////////////////////////////////////////////////////////////////////
2183
2184static proc scalarproduct (intvec w,intvec v)
2185"USAGE:      scalarproduct(w,v); w,v intvec
2186ASSUME:      w and v are integer vectors of the same length
2187RETURN:      int, the scalarproduct of v and w
2188NOTE:        the procedure is called by findOrientedBoundary"
2189{
2190  int sp;
2191  for (int i=1;i<=size(w);i++)
2192  {
2193    sp=sp+v[i]*w[i];
2194  }
2195  return(sp);
2196}
2197
2198static proc intmatcoldelete (intmat w,int i)
2199"USAGE:      intmatcoldelete(w,i); w intmat, i int
2200RETURN:      intmat, the integer matrix w with the ith comlumn deleted
2201NOTE:        the procedure is called by intmatsort and normalFan"
2202{
2203  if ((i<1) or (i>ncols(w)) or (ncols(w)==1))
2204  {
2205    return(w);
2206  }
2207  if (i==1)
2208  {
2209    intmat M[nrows(w)][ncols(w)-1]=w[1..nrows(w),2..ncols(w)];
2210    return(M);
2211  }
2212  if (i==ncols(w))
2213  {
2214    intmat M[nrows(w)][ncols(w)-1]=w[1..nrows(w),1..ncols(w)-1];
2215    return(M);
2216  }
2217  else
2218  {
2219    intmat M[nrows(w)][i-1]=w[1..nrows(w),1..i-1];
2220    intmat N[nrows(w)][ncols(w)-i]=w[1..nrows(w),i+1..ncols(w)];
2221    return(intmatconcat(M,N));
2222  }
2223}
2224
2225static proc intmatconcat (intmat M,intmat N)
2226"USAGE:      intmatconcat(M,N); M,N intmat
2227RETURN:      intmat, M and N concatenated
2228NOTE:        the procedure is called by intmatcoldelete and sortintmat"
2229{
2230  if (nrows(M)>=nrows(N))
2231  {
2232    int m=nrows(M);
2233
2234  }
2235  else
2236  {
2237    int m=nrows(N);
2238  }
2239  intmat P[m][ncols(M)+ncols(N)];
2240  P[1..nrows(M),1..ncols(M)]=M[1..nrows(M),1..ncols(M)];
2241  P[1..nrows(N),ncols(M)+1..ncols(M)+ncols(N)]=N[1..nrows(N),1..ncols(N)];
2242  return(P);
2243}
2244
2245static proc sortlist (list v,int pos)
2246"USAGE:      sortlist(v,pos); v list, pos int
2247RETURN:      list, the list L ordered in an ascending way according to the pos-th entries
2248NOTE:        called by tropicalCurve"
2249{
2250  if(size(v)==1)
2251  {
2252    return(v);
2253  }
2254  list w=minInList(v,pos);
2255  v=delete(v,w[2]);
2256  v=sortlist(v,pos);
2257  v=list(w[1])+v;
2258  return(v);
2259}
2260
2261static proc minInList (list v,int pos)
2262"USAGE:      minInList(v,pos); v list, pos int
2263RETURN:      list, (v[i],i) such that v[i][pos] is minimal
2264NOTE:        called by sortlist"
2265{
2266  int min=v[1][pos];
2267  int minpos=1;
2268  for (int i=2;i<=size(v);i++)
2269  {
2270    if (v[i][pos]<min)
2271    {
2272      min=v[i][pos];
2273      minpos=i;
2274    }
2275  }
2276  return(list(v[minpos],minpos));
2277}
2278
2279static proc stringdelete (string w,int i)
2280"USAGE:      stringdelete(w,i); w string, i int
2281RETURN:      string, the string w with the ith component deleted
2282NOTE:        the procedure is called by texnumber and choosegfanvector"
2283{
2284  if ((i>size(w)) or (i<=0))
2285  {
2286    return(w);
2287  }
2288  if ((size(w)==1) and (i==1))
2289  {
2290    return("");
2291
2292  }
2293  if (i==1)
2294  {
2295    return(w[2..size(w)]);
2296  }
2297  if (i==size(w))
2298  {
2299    return(w[1..size(w)-1]);
2300  }
2301  else
2302  {
2303    string erg=w[1..i-1],w[i+1..size(w)];
2304    return(erg);
2305  }
2306}
2307
2308static proc abs (def n)
2309"USAGE:      abs(n); n poly or int
2310RETURN:      poly or int, the absolute value of n"
2311{
2312  if (n>=0)
2313  {
2314    return(n);
2315  }
2316  else
2317  {
2318    return(-n);
2319  }
2320}
2321
2322static proc commondenominator (matrix M)
2323"USAGE:   commondenominator(M);  M matrix
2324ASSUME:   the base ring has characteristic zero
2325RETURN:   int, the lowest common multiple of the denominators of the leading coefficients
2326               of the entries in M
2327NOTE:        the procedure is called from polymakeToIntmat"
2328{
2329  int i,j;
2330  int kgV=1;
2331  // successively build the lowest common multiple of the denominators of the leading coefficients
2332  // of the entries in M
2333  for (i=1;i<=nrows(M);i++)
2334  {
2335    for (j=1;j<=ncols(M);j++)
2336    {
2337      kgV=lcm(kgV,int(denominator(leadcoef(M[i,j]))));
2338    }
2339  }
2340  return(kgV);
2341}
2342
2343static proc maxPosInIntvec (intvec v)
2344"USAGE:      maxPosInIntvec(v); v intvec
2345RETURN:      int, the first position of a maximal entry in v
2346NOTE:        called by sortintmat"
2347{
2348  int max=v[1];
2349  int maxpos=1;
2350  for (int i=2;i<=size(v);i++)
2351  {
2352    if (v[i]>max)
2353    {
2354      max=v[i];
2355      maxpos=i;
2356    }
2357  }
2358  return(maxpos);
2359}
2360
2361static proc maxPosInIntmat (intmat v)
2362"USAGE:      maxPosInIntmat(v); v intmat
2363ASSUME:      v has a unique maximal entry
2364RETURN:      intvec, the position (i,j) of the maximal entry in v
2365NOTE:        called by findOrientedBoundary"
2366{
2367  int max=v[1,1];
2368  intvec maxpos=1,1;
2369  int i,j;
2370  for (i=1;i<=nrows(v);i++)
2371  {
2372    for (j=1;j<=ncols(v);j++)
2373    {
2374      if (v[i,j]>max)
2375      {
2376        max=v[i,j];
2377        maxpos=i,j;
2378      }
2379    }
2380  }
2381  return(maxpos);
2382}
2383
2384static proc sortintvec (intvec w)
2385"USAGE:      sortintvec(v); v intvec
2386RETURN:      intvec, the entries of v are ordered in an ascending way
2387NOTE:        called from ellipticNF"
2388{
2389  int j,k,stop;
2390  intvec v=w[1];
2391  for (j=2;j<=size(w);j++)
2392  {
2393    k=1;
2394    stop=0;
2395    while ((k<=size(v)) and (stop==0))
2396    {
2397      if (v[k]<w[j])
2398      {
2399        k++;
2400      }
2401      else
2402      {
2403        stop=1;
2404      }
2405    }
2406    if (k==size(v)+1)
2407    {
2408      v=v,w[j];
2409    }
2410    else
2411    {
2412      if (k==1)
2413      {
2414        v=w[j],v;
2415      }
2416      else
2417      {
2418        v=v[1..k-1],w[j],v[k..size(v)];
2419      }
2420    }
2421  }
2422  return(v);
2423}
2424
2425static proc sortlistbyintvec (list L,intvec w)
2426"USAGE:      sortlistbyintvec(L,w); L list, w intvec
2427RETURN:      list, the entries of L are ordered such that the corresponding reordering of
2428                   w would order w in an ascending way
2429NOTE:        called from ellipticNF"
2430{
2431  int j,k,stop;
2432  intvec v=w[1];
2433  list LL=L[1];
2434  for (j=2;j<=size(w);j++)
2435  {
2436    k=1;
2437    stop=0;
2438    while ((k<=size(v)) and (stop==0))
2439    {
2440      if (v[k]<w[j])
2441      {
2442        k++;
2443      }
2444      else
2445      {
2446        stop=1;
2447      }
2448    }
2449    if (k==size(v)+1)
2450    {
2451      v=v,w[j];
2452      LL=insert(LL,L[j],size(LL));
2453    }
2454    else
2455    {
2456      if (k==1)
2457      {
2458        v=w[j],v;
2459        LL=insert(LL,L[j]);
2460      }
2461      else
2462      {
2463        v=v[1..k-1],w[j],v[k..size(v)];
2464        LL=insert(LL,L[j],k-1);
2465      }
2466    }
2467  }
2468  return(LL);
2469}
2470
2471static proc matrixtointmat (matrix MM)
2472"USAGE:      matrixtointmat(v); MM matrix
2473ASSUME:      MM is a matrix with only integers as entries
2474RETURN:      intmat, the matrix MM has been transformed to type intmat
2475NOTE:        called from ellipticNF"
2476{
2477  intmat M[nrows(MM)][ncols(MM)]=M;
2478  int i,j;
2479  for (i=1;i<=nrows(M);i++)
2480  {
2481    for (j=1;j<=ncols(M);j++)
2482    {
2483      execute("M["+string(i)+","+string(j)+"]="+string(MM[i,j])+";");
2484    }
2485  }
2486  return(M);
2487}
2488
2489//////////////////////////////////////////////////////////////////////////////
2490
2491static proc polygonToCoordinates (list points)
2492"USAGE:      polygonToCoordinates(points);   points list
2493ASSUME:      points is a list of integer vectors each of size two describing the
2494             marked points of a convex lattice polygon like the output of polygonDB
2495RETURN:      list, the first entry is a string representing the coordinates corresponding
2496                   to the latticpoints seperated by commata
2497                   the second entry is a list where the ith entry is a string representing
2498                   the coordinate of corresponding to the ith lattice point
2499                   the third entry is the latex format of the first entry
2500NOTE:        the procedure is called by fan"
2501{
2502  string coord;
2503  list coords;
2504  string latex;
2505  for (int i=1;i<=size(points);i++)
2506  {
2507    coords[i]="u"+string(points[i][1])+string(points[i][2]);
2508    coord=coord+coords[i]+",";
2509    latex=latex+"u_{"+string(points[i][1])+string(points[i][2])+"},";
2510  }
2511  coord=coord[1,size(coord)-1];
2512  latex=latex[1,size(latex)-1];
2513  return(list(coord,coords,latex));
2514}
2515
2516
Note: See TracBrowser for help on using the repository browser.