source: git/Singular/LIB/gfan.lib @ 61fbaf

///////////////////////////////////////////////////////////////////////////////
version="version gfan.lib 4.1.2.0 Feb_2019 "; // $Id$
category = "Convex Geometry";
info="
LIBRARY:  gfan.lib   Interface to gfan and gfanlib for computations in convex geometry
AUTHORS:  Anders N. Jensen, email: jensen@imf.au.dk
          Yue Ren,          email: ren@mathematik.uni-kl.de
          Frank Seelisch

PROCEDURES:
       fullSpace(n);       cone, the ambient space of dimension n
       origin(n);          cone, the origin in an ambient space of dimension n
       positiveOrthant(n); cone, the positive orthant of dimension n
       ambientDimension(c); the dimension of the ambient space the input lives in
       canonicalizeCone(c); a unique representation of the cone c
       codimension(c);      the codimension of the input
       coneViaPoints();     define a cone
       coneViaInequalities(); define a cone
       coneLink(c,w);       the link of c around w
       containsAsFace(c,d); is d a face of c
       containsInSupport(c,d); is d contained in c
       containsPositiveVector(c); contains a vector with only positive entries?
       containsRelatively(c,p); p in c?
       convexHull(c1,c2);   convex hull
       convexIntersection(c1,c2); convex hull
       dimension(c);         dimension of c
       dualCone(c);          the dual of c
       equations(c);         defining equations of c
       faceContaining(c,w);  the face of c containing w in its relative interior
       facets(c);            the facets of c
       generatorsOfLinealitySpace(c); generators of the lineality space of c
       generatorsOfSpan(c);  generators of the span of c
       getLinearForms(c);    linear forms previously stored in c
       getMultiplicity(c);   multiplicity previously stored in c
       inequalities(c);      inequalities of c
       isFullSpace(c);       is the entire ambient space?
       isOrigin(c);          is the origin?
       isSimplicial(c);      is simplicial?
       linealityDimension(c); the dimension of the lineality space of c
       linealitySpace(c);    the lineality space of c
       negatedCone(c);       the negative of c
       polytopeViaInequalities();
       polytopeViaPoints();
       quotientLatticeBasis(c);  basis of Z^n intersected with the span of c modulo Z^n intersected with the lineality space of c
       randomPoint(c);       a random point in the relative interior of c
       rays(c);              generators of the rays of c
       relativeInteriorPoint(c); point in the relative interior of c
       semigroupGenerator(c); generator of Z^n intersected with c modulo Z^n intersected with the lineality space of c
       setLinearForms(c);    stores linear forms in c
       setMultiplicity(c);   stores a multiplicity in c
       span(c);              unique irredundant equations of c
       uniquePoint(c);       a unique point in c stable under reflections at coordinate hyperplanes
       containsInCollection(f,c); f contains c?
       emptyFan(n);           empty fan in ambient dimension n
       fanViaCones(L);        fan generated by the cones in L
       fullFan(n);            full fan in ambient dimension n
       fVector(f);            the f-Vector of f
       getCone(f,d,i[,m]);  the i-th cone of dimension d in f
       insertCone(f,c[,b]);   inserts the cone c into f
       isCompatible(f,c);     f and c live in the same ambient space
       isPure(f);             all maximal cones of f are of the same dimension
       nmaxcones(f);          number of maximal cones in f
       ncones(f);             number of cones in f
       numberOfConesOfDimension(f,d[,m]); the number of cones in dimension d
       removeCone(f,c[,b]);   removes the cone c
       dualPolytope(p);       the dual of p
       newtonPolytope(f);     convex hull of all exponent vectors of f
       vertices(p);           vertices of p
       onesVector(n);         intvec of length n with all entries 1
";

///////////////////////////////////////////////////////////////////////////////

/////
// non gfanlib.so functions
/////

proc fullSpace(int n)
"USAGE:   fullSpace(n);  n int
RETURN:  cone, the ambient space of dimension n
EXAMPLE: example positiveOrthant; shows an example
"
{
  cone c = n;
  return (c);
}
example
{
  "EXAMPLE:"; echo = 2;
  cone c = fullSpace(2);
  generatorsOfLinealitySpace(c);
}

proc origin(int n)
"USAGE:   origin(n);  n int
RETURN:  cone, the origin in an ambient space of dimension n
EXAMPLE: example origin; shows an example
"
{
  bigintmat ineq[0][n];
  bigintmat eq[n][n];
  for (int i=1; i<=n; i++)
  {
    eq[i,i]=1;
  }
  cone o = coneViaInequalities(ineq,eq);
  return (o);
}
example
{
  "EXAMPLE:"; echo = 2;
  cone c = origin(2);
  equations(c);
}

proc positiveOrthant(int n)
"USAGE:   positiveOrthant(n);  n int
RETURN:  cone, the positive orthant of dimension n
EXAMPLE: example positiveOrthant; shows an example
"
{
  bigintmat ineq[n][n];
  for (int i=1; i<=n; i++)
  {
    ineq[i,i]=1;
  }
  cone posOrthant = coneViaInequalities(ineq);
  return (posOrthant);
}
example
{
  "EXAMPLE:"; echo = 2;
  cone c = positiveOrthant(2);
  rays(c);
}

/////
// gfan interface functions
/////

static proc intmatToGfanVectorConfiguration(intmat P)
{
  string gfanVectorConfiguration = "{";
  int c = ncols(P);
  for (int i=1; i<=nrows(P); i++)
  {
    gfanVectorConfiguration = gfanVectorConfiguration
      +"("+string(intvec(P[i,1..c]))+"),";
  }
  int k = size(gfanVectorConfiguration);
  gfanVectorConfiguration = gfanVectorConfiguration[1..k-1];
  gfanVectorConfiguration = gfanVectorConfiguration+"}";
  return (gfanVectorConfiguration);
}

proc secondaryFan(intmat P, def #)
"USAGE:   secondaryFan(P);   P intmat
         secondaryFan(P,s); P intmat, s string
RETURN:  fan, the secondary fan of the vector configuration P
NOTE:    s is a option string that is passed to gfan, possible options are
           '--log1' to '--log3' for output during the computation
EXAMPLE: example secondaryFan; shows an examplex
"
{
  string filename = "/tmp/gfanlib_secondaryFan_"+string(random(1,2147483647));
  string filenameIn = filename+".in.gfan";
  string filenameOut = filename+".out.gfan";

  string filestring = intmatToGfanVectorConfiguration(P);
  filestring;
  write(":w "+filenameIn,filestring);
  int dummy = system("sh","gfan_secondaryfan "+string(#)+" < "+filenameIn+" > "+filenameOut);
  string fanString = read(filenameOut);

  dummy = system("sh","rm "+filenameIn+" "+filenameOut);
  return (fanFromString(fanString));
}
example
{
  "EXAMPLE:"; echo = 2;
  intmat P[4][3] =
    1,0,0,
    1,1,0,
    1,0,1,
    1,1,1;
  fan F = secondaryFan(P);
  F;
}

/////
// cone related functions
/////

proc ambientDimension()
"USAGE:   ambientDimension(c);  c cone
         ambientDimension(f);  f fan
         ambientDimension(p);  p polytope
RETURN:  int, the dimension of the ambient space the input lives in
EXAMPLE: example ambientDimension; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M1[2][2]=
    1,0,
    0,1;
  cone c1=coneViaPoints(M1);
  ambientDimension(c1);
  intmat M2[2][3]=
    1,0,0,
    0,1,0;
  cone c2=coneViaPoints(M2);
  ambientDimension(c2);

  fan f = emptyFan(3);
  ambientDimension(f);
}

proc canonicalizeCone()
"USAGE:   canonicalizeCone(c);  c cone
RETURN:  cone, a unique representation of the cone c
EXAMPLE: example canonicalizeCone; shows an example
"
{
}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[5][3]=
    8,1,9,
    9,2,4,
    0,6,2,
    8,8,8,
    0,9,5;
  cone c=coneViaInequalities(M);
  c;
  cone cc=canonicalizeCone(c);
  cc;
  // computes a unique representation of c
  c == cc;
  // some procedures work with the known inequalities and equations
  // in order to obtain a unique output,
  // bring the cone in canonical form beforehand
  relativeInteriorPoint(c);
  relativeInteriorPoint(cc);
}

proc codimension()
"USAGE:   codimension(c);  c cone
         codimension(f);  f fan
         codimension(p);  p polytope
RETURN:  int, the codimension of the input
EXAMPLE: example codimension; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M1[1][2]=
    1,0;
  cone c1=coneViaPoints(M1);
  codimension(c1);
  intmat M2[2][2]=
    1,0,
    0,1;
  cone c2=coneViaPoints(M2);
  codimension(c2);

  fan f = emptyFan(2);
  codimension(f);
  insertCone(f,c1);
  codimension(f);
  insertCone(f,c2);
  codimension(f);
}

proc coneViaPoints()
"USAGE:   coneViaPoints(HL); intmat HL
          coneViaPoints(HL,L); intmat HL, intmat L
          coneViaPoints(HL,L,flags); intmat HL, intmat L, int flags
RETURN:  cone
PURPOSE: cone generated by half lines generated by the row vectors of HL
  and (if stated) by lines generated by the row vectors of L;
  flags may range between 0,..,3 defining an upper and lower bit
  (0=0*2+0, 1=0*2+1, 2=1*2+0, 3=1*2+1),
  if upper bit is 1, then program assumes that each row vector in HL
  generates a ray of the cone,
  if lower bit is 1, then program assumes that the span of the row
  vectors of L is the lineality space of the cone,
  if either bit is 0, then program computes the information itself.
EXAMPLE: example coneViaPoints; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  // Let's define a cone in R^3 generated by the following half lines:
  intmat HL[5][3]=
   1,0, 0,
  -1,0, 0,
   0,1, 1,
   0,1,-1,
   0,0, 1;
  cone c=coneViaPoints(HL);
  c;
  kill HL,c;
  // Note that (1,0,0) and (-1,0,0) form a line, hence also possible:
  intmat HL[3][3]=
  0,1, 1,
  0,1,-1,
  0,0, 1;
  intmat L[1][3]=
  1,0,0;
  cone c=coneViaPoints(HL,L);
  c;
  kill HL,L,c;
  // lineality space is exactly Lin(1,0,0)
  intmat HL[3][3]=
  0,1, 1,
  0,1,-1,
  0,0, 1;
  intmat L[1][3]=
  1,0,0;
  cone c=coneViaPoints(HL,L,1);
  c;
  kill HL,L,c;
  // and that (0,1,-1), (0,1,1) generate rays
  intmat HL[3][3]=
  0,1, 1,
  0,1,-1;
  intmat L[1][3]=
  1,0,0;
  cone c=coneViaPoints(HL,L,1);
  c;
  kill HL,L,c;
  // and that (0,1,-1), (0,1,1) generate rays
  intmat HL[3][3]=
  0,1, 1,
  0,1,-1;
  intmat L[1][3]=
  1,0,0;
  cone c=coneViaPoints(HL,L,3);
  c;
}

proc coneViaInequalities()
"USAGE:   coneViaInequalities(IE); intmat IE
          coneViaInequalities(IE,E); intmat IE, intmat E
          coneViaInequalities(IE,E,flags); intmat IE, intmat E, int flags
RETURN:  cone
PURPOSE: cone consisting of all points x, such that IE*x >= 0 in each component
  and (if stated) E*x = 0;
  inequalities and (if stated) equations will be transformed, getting rid of
  redundancies;
  flags may range between 0,..,3 defining an upper and lower bit
  (0=0*2+0, 1=0*2+1, 2=1*2+0, 3=1*2+1),
  if higher bit is 1, then program assumes each inequality yields a facet,
  if lower bit is 1, then program assumes the kernel of E is the span of the cone,
  if either bit is 0, then program computes the information itself.
EXAMPLE: example coneViaInequalities; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  // Let's define a cone in R^3 given by the following inequalities:
  intmat IE[6][3]=
  1,3,5,
  1,5,3,
  0,1,-1,
  0,1,1,
  1,0,0,
  -1,0,0;
  cone c=coneViaInequalities(IE);
  c;
  // Note that the last two inequalities yield x1 = 0, hence also possible:
  intmat IE[4][3]=
  0,1,-1,
  0,1,1;
  intmat E[1][3]=
  1,0,0;
  cone c=coneViaInequalities(IE,E);
  c;
  // each inequalities gives rise to a facet
  intmat IE[2][3]=
  0,1,-1,
  0,1,1;
  intmat E[1][3]=
  1,0,0;
  cone c=coneViaInequalities(IE,E,1);
  c;
  // and the kernel of E is the span of the cone
  intmat IE[2][3]=
  0,1,-1,
  0,1,1;
  intmat E[1][3]=
  1,0,0;
  cone c=coneViaInequalities(IE,E,3);
  c;
}

proc coneLink()
"USAGE:   coneLink(c,w);  c cone, w intvec/bigintmat
RETURN:  cone, the link of c around w
EXAMPLE: example coneLink; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[3][3]=
    1,0,0,
    0,1,0,
    0,0,1;
  cone c=coneViaPoints(M);
  intvec v=1,0,0;
  cone cv=coneLink(c,v);
  rays(cv);
  generatorsOfLinealitySpace(cv);
  intvec w=1,1,1;
  cone cw=coneLink(c,w);
  rays(cw);
  generatorsOfLinealitySpace(cw);
}

proc containsAsFace()
"USAGE:   containsAsFace(c,d);  c cone, d cone
RETURN:  1, if d is a face of c; 0 otherwise
EXAMPLE: example containsAsFace; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[2][2]=
    1,0,
    0,1;
  cone c=coneViaPoints(M);
  intmat N1[1][2]=
    1,1;
  cone d1=coneViaPoints(N1);
  containsInSupport(c,d1);
  containsAsFace(c,d1);
  intmat N2[1][2]=
    0,1;
  cone d2=coneViaPoints(N2);
  containsInSupport(c,d2);
  containsAsFace(c,d2);
}

proc containsRelatively()
"USAGE:   containsRelatively(c,p);  c cone, intvec p
RETURN:  1 iff the given cone contains the given point in its relative interior; 0 otherwise
EXAMPLE: example containsRelatively; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[2][2]=
  1,0,
  0,1;
  cone c=coneViaPoints(M);
  intvec p1=1,1;
  containsRelatively(c,p1);
  intvec p2=0,1;
  containsRelatively(c,p2);
}

proc containsInSupport()
"USAGE:   containsInSupport(c,d);  c cone, d cone
         containsInSupport(c,p);   c cone, p intvec/bigintmat
RETURN:  1, if d resp. p is contained in c; 0 otherwise
EXAMPLE: example containsInSupport; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[2][2]=
    1,0,
    0,1;
  cone c=coneViaPoints(M);
  containsInSupport(c,c);
  intmat N1[2][2]=
    1,1,
    0,1;
  cone d1=coneViaPoints(N1);
  containsInSupport(c,d1);
  intmat N2[2][2]=
    1,1,
    1,-1;
  cone d2=coneViaPoints(N2);
  containsInSupport(c,d2);
  intvec p1=0,1;
  containsInSupport(c,p1);
  intvec p2=1,-1;
  containsInSupport(c,p2);
}

proc containsPositiveVector()
"USAGE:   containsPositiveVector(c);  c cone
RETURN:  1, if c contains a vector with only positive entries in its relative interior
EXAMPLE: example containsPositiveVector; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M1[2][2]=
    1,1,
    1,-1;
  cone c1=coneViaPoints(M1);
  containsPositiveVector(c1);
  intmat M2[2][2]=
    0,1,
    -1,0;
  cone c2=coneViaPoints(M2);
  containsPositiveVector(c2);
}

proc convexHull()
"USAGE:   convexHull(c1,c2);   c1 cone, c2 cone
         convexHull(c1,p1);   c1 cone, p1 polytope
         convexHull(p1,c1);   p1 cone, c1 polytope
         convexHull(p1,p2);   p1 polytope, p2 polytope
RETURN:  cone resp polytope, the convex hull of its two input objects
EXAMPLE: example convexHull; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M1[2][2]=
    1,0,
    0,1;
  cone c1=coneViaPoints(M1);
  intmat M2[2][2]=
    1,1,
    1,-1;
  cone c2=coneViaPoints(M2);
  intmat M3[2][2]=
    1,0,
    0,-1;
  cone c3=coneViaPoints(M3);
  cone c12=convexHull(c1,c2);
  c12;
  print(rays(c12));
  cone c23=convexHull(c2,c3);
  c23;
  print(rays(c23));
  cone c13=convexHull(c1,c3);
  c13;
  print(rays(c13));
}

proc convexIntersection()
"USAGE:   convexIntersection(c1,c2);   c1 cone, c2 cone
         convexIntersection(c1,p1);   c1 cone, p1 polytope
         convexIntersection(p1,c1);   p1 cone, c1 polytope
         convexIntersection(p1,p2);   p1 polytope, p2 polytope
RETURN:  cone resp polytope, the convex hull of its two input objects
EXAMPLE: example convexIntersection; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M1[2][2]=
    1,0,
    0,1;
  cone c1=coneViaPoints(M1);
  intmat M2[2][2]=
    1,1,
    1,-1;
  cone c2=coneViaPoints(M2);
  intmat M3[2][2]=
    1,0,
    0,-1;
  cone c3=coneViaPoints(M3);
  cone c12=convexIntersection(c1,c2);
  c12;
  print(rays(c12));
  cone c23=convexIntersection(c2,c3);
  c23;
  print(rays(c23));
  cone c13=convexIntersection(c1,c3);
  c13;
  print(rays(c13));
}

proc dimension()
"USAGE:   dimension(c);  c cone
         dimension(f);  f fan
         dimension(p);  p polytope
RETURN:  int, the dimension of the input
EXAMPLE: example dimension; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M1[1][2]=
    1,0;
  cone c1=coneViaPoints(M1);
  dimension(c1);
  intmat M2[2][2]=
    1,0,
    0,1;
  cone c2=coneViaPoints(M2);
  dimension(c2);

  fan f = emptyFan(2);
  dimension(f);
  insertCone(f,c1);
  dimension(f);
  insertCone(f,c2);
  dimension(f);
}

proc dualCone()
"USAGE:   dualCone(c);  c cone
RETURN:  cone, the dual of c
EXAMPLE: example dualCone; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M1[2][2]=
    1,0,
    0,1;
  cone c1=coneViaPoints(M1);
  cone d1=dualCone(c1);
  d1;
  print(rays(d1));
  intmat M2[2][2]=
    1,1,
    0,1;
  cone c2=coneViaPoints(M2);
  cone d2=dualCone(c2);
  d2;
  print(rays(d2));
}

proc equations()
"USAGE:   equations(c);  c cone
         equations(p);  p polytope
RETURN:  bigintmat, defining equations of c resp p
NOTE:    neither unique nor complete, unless c resp p in canonical form
EXAMPLE: example equations; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M1[2][2]=
    1,0,
    0,1;
  cone c1=coneViaPoints(M1);
  bigintmat E1=equations(c1);
  print(E1);
  intmat M2[1][2]=
    1,0;
  cone c2=coneViaPoints(M2);
  bigintmat E2=equations(c2);
  print(E2);
}

proc faceContaining()
"USAGE:   faceContaining(c,w);  c cone, w intvec/bigintmat
ASSUME:  containsInSupport(c,w)==1
RETURN:  cone, the face of c containing w in its relative interior
EXAMPLE: example faceContaining; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[2][2]=
    1,0,
    0,1;
  cone c=coneViaPoints(M);
  faceContaining(c,intvec(1,0));
  faceContaining(c,intvec(0,1));
  faceContaining(c,intvec(1,1));
  faceContaining(c,intvec(0,0));
}

proc facets()
"USAGE:   facets(c);  c cone
         facets(p);  p polytope
RETURN:  bigintmat, the facets of c resp p
EXAMPLE: example facets; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M1[2][2]=
    1,0,
    0,1;
  cone c1=coneViaPoints(M1);
  bigintmat F1=facets(c1);
  print(F1);
  intmat M2[2][2]=
    1,1,
    0,-1;
  cone c2=coneViaPoints(M2);
  bigintmat F2=facets(c2);
  print(F2);
}

proc generatorsOfLinealitySpace()
"USAGE:   generatorsOfLinealitySpace(c);  c cone
RETURN:  bigintmat, generators of the lineality space of c
EXAMPLE: example generatorsOfLinealitySpace; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[5][3]=
    1,0,0,
    0,1,0,
    0,0,1,
    -1,0,0,
    0,-1,0;
  cone c=coneViaPoints(M);
  bigintmat L=generatorsOfLinealitySpace(c);
  print(L);
}

proc generatorsOfSpan()
"USAGE:   generatorsOfSpan(c);  c cone
RETURN:  bigintmat, generators of the span of c
EXAMPLE: example generatorsOfSpan; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[3][5]=
    1,0,0,0,0,
    0,1,0,0,0,
    0,0,1,0,0;
  cone c=coneViaPoints(M);
  bigintmat S=generatorsOfSpan(c);
  print(S);
}

proc getLinearForms()
"USAGE:   getLinearForms(c);  c cone
         getLinearForms(p);  p polytope
RETURN:  bigintmat, linear forms previously stored in c resp p
EXAMPLE: example getLinearForms; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[2][3]=
    -1,0,0,
    0,-1,0;
  cone c=coneViaPoints(M);
  getLinearForms(c);
  intvec v=1,1,1;
  setLinearForms(c,v);
  getLinearForms(c);
}

proc getMultiplicity()
"USAGE:   getMultiplicity(c);  c cone
         getMultiplicity(p);  p polytope
RETURN:  bigint, 1 or a multiplicity previously stored in c resp p
EXAMPLE: example getMultiplicity; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[2][3]=
    -1,0,0,
    0,-1,0;
  cone c=coneViaPoints(M);
  getMultiplicity(c);
  setMultiplicity(c,3);
  getMultiplicity(c);
}

proc inequalities()
"USAGE:   inequalities(c);  c cone
         inequalities(p);  p polytope
RETURN:  bigintmat, the inequalities of c resp p
NOTE:    neither unique nor irredundant, unless c resp p in canonical form
EXAMPLE: example inequalities; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M1[2][2]=
    1,0,
    0,1;
  cone c1=coneViaPoints(M1);
  bigintmat I1=inequalities(c1);
  print(I1);
  intmat M2[2][2]=
    1,1,
    0,-1;
  cone c2=coneViaPoints(M2);
  bigintmat I2=inequalities(c2);
  print(I2);
}

proc isFullSpace()
"USAGE:   isFullSpace(c);  c cone
RETURN:  1, if c is the entire ambient space; 0 otherwise
EXAMPLE: example isFullSpace; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  cone c1;
  isFullSpace(c1);
  intmat M2[2][2]=
    1,0,
    0,1;
  cone c2=coneViaPoints(M2);
  isFullSpace(c2);
  intmat M3[4][2]=
    1,0,
    0,1,
    -1,0,
    0,-1;
  cone c3=coneViaPoints(M3);
  isFullSpace(c3);
}

proc isOrigin()
"USAGE:   isOrigin(c);  c cone
RETURN:  1, if c is the origin; 0 otherwise
EXAMPLE: example isOrigin; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  cone c1;
  isOrigin(c1);
  intmat M2[2][2]=
    1,0,
    0,1;
  cone c2=coneViaPoints(M2);
  isOrigin(c2);
  intmat M3[4][2]=
    1,0,
    0,1,
    -1,0,
    0,-1;
  cone c3=coneViaPoints(M3);
  isOrigin(c3);
}

proc isSimplicial()
"USAGE:   isSimplicial(c);  c cone
         isSimplicial(f);  f fan
RETURN:  1, if c resp f is simplicial; 0 otherwise
EXAMPLE: example isSimplicial; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M1[3][3]=
    1,0,0,
    0,1,0,
    0,0,1;
  cone c1=coneViaPoints(M1);
  isSimplicial(c1);
  intmat M2[4][3]=
    1,0,0,
    0,1,0,
    0,0,1,
    1,1,-1;
  cone c2=coneViaPoints(M2);
  isSimplicial(c2);
  /***********************/
  fan f=emptyFan(3);
  isSimplicial(f);
  intmat N1[3][3]=
    1,0,0,
    0,1,0,
    0,0,1;
  cone d1=coneViaPoints(N1);
  insertCone(f,d1);
  isSimplicial(f);
  intmat N2[4][3]=
    1,0,0,
    0,1,0,
    1,0,-1,
    0,1,-1;
  cone d2=coneViaPoints(N2);
  insertCone(f,d2);
  isSimplicial(f);
}

proc linealityDimension()
"USAGE:   linealityDimension(c);  c cone
         linealityDimension(f);  f fan
RETURN:  int, the dimension of the lineality space of c resp f
EXAMPLE: example linealityDimension; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M1[3][3]=
    1,0,0,
    0,1,0,
    0,0,1;
  cone c1=coneViaPoints(M1);
  linealityDimension(c1);
  intmat M2[4][3]=
    1,0,0,
    0,1,0,
    0,0,1,
    -1,0,0;
  cone c2=coneViaPoints(M2);
  linealityDimension(c2);
}

proc linealitySpace()
"USAGE:   linealitySpace(c);  c cone
RETURN:  cone, the lineality space of c
EXAMPLE: example linealitySpace; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M1[3][3]=
    1,0,0,
    0,1,0,
    0,0,1;
  cone c1=coneViaPoints(M1);
  cone l1=linealitySpace(c1);
  l1;
  intmat M2[4][3]=
    1,0,0,
    0,1,0,
    0,0,1,
    -1,0,0;
  cone c2=coneViaPoints(M2);
  cone l2=linealitySpace(c2);
  l2;
}

proc negatedCone()
"USAGE:   negatedCone(c);  c cone
RETURN:  cone, the negative of c
EXAMPLE: example negatedCone; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[2][2]=
    1,0,
    0,1;
  cone c=coneViaPoints(M);
  cone cn=negatedCone(c);
  print(rays(cn));
}

proc quotientLatticeBasis()
"USAGE:   quotientLatticeBasis(c);  c cone
RETURN:  bigintmat, a basis of Z^n intersected with the span of c modulo Z^n intersected with the lineality space of c
EXAMPLE: example quotientLatticeBasis; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[3][2]=
    1,0,
    0,1,
    -1,0;
  cone c=coneViaPoints(M);
  bigintmat Q=quotientLatticeBasis(c);
  print(Q);
}

proc randomPoint()
"USAGE:   randomPoint(c);  c cone
         randomPoint(c,b);  c cone, b int
RETURN:  bigintmat, a random point in the relative interior of c
NOTE:    returns a weighted sum over all its rays
         if b is given and b>0, only chooses weights between 1 and b
EXAMPLE: example randomPoint; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[2][2]=
    1,0,
    0,1;
  cone c=coneViaPoints(M);
  bigintmat Q=randomPoint(c);
  print(Q);
  bigintmat P=randomPoint(c,5);
  print(P);
}

proc rays()
"USAGE:   rays(c);  c cone
RETURN:  bigintmat, generators of the rays of c, orthogonal to its lineality space
EXAMPLE: example rays; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M1[2][2]=
    1,0,
    0,1;
  cone c1=coneViaPoints(M1);
  bigintmat R1=rays(c1);
  print(R1);
  intmat M2[3][2]=
    1,0,
    0,1,
    -1,0;
  cone c2=coneViaPoints(M2);
  bigintmat R2=rays(c2);
  print(R2);
}

proc relativeInteriorPoint()
"USAGE:   relativeInteriorPoint(c);  c cone
RETURN:  bigintmat, a point in the relative interior of c
NOTE:    not unique, unless c is in its canonical form
EXAMPLE: example relativeInteriorPoint; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M1[2][2]=
    1,0,
    0,1;
  cone c1=coneViaPoints(M1);
  relativeInteriorPoint(c1);
  intmat M2[2][2]=
    1,0,
    1,1;
  cone c2=coneViaPoints(M2);
  relativeInteriorPoint(c2);
}

proc semigroupGenerator()
"USAGE:   semigroupGenerator(c);  c cone
RETURN:  bigintmat, the generator of Z^n intersected with c modulo Z^n intersected with the lineality space of c
ASSUME:  dimension(c) == linealityDimension(c)+1
EXAMPLE: example semigroupGenerator; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[3][2]=
    1,0,
    0,1,
    -1,0;
  cone c=coneViaPoints(M);
  semigroupGenerator(c);
}

proc setLinearForms()
"USAGE:   setLinearForms(c);  c cone
         setLinearForms(p);  p polytope
RETURN:  none, stores linear forms in c resp p
EXAMPLE: example setLinearForms; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[2][3]=
    -1,0,0,
    0,-1,0;
  cone c=coneViaPoints(M);
  getLinearForms(c);
  intvec v=1,1,1;
  setLinearForms(c,v);
  getLinearForms(c);
}

proc setMultiplicity()
"USAGE:   setMultiplicity(c);  c cone
         setMultiplicity(p);  p polytope
RETURN:  none, stores a multiplicity in c resp p
EXAMPLE: example setMultiplicity; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[2][3]=
    -1,0,0,
    0,-1,0;
  cone c=coneViaPoints(M);
  getMultiplicity(c);
  setMultiplicity(c,3);
  getMultiplicity(c);
}

proc span()
"USAGE:   span(c);  c cone
RETURN:  bigintmat, unique irredundant equations of c
NOTE:    the name 'span' was chosen to be in line with polymake's nomenclature
EXAMPLE: example span; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[3][5]=
    1,0,0,0,0,
    0,1,0,0,0,
    0,0,1,0,0;
  cone c=coneViaPoints(M);
  bigintmat Eq=span(c);
  print(Eq);
}


proc uniquePoint()
"USAGE:   uniquePoint(c);  c cone
RETURN:  bigintmat, a unique point in c stable under reflections at coordinate hyperplanes
EXAMPLE: example uniquePoint; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M1[2][2]=
    1,0,
    0,1;
  cone c1=coneViaPoints(M1);
  uniquePoint(c1);
  intmat M2[2][2]=
    -1,0,
    0,1;
  cone c2=coneViaPoints(M2);
  uniquePoint(c2);
}



/////
// fan related functions
/////

proc containsInCollection()
"USAGE:   containsInCollection(f,c);  f fan, c cone
RETURN:  1, if f contains c; 0 otherwise
EXAMPLE: example containsInCollection; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  fan f=emptyFan(2);
  intmat M[2][2]=
    1,0,
    0,1;
  cone c=coneViaPoints(M);
  containsInCollection(f,c);
  insertCone(f,c);
  containsInCollection(f,c);
}

proc emptyFan()
"USAGE:   emptyFan(n);  n int
RETURN:  fan, an empty fan in ambient dimension n
EXAMPLE: example emptyFan; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  fan f=emptyFan(2);
  f;
}

proc fanViaCones()
"USAGE:   fanViaCones(L);  L list
         fanViaCones(c1[,...,ck]);  c1,...,ck cones
RETURN:  fan, creates a fan generated by the cones in L resp c1,...,ck
EXAMPLE: example fanViaCones; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[2][2]=1,0,0,1;
  cone c=coneViaPoints(M);
  intmat N[2][2]=1,0,0,-1;
  cone d=coneViaPoints(N);
  fan f=fanViaCones(c,d);
  f;
  list L=c,d;
  fan g=fanViaCones(L);
  g;
}

proc fullFan()
"USAGE:   fullFan(n);  n int
RETURN:  fan, an full fan in ambient dimension n
EXAMPLE: example fullFan; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  fan f=fullFan(2);
  f;
}

proc fVector()
"USAGE:   fVector(f);  f fan
RETURN:  bigintmat, the f-Vector of f
EXAMPLE: example fVector; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  fan f=emptyFan(2);
  fVector(f);
  intmat M[2][2]=1,0,0,1;
  cone c=coneViaPoints(M);
  insertCone(f,c);
  fVector(f);
}

proc getCone()
"USAGE:   getCone(f,d,i[,m]);  f fan, d int, i int, m int
ASSUME:  d is between 0 and ambientDimension(f)
         i is between 1 and numberOfConesOfDimension(f,d,o,m)
RETURN:  cone, returns in the fan f of all cones in dimension d the i-th cone
         if m!=0, it will enumerate over maximal cones only
EXAMPLE: example getCone; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[3][3]=
    1,0,0,
    0,1,0,
    0,0,1;
  cone c=coneViaPoints(M);
  fan f=emptyFan(3);
  insertCone(f,c);
  getCone(f,2,1,0);
  getCone(f,2,2,0);
}

proc insertCone()
"USAGE:   insertCone(f,c[,b]);  f fan, c cone, b int
ASSUME:  isCompatible(f,c)=1
RETURN:  none, inserts the cone c into f
         if b=0, then skips check whether f and c are compatible
EXAMPLE: example insertCone; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  fan f=emptyFan(3);
  f;
  intmat M[3][3]=
    1,0,0,
    0,1,0,
    0,0,1;
  cone c=coneViaPoints(M);
  insertCone(f,c);
  f;
}

proc isCompatible()
"USAGE:   isCompatible(f,c);  f fan, c cone
RETURN:  1 if f and c live in the same ambient space and
           if the intersection of c with any cone of f is a face of each;
         0 otherwise
EXAMPLE: example isCompatible; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  fan f=emptyFan(3);
  intmat M1[3][3]=
    1,0,0,
    0,1,0,
    0,0,1;
  cone c1=coneViaPoints(M1);
  isCompatible(f,c1);
  insertCone(f,c1);
  intmat M2[3][3]=
    1,1,1,
    1,0,0,
    0,1,0;
  cone c2=coneViaPoints(M2);
  isCompatible(f,c2);
  intmat M3[3][3]=
    1,0,0,
    0,1,0,
    0,0,-1;
  cone c3=coneViaPoints(M3);
  isCompatible(f,c3);
}

proc isPure()
"USAGE:   isPure(f);  f fan
RETURN:  1 if all maximal cones of f are of the same dimension
         0 otherwise
EXAMPLE: example isPure; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  fan f=fullFan(2);
  isPure(f);
  fan g=emptyFan(2);
  intmat M1[2][2]=
    1,0,
    0,1;
  cone c1=coneViaPoints(M1);
  insertCone(g,c1);
  isPure(g);
  intmat M2[1][2]=
    0,-1;
  cone c2=coneViaPoints(M2);
  insertCone(g,c2);
  isPure(g,c2);
}

proc nmaxcones()
"USAGE:   nmaxcones(f);  f fan
RETURN:  int, the number of maximal cones in f
EXAMPLE: example nmaxcones; shows an example
"
{

}
example
{
  fan f=emptyFan(3);
  nmaxcones(f);
  intmat M1[3][3]=
    1,0,0,
    0,1,0,
    0,0,1;
  cone c1=coneViaPoints(M1);
  insertCone(f,c1);
  nmaxcones(f);
  intmat M2[2][3]=
    1,0,0,
    0,-1,0;
  cone c2=coneViaPoints(M2);
  insertCone(f,c2);
  nmaxcones(f);
}

proc ncones()
"USAGE:   ncones(f);  f fan
RETURN:  int, the number of cones in f
EXAMPLE: example ncones; shows an example
"
{

}
example
{
  fan f=emptyFan(3);
  ncones(f);
  intmat M1[3][3]=
    1,0,0,
    0,1,0,
    0,0,1;
  cone c1=coneViaPoints(M1);
  insertCone(f,c1);
  ncones(f);
  intmat M2[2][3]=
    1,0,0,
    0,-1,0;
  cone c2=coneViaPoints(M2);
  insertCone(f,c2);
  ncones(f);
}

proc numberOfConesOfDimension()
"USAGE:   numberOfConesOfDimension(f,d[,m]);  f fan, d int, m int
ASSUME:  d is between 0 and ambientDimension(f)
RETURN:  cone, returns in the fan f the number of cones in dimension d
         if m!=0, it will only count maximal cones
EXAMPLE: example numberOfConesOfDimension; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  fan f=emptyFan(3);
  ncones(f);
  intmat M[3][3]=
    1,0,0,
    0,1,0,
    0,0,1;
  cone c=coneViaPoints(M);
  insertCone(f,c);
  numberOfConesOfDimension(f,0,0);
  numberOfConesOfDimension(f,0,1);
  numberOfConesOfDimension(f,1,0);
  numberOfConesOfDimension(f,0,1);
  numberOfConesOfDimension(f,2,0);
  numberOfConesOfDimension(f,2,1);
  numberOfConesOfDimension(f,3,0);
  numberOfConesOfDimension(f,3,1);
}

proc removeCone()
"USAGE:   removeCone(f,c[,b]);  f fan, c cone, b int
ASSUME:  containsInCollection(f,c)=1
RETURN:  none, removes the cone c from f
         if b=0, skips the check whether c is contained in f
EXAMPLE: example removeCone; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[2][2]=1,0,0,1;
  intmat N[2][2]=1,0,1,-1;
  cone c=coneViaPoints(M);
  cone d=coneViaPoints(N);
  fan f=emptyFan(2);
  insertCone(f,c);
  insertCone(f,d);
  f;
  removeCone(f,c);
  f;
}



/////
// polytope related functions
/////

proc dualPolytope()
"USAGE:   dualPolytope(p);  p polytope
RETURN:  polytope, the dual of p
EXAMPLE: example dualPolytope; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[4][2]=
    0,0,
    1,0,
    0,1,
    1,1;
  polytope p=polytopeViaPoints(M);
  dualPolytope(p);
}


proc newtonPolytope()
"USAGE:   newtonPolytope(f);  f poly
RETURN:  polytope, the convex hull of all exponent vectors of f
EXAMPLE: example newtonPolytope; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  ring r;
  poly f=x+y+z;
  polytope p=newtonPolytope(f);
  p;
}

proc polytopeViaPoints()
"USAGE: polytopeViaPoints(V [, flags]);  intmat V, int flags
RETURN: polytope which is the intersection of the cone generated by the row vectors
of V with the hyperplane, in which the first coordinate equals 1;
flags may be 0 or 1,@*
if flags is 1, then program assumes that each row vector of M generates a ray in the cone,
if flags is 0, then program computes that information itself
EXAMPLE: example polytopeViaPoints; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  // This is a polytope in R^2 generated by (0,0), (1,0), (0,1), (0,0);
  intmat V[4][3]=
  1,0,0,
  1,1,0,
  1,0,1,
  1,1,1;
  polytope p1=polytopeViaPoints(V);
  p1;
  // This is a polytope in R^2 generated by (1/2,2/3), (3/4,4/5), (5/6,6/7):
  intmat V[3][3]=
  6,3,4,
  20,15,16,
  42,35,36;
  polytope p2=polytopeViaPoints(V);
  p2;
  // This polytope is the positive orthant in R^2:
  // (0,1,0) and (0,0,1) imply that the polytope is unbounded in that direction
  intmat V[3][3]=
  1,0,0,
  0,1,0,
  0,0,1;
  polytope p3=polytopeViaPoints(V);
  p3;
}

proc polytopeViaInequalities()
"USAGE: polytopeViaInequalities(EV [, E [, flags]]);  intmat EV,E, int flags
RETURN: polytope consisting of all points x, such that IE*x >= 0 in each component
and (if stated) E*x = 0;
flags may range between 0,..,3 defining an upper and lower bit
(0=0*2+0, 1=0*2+1, 2=1*2+0, 3=1*2+1),
if higher bit is 1, then program assumes each inequality yields a facet,
if lower bit is 1, then program assumes the kernel of E is the span of the cone,
if either bit is 0, then program computes the information itself.
EXAMPLE: example polytopeViaPoints; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat IE[2][3]=
  1,0,0,
  0,1,0;
  intmat E[1][3]=
  0,0,1;
  polytope p=polytopeViaInequalities(IE,E);
  p;
}

proc vertices()
"USAGE:   vertices(p);  p polytope
RETURN:  bigintmat, the vertices of p modulo its lineality space
EXAMPLE: example vertices; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intmat M[4][3]=
    1,0,0,
    1,2,0,
    1,0,2,
    1,2,2,
    1,1,1;
  polytope p=polytopeViaPoints(M);
  vertices(p);
}

proc onesVector()
"USAGE:   onesVector(n);  n int
RETURN:  intvec, intvec of length n with all entries 1
EXAMPLE: example onesVector; shows an example
"
{

}
example
{
  "EXAMPLE:"; echo = 2;
  intvec w = onesVector(3);
  w;
}

static proc mod_init()
{
  intvec save=option(get);
  option(noredefine);
  LIB "customstd.so";
  LIB "gfanlib.so";
  option(set,save);
}