source: git/Singular/LIB/polymake.lib @ 6391eb

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