source: git/Singular/LIB/polymake.lib @ 8275a9

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