source: git/Singular/LIB/polymake.lib @ 334c21f

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