singular-data 1:4.1.0-p3+ds-2build1 / usr / share / singular / LIB / tropicalNewton.lib

///////////////////////////////////////////////////////////////////////////////
version="version tropical.lib 4.0.3.3 Nov_2016 "; // $Id: 5cf6b6a19bfc3f72916c42f0d8b9890c334149f5 $
category="Tropical Geometry";
info="
LIBRARY:         tropicalNewton.lib  Computations in Tropical Geometry
                                     using Newton Polygon methods

AUTHORS:         Tommy Hofman,       email: thofmann@mathematik.uni.kl.de
                 Yue Ren,            email: reny@cs.bgu.ac.il

OVERVIEW:
  This libraries contains algorithms for computing
  - non-trivial points on tropical varieties,
  - zero-dimensional tropical varieties,
  - one-codimensional links of tropical varieties
  based on Newton polygon methods.

REFERENCES: Hofmann, Ren: Computing tropical points and tropical links, arXiv:1611.02878
            (WARNING: this library follows the max convention instead
                      and triangular sets follow the definition of the Singular book)

SEE ALSO: tropicalVariety, tropical_lib

KEYWORDS: tropical geometry; tropical varieties; Newton polygons

PROCEDURES:
  setUniformizingParameter()  sets the uniformizingParameter
  val()                       returns valuation of element in ground field
  newtonPolygonNegSlopes()    returns negative of the Newton Polyong slopes
  cccMatrixToPositiveIntvec() helper function to turn a computed valuation vector
                              into a usable weight vector in Singular
  tropicalPointNewton()       computes point on tropical variety
  switchRingsAndComputeInitialIdeal()
                              switches rings and computes initial ideal
  tropicalVarietyNewton()     computes tropical variety of zero-dimensional ideal
  tropicalLinkNewton()        computes tropical variety that is polyhedral fan
                                and has codimension one lineality space
";

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

LIB "gfan.lib";
LIB "triang.lib";
LIB "linalg.lib";

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

proc setUniformizingParameter(number p)
"USAGE:   setUniformizingParameter(p); p number
RETURN:  none, sets the uniformizing parameter as p
ASSUME:  char(K)==0 and p prime number or
         trdeg(K)>0 and p transcendental variable or
         p==0
EXAMPLE: example setUniformizingParameter; shows an example
"
{
  // kill uniformizingParemeter if previously defined
  if (defined(uniformizingParameter))
  {
    kill uniformizingParameter;
  }

  // input sanity check
  if (p!=0)
  {
    if (pardeg(p)==0)
    {
      if(char(basering)!=0 && prime(int(p))!=p)
      {
        ERROR("setUniformizingParameter: unexpected parameters");
      }
    }
    else
    {
      if (size(p)!=2 || pardeg(p)!=1)
      {
        ERROR("setUniformizingParameter: unexpected parameters");
      }
    }
  }

  // set uniformizingParameter and export it
  number uniformizingParameter = p;
  export(uniformizingParameter);
}
example
{ "EXAMPLE:"; echo = 2;
  // poor man's polynomials over Puiseux series:
  ring r = (0,t),x,dp;
  setUniformizingParameter(t);
  val(t2+t3);
  val(t^-2+t^-3);
  // poor man's polynomials over p-adic numbers:
  ring s = 0,x,dp;
  setUniformizingParameter(2);
  val(12);
  val(1/12);
}

proc val(number c)
"USAGE:   val(c); c number
RETURN:  int, the valuation of a element in the ground field
ASSUME:  uniformizingParameter is set and c!=0
EXAMPLE: example val; shows an example
"
{
  if (defined(uniformizingParameter)==0)
  {
    ERROR("val: uniformizingParameter not set");
  }
  if (c==0)
  {
    ERROR("val: input is 0");
  }

  if (uniformizingParameter==0)
  {
    return (0);
  }

  int vc = 0;
  if (pardeg(uniformizingParameter)>0)
  {
    def origin = basering;
    number c_denom = denominator(c);
    number c_num = numerator(c);
    execute("ring r_Val=0,"+string(uniformizingParameter)+",ds");
    poly c_denom = imap(origin,c_denom);
    poly c_num = imap(origin,c_num);
    vc = ord(c_num)-ord(c_denom);
    setring origin;
  }
  else
  {
    int c_denom = int(denominator(c));
    int c_num = int(numerator(c));
    int p = int(uniformizingParameter);

    while (c_num mod p==0)
    {
      c_num = c_num div p;
      vc++;
    }
    while (c_denom mod p==0)
    {
      c_denom = c_denom div p;
      vc--;
    }
  }

  return (vc);
}
example
{ "EXAMPLE:"; echo = 2;
  // poor man's polynomials over Puiseux series:
  ring r = (0,t),x,dp;
  setUniformizingParameter(t);
  val(t2+t3);
  val(t^-2+t^-3);
  // poor man's polynomials over p-adic numbers:
  ring s = 0,x,dp;
  setUniformizingParameter(2);
  val(12);
  val(1/12);
}

proc newtonPolygonNegSlopes(poly g, list #)
"USAGE:   newtonPolygonNegSlopes(g,b); g poly, b int
RETURN:  list, the negative slopes of the Newton Polygon of g
         if b==1, computes root (type poly) instead if (easily) possible
ASSUME:  uniformizingParameter is set and g univariate
EXAMPLE: example newtonPolygonNegSlopes; shows an example
"
{
  if (size(#)>0)
  {
    if (#[1]==1 && size(g)==2)
    {
      if (deg(g[1])==1 && deg(g[2])==0)
      {
        poly root = -leadcoef(g[2])/leadcoef(g[1]);
        return (list(root));
      }
    }
  }

  if (defined(uniformizingParameter)==0)
  {
    ERROR("newtonPolygonNegSlopes: uniformizingParameter not set");
  }
  int k = size(g);
  intmat M [k+1][3];
  for (int i=1; i<=k; i++)
  {
    M[i,1] = 1;
    M[i,2] = deg(g[i]);
    M[i,3] = val(leadcoef(g[i]));
  }
  M[k+1,1] = 0;
  M[k+1,2] = 0;
  M[k+1,3] = 1;
  cone Delta = coneViaPoints(M);
  bigintmat F = facets(Delta);

  list slopes;
  number slope;
  for (i=1; i<=nrows(F); i++)
  {
    if (F[i,3]!=0)
    {
      slope = number(F[i,2])/number(F[i,3]);
      slopes[size(slopes)+1] = slope;
    }
  }
  return (slopes);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = (0,t),x,dp;
  setUniformizingParameter(t);
  poly g = tx2+x+1;
  newtonPolygonNegSlopes(g);
  // poor man's polynomials over p-adic numbers:
  ring s = 0,x,dp;
  setUniformizingParameter(3);
  poly g = x2+9x+1;
  newtonPolygonNegSlopes(g);
}


static proc listDot(list VV, intvec expv)
{
  number ld;
  for (int i=1; i<=size(expv); i++)
  {
    if (expv[i]>0)
    {
      if (typeof(VV[i])=="number")
      {
        ld = ld + VV[i]*expv[i];
      }
      else
      {
        ERROR("listDot: unkown or invalid relevant entry");
      }
    }
  }
  return (ld);
}


static proc lexSmallestVariableIndex(poly g)
{
  intvec alpha = leadexp(g);
  for (int i=1; i<=nvars(basering); i++)
  {
    if (alpha[i]>0)
    {
      return (i);
    }
  }
  return (-1);
}

static proc numbersLessThanAsRationalNumbers(number a, number b)
{
  int a_num = int(numerator(a));
  int a_denom = int(denominator(a));
  int b_num = int(numerator(b));
  int b_denom = int(denominator(b));
  return (a_num*b_denom-b_num*a_denom<0);
}

static proc expectedValuation(poly h, list VV)
{
  number vh = val(leadcoef(h[1]))+listDot(VV,leadexp(h[1]));
  number vhh;
  int amb = 0;
  for (int i=2; i<=size(h); i++)
  {
    vhh = val(leadcoef(h[i]))+listDot(VV,leadexp(h[i]));
    if (vh == vhh)
    {
      amb = 1;
    }
    if (numbersLessThanAsRationalNumbers(vhh,vh))
    {
      vh = vhh;
      amb = 0;
    }
  }
  if (amb > 0)
  {
    "h:"; h;
    "VV:"; VV;
    ERROR("expectedValuation: valuation ambiguous");
  }

  return (vh);
}

static proc expectedNewtonPolygonNegSlopes(poly g, list VV, list #)
{
  int k = size(VV);
  for (int i=1; i<=size(VV); i++)
  {
    if (typeof(VV[i])=="poly")
    {
      g = subst(g,var(i),VV[i]);
    }
  }

  if (size(#)>0)
  {
    if (#[1]==1 && size(g)==2)
    {
      if (deg(g[1])==1 && deg(g[2])==0)
      {
        poly root = -leadcoef(g[2])/leadcoef(g[1]);
        return (list(root));
      }
    }
  }

  int n = lexSmallestVariableIndex(g);
  matrix G = coef(g,var(n));

  k = ncols(G);
  intmat M[k+1][3];
  number vh;
  for (i=1; i<=k; i++)
  {
    vh = expectedValuation(G[2,i],VV);
    M[i,1] = int(denominator(vh));
    M[i,2] = deg(G[1,i]);
    M[i,3] = int(numerator(vh));
  }
  M[k+1,1] = 0;
  M[k+1,2] = 0;
  M[k+1,3] = 1;

  cone Delta = coneViaPoints(M);
  bigintmat F = facets(Delta);

  list slopes;
  number slope;
  for (i=1; i<=nrows(F); i++)
  {
    if (F[i,3]!=0)
    {
      slope = number(F[i,2])/number(F[i,3]);
      slopes[size(slopes)+1] = slope;
    }
  }
  return (slopes);
}


static proc randomValuedNumber(list #)
{
  int n = 99;
  if (size(#)>0)
  {
    n = #[1];
  }
  int v = random(0,n);

  number c = random(1,99)*uniformizingParameter^v;
  c = c+random(0,99)*uniformizingParameter^(v+random(0,n div 3));
  c = c+random(0,99)*uniformizingParameter^(v+random(0,n div 3));

  return (c);
}


static proc listOfNumbersToMatrix(list L)
{
  int k = size(L);
  matrix Lmat[1][k];
  for (int i=1; i<=k; i++)
  {
    if (typeof(L[i]) == "number")
    {
      Lmat[1,i] = L[i];
    }
    if (typeof(L[i]) == "poly")
    {
      Lmat[1,i] = val(number(L[i]));
    }
  }

  return (Lmat);
}

proc cccMatrixToPositiveIntvec(matrix L)
"USAGE:   cccMatrixToPositiveIntvec(M); M matrix
RETURN:  intvec, strictly positive equivalent as weight vector to row vector in M
ASSUME:  constant coefficient case only,
         will scale weight vector and add vectors of ones to it
EXAMPLE: example cccMatrixToPositiveIntvec; shows an example
"
{
  int k = ncols(L);
  intvec denoms;
  intvec enums;
  int denoms_gcd = int(denominator(number(L[1,1])));
  int denoms_prod = int(denominator(number(L[1,1])));
  for (int i=2; i<=k; i++)
  {
    denoms_gcd = gcd(denoms_gcd,int(denominator(number(L[1,i]))));
    denoms_prod = denoms_prod*int(denominator(number(L[1,i])));
  }
  int denoms_lcm = denoms_prod div denoms_gcd;

  intvec w;
  for (i=1; i<=k; i++)
  {
    w[i] = int(number(L[1,i])*denoms_lcm);
  }
  return (addOneVectorUntilPositive(w));
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = (0,t),(p01,p02,p12,p03,p13,p23,p04,p14,p24,p34),dp;
  number uniformizingParameter = t;
  export(uniformizingParameter);
  ideal I =
    p23*p14-p13*p24+p12*p34,
    p23*p04-p03*p24+p02*p34,
    p13*p04-p03*p14+p01*p34,
    p12*p04-p02*p14+p01*p24,
    p12*p03-p02*p13+p01*p23;
  system("--random",1337);
  matrix p = tropicalPointNewton(I);
  print(p);
  intvec w = cccMatrixToPositiveIntvec(p);
  print(w);
  def s = switchRingsAndComputeInitialIdeal(I,w);
  kill uniformizingParameter;
}

static proc oneVector(int n)
{
  intvec w;
  for (int i=1; i<=n; i++)
  {
    w[i] = 1;
  }
  return (w);
}

static proc addOneVectorUntilPositive(intvec w)
{
  int n = size(w);
  int w_min = w[1];
  for (int i=2; i<=n; i++)
  {
    if (w[i]<w_min)
    {
      w_min = w[i];
    }
  }
  return (w-(w_min-1)*oneVector(n));
}

proc tropicalPointNewton(ideal I, list #)
"USAGE:   tropicalPointsLasVegas(I); I ideal
RETURN:  matrix, a matrix containing a tropical point as row vector
ASSUME:  uniformizingParameter is set and I monomial free
NOTE:    if printlevel sufficiently high will print intermediate data and timings
         returns error if randomly chosen hyperplanes are not generic
EXAMPLE: example tropicalPointNewton; shows an example
"
{
  if (defined(uniformizingParameter)==0)
  {
    ERROR("tropicalPointNewton: uniformizingParameter not set");
  }
  int ttotal = timer;
  ring origin = basering;

  int tindepSet = timer;
  ideal stdI = std(I);
  attrib(stdI,"isSB",1);
  intvec u = indepSet(stdI);
  tindepSet = timer - tindepSet;
  dbprint("maximal independent set: "+string(u));

  int n = nvars(origin);
  string vars;
  string subststr = "ideal J = subst(I,";
  list VV;
  number c;
  for (int i=1; i<=n; i++)
  {
    if (u[i]==0)
    {
      vars = vars+varstr(i)+",";
    }
    else
    {
      c = randomValuedNumber(#);
      VV[i]=number(val(c));
      dbprint("substituting "+string(var(i))+" with number of valuation "+string(VV[i]));
      subststr = subststr+varstr(i)+","+string(c)+",";
    }
  }
  vars = vars[1..size(vars)-1];
  subststr = string(subststr[1..size(subststr)-1])+");";
  execute(subststr);

  execute("ring s = ("+charstr(origin)+"),("+vars+"),lp;");
  ideal J = imap(origin,J);
  int tstdlp = timer;
  dbprint("computing triangular decomposition (picking first factor)");
  J = satstd(J);
  attrib(J,"isSB",1);
  if (dim(J)>0)
  {
    dbprint("not zero-dimensional, please restart algorithm");
    bigintmat O[1][1];
    return (O);
  }
  if (!isTriangularSet(J))
  {
    J = triangL(J)[1];
  }
  tstdlp = timer - tstdlp;

  dbprint("starting analysis of Newton polygons");
  setring origin;
  J = imap(s,J);

  int tnewton = timer;
  list Vlocal = newtonPolygonNegSlopes(J[1],1);
  int l = lexSmallestVariableIndex(J[1]);
  VV[l] = Vlocal[1];
  dbprint("possible valuations for "+string(var(l))+" (picking first): "+string(Vlocal));

  for (i=2; i<=size(J); i++)
  {
    Vlocal = expectedNewtonPolygonNegSlopes(J[i],VV,1);
    l = lexSmallestVariableIndex(J[i]);
    dbprint("possible valuations for "+string(var(l))+" (picking first): "+string(Vlocal));
    VV[l] = Vlocal[1];
  }
  tnewton = timer-tnewton;
  matrix w = -listOfNumbersToMatrix(VV);

  ttotal = timer-ttotal;
  dbprint("time used total: "+string(ttotal));
  dbprint("computing independent set: "+string(tindepSet));
  dbprint("computing triangular decomposition: "+string(tstdlp));
  dbprint("analyzing newton polygons: "+string(tnewton));

  return (w);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = (0,t),(p01,p02,p12,p03,p13,p23,p04,p14,p24,p34),dp;
  number uniformizingParameter = t;
  export(uniformizingParameter);
  ideal I =
    p23*p14-p13*p24+p12*p34,
    p23*p04-p03*p24+p02*p34,
    p13*p04-p03*p14+p01*p34,
    p12*p04-p02*p14+p01*p24,
    p12*p03-p02*p13+p01*p23;
  system("--random",1337);
  printlevel = 3;
  matrix p = tropicalPointNewton(I);
  print(p);
  intvec w = cccMatrixToPositiveIntvec(p);
  print(w);
  def s = switchRingsAndComputeInitialIdeal(I,w);
  kill uniformizingParameter;
}


static proc sumOfLeadExpV(poly f)
{
  intvec expvSum;
  expvSum = leadexp(f[1]);
  for (int i=2; i<=size(f); i++)
  {
    expvSum = expvSum + leadexp(f[i]);
  }
  return (expvSum);
}


static proc isTriangularSet(ideal F)
{
  // checks whether F has the right amount of elements
  int n = size(F);
  if (n!=nvars(basering))
  {
    return (0);
  }
  int i,j;
  poly fi;
  intvec expvSum;
  for (i=1; i<=n; i++)
  {
    fi = F[i];
    expvSum = sumOfLeadExpV(fi);
    // checks whether fi has monomial containing x_i
    if (expvSum[n-i+1]==0)
    {
      return (0);
    }
    // checks whether fi has no monomial containing x_j, j>i
    for (j=i+1; j<=n; j++)
    {
      if (expvSum[n-j+1]>0)
      {
        return (0);
      }
    }
  }
  return (1);
}


proc tropicalVarietyNewton(ideal I)
"USAGE:   tropicalVarietyNewton(I); I ideal
RETURN:  matrix, a matrix containing all elements of the tropical variety
ASSUME:  uniformizingParameter is set, I monomial free and zero-dimensional
EXAMPLE: example tropicalVarietyNewton; shows an example
"
{
  if (defined(uniformizingParameter)==0)
  {
    ERROR("tropicalVarietyNewton: uniformizingParameter not set");
  }
  if(ordstr(basering)[1,2]!="lp")
  {
    ERROR("tropicalVarietyNewton: order must be lp");
  }

  if (!isTriangularSet(I))
  {
    I = satstd(I);
    attrib(I,"isSB",1);
    if (!isTriangularSet(I))
    {
      list triangI = triangL(I);
    }
    else
    {
      list triangI = I;
    }
  }
  else
  {
    list triangI = I;
  }
  int i,ii,j,jj,k,l;
  list Vlocal;
  ideal J;
  list VV;
  for (ii=1; ii<=size(triangI); ii++)
  {
    J = triangI[ii];
    Vlocal = newtonPolygonNegSlopes(J[1],1);
    for (jj=1; jj<=size(Vlocal); jj++)
    {
      list Vjj;
      l = lexSmallestVariableIndex(J[1]);
      Vjj[l] = Vlocal[jj];
      VV[jj] = Vjj;
      kill Vjj;
    }

    for (i=2; i<=size(J); i++)
    {
      k = size(VV);
      for (jj=1; jj<=k; jj++)
      {
        list Vlocal(jj) = expectedNewtonPolygonNegSlopes(J[i],VV[jj],1);
      }
      l = lexSmallestVariableIndex(J[i]);
      for (jj=1; jj<=k; jj++)
      {
        for (j=1; j<=size(Vlocal(jj)); j++)
        {
          if (j==1)
          {
            VV[jj][l] = Vlocal(jj)[1];
          }
          else
          {
            VV[size(VV)+1] = VV[jj];
            VV[size(VV)][l] = Vlocal(jj)[j];
          }
        }
      }
      for (jj=1; jj<=k; jj++)
      {
        kill Vlocal(jj);
      }
    }
  }
  list TI;
  list VVjj;
  for (jj=1; jj<=size(VV); jj++)
  {
    VVjj = VV[jj];
    matrix Vjj[1][size(VVjj)] = VVjj[1..size(VVjj)];
    TI[jj] = -Vjj;
    kill Vjj;
  }

  return (TI);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = (0,t),(z,y,x),lp;
  number uniformizingParameter = t;
  export(uniformizingParameter);
  ideal I = tx2+x+1,txy2+xy+1,xyz+1;
  list TI = tropicalVarietyNewton(I);
  for (int i=1; i<=size(TI); i++)
  { print(TI[i]); }
  kill uniformizingParameter;
}


proc switchRingsAndComputeInitialIdeal(ideal I, intvec w)
"USAGE:   switchRingsAndComputeInitialIdeal(I,w); I ideal, w intvec
RETURN:  ring, a ring containing the initial ideal with respect to w
ASSUME:  constant coefficient case and w strictly positive integer
NOTE:    if printlevel sufficiently high will print timing
EXAMPLE: example switchRingsAndComputeInitialIdeal; shows an example
"
{
  def origin = basering;
  intvec wOne = oneVector(nvars(origin));
  execute("ring rInitialIdeal = ("+charstr(origin)+"),("+varstr(origin)+"),(a(wOne),wp(w));");
  ideal I = imap(origin,I);
  int tinI = timer;
  option(redSB);
  ideal stdI = satstd(I);
  ideal inI = initial(stdI,w);
  tinI = timer-tinI;
  dbprint("time used computing initial ideal: "+string(tinI));

  export(I);
  export(stdI);
  export(inI);
  return (rInitialIdeal);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = (0,t),(p01,p02,p12,p03,p13,p23,p04,p14,p24,p34),dp;
  number uniformizingParameter = t;
  export(uniformizingParameter);
  ideal I =
    p23*p14-p13*p24+p12*p34,
    p23*p04-p03*p24+p02*p34,
    p13*p04-p03*p14+p01*p34,
    p12*p04-p02*p14+p01*p24,
    p12*p03-p02*p13+p01*p23;
  system("--random",1337);
  printlevel = 3;
  matrix p = tropicalPointNewton(I);
  print(p);
  intvec w = cccMatrixToPositiveIntvec(p);
  print(w);
  def s = switchRingsAndComputeInitialIdeal(I,w);
  kill uniformizingParameter;
}

static proc pivotIndices(matrix H)
{
  intvec p;
  p[ncols(H)]=0;
  int pp;
  int i,j;
  for (i=1; i<=ncols(H); i++)
  {
    for (j=nrows(H); j>=0; j--)
    {
      if (H[j,i]!=0)
      {
        if (j>pp)
        {
          p[i] = 1;
          pp = j;
        }
        break;
      }
    }
  }
  return (p);
}


static proc varstrIntvec(intvec p)
{
  string s;
  for (int i=1; i<=size(p); i++)
  {
    s = s+varstr(p[i])+",";
  }
  s = s[1..size(s)-1];
  return (s);
}


static proc substRing(int i, string orderstring)
{
  int n = nvars(basering);
  if (i==1)
  {
    intvec p = 2..n;
  }
  else
  {
    if (i==n)
    {
      intvec p = 1..n-1;
    }
    else
    {
      intvec p = 1..i-1,i+1..n;
    }
  }
  execute("ring ssub = (0,t),("+varstrIntvec(p)+"),"+orderstring+";");
  return (ssub);
}


static proc extendTropNewton(matrix Ti, int i, number toAdd, intvec toFill)
{
  // extend TI by one, inserting toAdd in position i
  int n = ncols(Ti)+1;
  matrix Tii[1][n];
  Tii[1,i] = toAdd;
  if (i==1)
  {
    Tii[1,2..n] = Ti[1,1..n-1];
  }
  else
  {
    if (i==n)
    {
      Tii[1,1..n-1] = Ti[1,1..n-1];
    }
    else
    {
      Tii[1,1..i-1] = Ti[1,1..i-1];
      Tii[1,i+1..n] = Ti[1,i..n-1];
    }
  }
  // extend Tii, inserting 0 in the 1 positions of toFill
  n = n+sum(toFill);
  matrix Tiii[1][n];
  int TiiCounter=1;
  for (int j=1; j<=n; j++)
  {
    if (toFill[j]==1)
    {
      Tiii[1,j]=0;
    }
    else
    {
      Tiii[1,j]=Tii[1,TiiCounter];
      TiiCounter++;
    }
  }
  return (Tiii);
}


static proc mergeTropNewton(list T, list Ti, int i, number toAdd, intvec toFill)
{
  int ii,j;
  matrix Tii;
  for (ii=1; ii<=size(Ti); ii++)
  {
    Tii=extendTropNewton(Ti[ii],i,toAdd,toFill);
    for (j=1; j<=size(T); j++)
    {
      if (T[j]==Tii)
      {
        break;
      }
    }
    if (j>size(T))
    {
      T[size(T)+1] = Tii;
    }
  }
  return (T);
}


proc tropicalLinkNewton(ideal inI)
"USAGE:   tropicalLinkNewton(inI); inI ideal
RETURN:  matrix, a matrix containing generators of all rays of the tropical variety
ASSUME:  constant coefficient case, inI is monomial free,
         its tropical variety has codimension one lineality space and
         is a polyhedral fan
NOTE:    if printlevel sufficiently high will print intermediate results
EXAMPLE: example tropicalLinkNewton; shows an example
"
{
  ring origin=basering;
  dbprint("reducing to one-dimensional fan");
  cone C0 = homogeneitySpace(inI);
  intmat HH = intmat(generatorsOfLinealitySpace(C0));
  matrix H = gauss_nf(HH);

  intvec p = pivotIndices(H);
  string resVars;
  ideal resImage;
  for (int i=1; i<=nvars(basering); i++)
  {
    if (p[i]==1)
    {
      resImage[i]=1;
    }
    else
    {
      resImage[i]=var(i);
      resVars = resVars+varstr(i)+",";
    }
  }
  resVars = resVars[1..size(resVars)-1];

  execute("ring srestr = ("+string(char(origin))+",t),("+resVars+"),dp;");
  number uniformizingParameter = t;
  map resMap = origin,imap(origin,resImage);
  ideal inI = resMap(inI);
  inI = satstd(inI);

  dbprint("intersecting with pairs of affine hyperplanes");
  ideal substImagePos;
  ideal substImageNeg;
  int n = nvars(srestr);
  list T;
  int j;
  for (i=1; i<=n; i++)
  {
    setring srestr;
    substImagePos = maxideal(1);
    substImagePos[i] = t;
    substImageNeg = maxideal(1);
    substImageNeg[i] = t^(-1);

    def ssubstDP = substRing(i,"dp");
    setring ssubstDP;
    ideal substImagePos = imap(srestr,substImagePos);
    map substMapPos = srestr,substImagePos;
    ideal inIPos = substMapPos(inI);
    if (dim(std(inIPos))<0)
    {
      setring srestr;
      kill ssubstDP;
      dbprint(string(i)+": empty");
      i++;
      continue;
    }
    dbprint(string(i)+": non-empty, computing tropical variety");

    setring srestr;
    kill ssubstDP;
    def ssubstLP = substRing(i,"lp");
    setring ssubstLP;
    ideal substImagePos = imap(srestr,substImagePos);
    map substMapPos = srestr,substImagePos;
    ideal inIPos = std(substMapPos(inI));
    number uniformizingParameter=t;
    export(uniformizingParameter);
    list Tpos = tropicalVarietyNewton(inIPos);
    ideal substImageNeg = imap(srestr,substImageNeg);
    map substMapNeg = srestr,substImageNeg;
    ideal inINeg = std(substMapNeg(inI));
    list Tneg = tropicalVarietyNewton(inINeg);
    setring srestr;
    list Tpos = imap(ssubstLP,Tpos);
    list Tneg = imap(ssubstLP,Tneg);
    kill ssubstLP;
    T = mergeTropNewton(T,Tpos,i,number(1),p);
    T = mergeTropNewton(T,Tneg,i,number(-1),p);
    dbprint("total number of rays: "+string(size(T)));
    kill Tpos;
    kill Tneg;
  }
  setring origin;
  return(imap(srestr,T));
}
example
{ "EXAMPLE:"; echo = 2;
  // a 10 valent facet in tropical Grass(3,7)
  ring r = (0,t),
           (p012,p013,p023,p123,p014,p024,p124,p034,p134,p234,
            p015,p025,p125,p035,p135,p235,p045,p145,p245,p345,
            p016,p026,p126,p036,p136,p236,p046,p146,p246,p346,
            p056,p156,p256,p356,p456),
           wp(4,7,5,7,4,4,4,7,5,7,2,1,2,4,4,4,2,1,2,4,7,5,7,7,
              5,7,7,5,7,4,4,4,4,4,4);
  number uniformizingParameter = t;
  export(uniformizingParameter);
  ideal inI =
    p345*p136+p134*p356,  p125*p045+p015*p245,  p124*p015-p014*p125,
    p135*p245-p125*p345,  p135*p045+p015*p345,  p124*p045+p014*p245,
    p024*p125-p012*p245,  p145*p236-p124*p356,  p124*p135-p123*p145,
    p024*p015+p012*p045,  p134*p026+p023*p146-p024*p136,
    p145*p036+p014*p356,  p014*p135-p013*p145,  p234*p145+p124*p345,
    p034*p145-p014*p345,  p024*p135-p012*p345,  p125*p035+p015*p235,
    p235*p045-p035*p245,  p234*p136-p134*p236,  p134*p036-p034*p136,
    p146*p356-p136*p456,  p135*p146-p134*p156,
    p135*p026+p023*p156+p012*p356,  p124*p035+p014*p235,
    p123*p025+p012*p235,  p013*p025-p012*p035,  p345*p146+p134*p456,
    p125*p036+p015*p236,  p345*p026-p023*p456+p024*p356,
    p123*p015-p013*p125,  p234*p025-p024*p235,  p034*p025-p024*p035,
    p234*p125+p123*p245,  p245*p036-p045*p236,  p123*p045+p013*p245,
    p034*p125-p013*p245,  p234*p015+p013*p245,  p245*p156+p125*p456,
    p034*p015+p013*p045,  p045*p156-p015*p456,  p135*p236-p123*p356,
    p235*p146-p134*p256,  p135*p036+p013*p356,  p124*p036+p014*p236,
    p123*p014-p013*p124,  p035*p146-p134*p056,  p145*p126+p124*p156,
    p234*p045-p034*p245,  p235*p026+p023*p256-p025*p236,
    p145*p016+p014*p156,  p035*p026+p023*p056-p025*p036,
    p345*p236+p234*p356,  p234*p135+p123*p345,  p345*p036+p034*p356,
    p034*p135-p013*p345,  p345*p156+p135*p456,  p124*p034+p014*p234,
    p145*p246-p124*p456,  p123*p024+p012*p234,  p145*p046+p014*p456,
    p013*p024-p012*p034,  p024*p156+p012*p456,  p125*p056+p015*p256,
    p245*p056-p045*p256,  p236*p146-p136*p246,  p134*p126+p123*p146,
    p136*p046-p036*p146,  p235*p036-p035*p236,  p134*p016+p013*p146,
    p123*p035+p013*p235,  p235*p156-p135*p256,
    p123*p026-p023*p126+p012*p236,  p135*p056-p035*p156,
    p023*p016-p013*p026+p012*p036,  p124*p056+p014*p256,
    p234*p146-p134*p246,  p025*p126-p012*p256,  p134*p046-p034*p146,
    p025*p016+p012*p056,  p234*p035-p034*p235,  p345*p256+p235*p456,
    p234*p026+p023*p246-p024*p236,  p345*p056+p035*p456,
    p034*p026+p023*p046-p024*p036,  p125*p016-p015*p126,
    p025*p246-p024*p256,  p025*p046-p024*p056,  p245*p126-p125*p246,
    p125*p046+p015*p246,  p045*p126+p015*p246,  p245*p016-p015*p246,
    p045*p016-p015*p046,  p123*p036+p013*p236,  p236*p156+p126*p356,
    p135*p126+p123*p156,  p036*p156-p016*p356,  p135*p016+p013*p156,
    p124*p016-p014*p126,  p235*p056-p035*p256,  p245*p046-p045*p246,
    p234*p036-p034*p236,  p123*p034+p013*p234,  p246*p356-p236*p456,
    p234*p156-p123*p456,  p135*p246-p123*p456,  p345*p126-p123*p456,
    p046*p356-p036*p456,  p034*p156+p013*p456,  p135*p046+p013*p456,
    p345*p016-p013*p456,  p124*p046+p014*p246,  p024*p126-p012*p246,
    p024*p016+p012*p046,  p345*p246+p234*p456,  p345*p046+p034*p456,
    p235*p126+p123*p256,  p236*p056-p036*p256,  p123*p056+p013*p256,
    p035*p126-p013*p256,  p235*p016+p013*p256,  p035*p016+p013*p056,
    p235*p246-p234*p256,  p234*p056-p034*p256,  p035*p246-p034*p256,
    p235*p046-p034*p256,  p035*p046-p034*p056,  p126*p036+p016*p236,
    p123*p016-p013*p126,  p234*p126+p123*p246,  p236*p046-p036*p246,
    p123*p046+p013*p246,  p034*p126-p013*p246,  p234*p016+p013*p246,
    p246*p156+p126*p456,  p034*p016+p013*p046,  p046*p156-p016*p456,
    p234*p046-p034*p246,  p126*p056+p016*p256,  p246*p056-p046*p256,
    p126*p046+p016*p246,  p024*p235*p145+p124*p025*p345,
    p024*p035*p145-p014*p025*p345,  p123*p145*p245-p124*p125*p345,
    p013*p145*p245-p014*p125*p345,  p013*p045*p145+p014*p015*p345,
    p024*p235*p136-p134*p025*p236,  p123*p245*p136+p134*p125*p236,
    p013*p245*p136+p134*p015*p236,  p034*p245*p136-p134*p045*p236,
    p134*p156*p356-p135*p136*p456,  p123*p145*p146-p124*p134*p156,
    p013*p145*p146-p014*p134*p156,  p013*p145*p026+p023*p014*p156+p012*p014*p356,
    p124*p025*p156+p012*p145*p256,  p012*p145*p056-p014*p025*p156,
    p024*p145*p256-p124*p025*p456,  p024*p145*p056+p014*p025*p456,
    p034*p235*p136-p134*p035*p236,  p134*p256*p356-p235*p136*p456,
    p134*p056*p356-p035*p136*p456,  p025*p036*p146-p024*p136*p056,
    p013*p125*p026-p023*p015*p126+p012*p015*p236,
    p123*p245*p146+p134*p125*p246,  p013*p245*p146+p134*p015*p246,
    p013*p245*p026-p023*p015*p246-p012*p045*p236,
    p013*p045*p026-p023*p015*p046-p012*p045*p036,
    p034*p245*p146-p134*p045*p246,  p013*p124*p026-p023*p014*p126+p012*p014*p236,
    p013*p145*p056-p014*p035*p156,  p024*p256*p356-p025*p236*p456,
    p024*p056*p356-p025*p036*p456,  p234*p256*p356-p235*p236*p456,
    p034*p256*p356-p035*p236*p456,  p034*p056*p356-p035*p036*p456,
    p012*p235*p145*p245+p124*p025*p125*p345,
    p012*p035*p145*p245-p014*p025*p125*p345,
    p012*p035*p045*p145+p014*p015*p025*p345,
    p012*p235*p245*p136-p134*p025*p125*p236,
    p012*p035*p245*p136+p134*p015*p025*p236,
    p024*p035*p245*p136-p134*p025*p045*p236,
    p014*p025*p125*p156+p012*p015*p145*p256,
    p012*p145*p245*p256-p124*p025*p125*p456,
    p012*p045*p145*p256+p014*p025*p125*p456,
    p012*p245*p256*p356-p025*p125*p236*p456,
    p012*p045*p256*p356+p015*p025*p236*p456,
    p012*p045*p056*p356+p015*p025*p036*p456,
    p123*p245*p256*p356+p125*p235*p236*p456,
    p013*p245*p256*p356+p015*p235*p236*p456,
    p013*p045*p256*p356+p015*p035*p236*p456,
    p013*p045*p056*p356+p015*p035*p036*p456;
  system("--random",1337);
  printlevel = 3;
  list TinI = tropicalLinkNewton(inI);
  for (int i=1; i<=size(TinI); i++)
  { print(TinI[i]); }
}


// disabled routines to check characteristic-freeness of tropical points

// proc saturateWithRespectToVariable(ideal I, int k)
// {
//   ASSUME(1,k>=1);
//   ASSUME(1,k<=nvars(basering));

//   def origin = basering;
//   int n = nvars(basering);
//   intvec weightVector = ringlist(origin)[3][1][2];

//   string newVars;
//   for (int i=1; i<k; i++)
//   {
//     newVars = newVars+string(var(i))+",";
//   }
//   for (i=k+1; i<=n; i++)
//   {
//     newVars = newVars+string(var(i))+",";
//   }
//   newVars = newVars+string(var(k));
//   execute("ring ringForSaturation = ("+charstr(origin)+"),("+newVars+"),dp;");

//   ideal I = satstd(imap(origin,I));
//   if (I==-1)
//   {
//     return (-1);
//   }
//   I = simplify(I,2+4+32);

//   setring origin;
//   I = imap(ringForSaturation,I);
//   return (I);
// }

// proc stepwiseSaturation(ideal I)
// {
//   if (I!=1)
//   {
//     list variablesToBeSaturated;
//     int l = nvars(basering);
//     for (int i=1; i<=l; i++)
//     { variablesToBeSaturated[i]=l-i+1; }

//     while (size(variablesToBeSaturated)>0)
//     {
//       dbprint("variablesToBeSaturated: "+string(variablesToBeSaturated));
//       I = saturateWithRespectToVariable(I,variablesToBeSaturated[1]);
//       variablesToBeSaturated = delete(variablesToBeSaturated,1);
//       if ((I==1) || (I==-1))
//       {
//         break;
//       }
//     }
//   }

//   return (I);
// }


// proc checkForContainmentInTropicalVariety(ideal I, intvec w, int charInt)
// {
//   def origin = basering;
//   intvec wOne = oneVector(nvars(origin));
//   execute("ring rInitialIdeal = ("+string(charInt)+"),("+varstr(origin)+"),(a(wOne),wp(w));");

//   ideal I = imap(origin,I);
//   int tinI = timer;
//   option(redSB);
//   ideal stdI = satstd(I);
//   attrib(stdI,"isSB",1);
//   ideal inI = initial(stdI,w);
//   tinI = timer-tinI;
//   dbprint("time used computing initial ideal: "+string(tinI));

//   int tsat = timer;
//   ideal satinI = stepwiseSaturation(inI);
//   tsat = timer-tsat;
//   dbprint("time used computing saturation: "+string(tsat));


//   export(I);
//   export(stdI);
//   export(inI);
//   export(satinI);
//   return (rInitialIdeal);
// }