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

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