source: git/Singular/LIB/oldpolymake.lib @ 1b2216

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