source: git/Singular/LIB/tropicalNewton.lib @ 4639d07

fieker-DuValspielwiese
Last change on this file since 4639d07 was e4805bd, checked in by Hans Schoenemann <hannes@…>, 7 years ago
doc: tropicalNewton.lib
  • Property mode set to 100644
File size: 30.2 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="version tropical.lib 4.0.3.3 Nov_2016 "; // $Id$
3category="Tropical Geometry";
4info="
5LIBRARY:         tropicalNewton.lib  Computations in Tropical Geometry
6                                     using Newton Polygon methods
7
8AUTHORS:         Tommy Hofman,       email: thofmann@mathematik.uni.kl.de
9                 Yue Ren,            email: reny@cs.bgu.ac.il
10
11OVERVIEW:
12  This libraries contains algorithms for computing
13  - non-trivial points on tropical varieties,
14  - zero-dimensional tropical varieties,
15  - one-codimensional links of tropical varieties
16  based on Newton polygon methods.
17
18REFERENCES: Hofmann, Ren: Computing tropical points and tropical links, arXiv:1611.02878
19            (WARNING: this library follows the max convention instead
20                      and triangular sets follow the definition of the Singular book)
21
22SEE ALSO: tropicalVariety, tropical_lib
23
24KEYWORDS: tropical geometry; tropical varieties; Newton polygons
25
26PROCEDURES:
27  setUniformizingParameter()  sets the uniformizingParameter
28  val()                       returns valuation of element in ground field
29  newtonPolygonNegSlopes()    returns negative of the Newton Polyong slopes
30  cccMatrixToPositiveIntvec() helper function to turn a computed valuation vector
31                              into a usable weight vector in Singular
32  tropicalPointNewton()       computes point on tropical variety
33  switchRingsAndComputeInitialIdeal()
34                              switches rings and computes initial ideal
35  tropicalVarietyNewton()     computes tropical variety of zero-dimensional ideal
36  tropicalLinkNewton()        computes tropical variety that is polyhedral fan
37                                and has codimension one lineality space
38";
39
40///////////////////////////////////////////////////////////////////////////////
41
42LIB "gfan.lib";
43LIB "triang.lib";
44LIB "linalg.lib";
45
46///////////////////////////////////////////////////////////////////////////////
47
48proc setUniformizingParameter(number p)
49"USAGE:   setUniformizingParameter(p); p number
50RETURN:  none, sets the uniformizing parameter as p
51ASSUME:  char(K)==0 and p prime number or
52         trdeg(K)>0 and p transcendental variable or
53         p==0
54EXAMPLE: example setUniformizingParameter; shows an example
55"
56{
57  // kill uniformizingParemeter if previously defined
58  if (defined(uniformizingParameter))
59  {
60    kill uniformizingParameter;
61  }
62
63  // input sanity check
64  if (p!=0)
65  {
66    if (pardeg(p)==0)
67    {
68      if(char(basering)!=0 && prime(int(p))!=p)
69      {
70        ERROR("setUniformizingParameter: unexpected parameters");
71      }
72    }
73    else
74    {
75      if (size(p)!=2 || pardeg(p)!=1)
76      {
77        ERROR("setUniformizingParameter: unexpected parameters");
78      }
79    }
80  }
81
82  // set uniformizingParameter and export it
83  number uniformizingParameter = p;
84  export(uniformizingParameter);
85}
86example
87{ "EXAMPLE:"; echo = 2;
88  // poor man's polynomials over Puiseux series:
89  ring r = (0,t),x,dp;
90  setUniformizingParameter(t);
91  val(t2+t3);
92  val(t^-2+t^-3);
93  // poor man's polynomials over p-adic numbers:
94  ring s = 0,x,dp;
95  setUniformizingParameter(2);
96  val(12);
97  val(1/12);
98}
99
100proc val(number c)
101"USAGE:   val(c); c number
102RETURN:  int, the valuation of a element in the ground field
103ASSUME:  uniformizingParameter is set and c!=0
104EXAMPLE: example val; shows an example
105"
106{
107  if (defined(uniformizingParameter)==0)
108  {
109    ERROR("val: uniformizingParameter not set");
110  }
111  if (c==0)
112  {
113    ERROR("val: input is 0");
114  }
115
116  if (uniformizingParameter==0)
117  {
118    return (0);
119  }
120
121  int vc = 0;
122  if (pardeg(uniformizingParameter)>0)
123  {
124    def origin = basering;
125    number c_denom = denominator(c);
126    number c_num = numerator(c);
127    execute("ring r_Val=0,"+string(uniformizingParameter)+",ds");
128    poly c_denom = imap(origin,c_denom);
129    poly c_num = imap(origin,c_num);
130    vc = ord(c_num)-ord(c_denom);
131    setring origin;
132  }
133  else
134  {
135    int c_denom = int(denominator(c));
136    int c_num = int(numerator(c));
137    int p = int(uniformizingParameter);
138
139    while (c_num mod p==0)
140    {
141      c_num = c_num div p;
142      vc++;
143    }
144    while (c_denom mod p==0)
145    {
146      c_denom = c_denom div p;
147      vc--;
148    }
149  }
150
151  return (vc);
152}
153example
154{ "EXAMPLE:"; echo = 2;
155  // poor man's polynomials over Puiseux series:
156  ring r = (0,t),x,dp;
157  setUniformizingParameter(t);
158  val(t2+t3);
159  val(t^-2+t^-3);
160  // poor man's polynomials over p-adic numbers:
161  ring s = 0,x,dp;
162  setUniformizingParameter(2);
163  val(12);
164  val(1/12);
165}
166
167proc newtonPolygonNegSlopes(poly g, list #)
168"USAGE:   newtonPolygonNegSlopes(g,b); g poly, b int
169RETURN:  list, the negative slopes of the Newton Polygon of g
170         if b==1, computes root (type poly) instead if (easily) possible
171ASSUME:  uniformizingParameter is set and g univariate
172EXAMPLE: example newtonPolygonNegSlopes; shows an example
173"
174{
175  if (size(#)>0)
176  {
177    if (#[1]==1 && size(g)==2)
178    {
179      if (deg(g[1])==1 && deg(g[2])==0)
180      {
181        poly root = -leadcoef(g[2])/leadcoef(g[1]);
182        return (list(root));
183      }
184    }
185  }
186
187  if (defined(uniformizingParameter)==0)
188  {
189    ERROR("newtonPolygonNegSlopes: uniformizingParameter not set");
190  }
191  int k = size(g);
192  intmat M [k+1][3];
193  for (int i=1; i<=k; i++)
194  {
195    M[i,1] = 1;
196    M[i,2] = deg(g[i]);
197    M[i,3] = val(leadcoef(g[i]));
198  }
199  M[k+1,1] = 0;
200  M[k+1,2] = 0;
201  M[k+1,3] = 1;
202  cone Delta = coneViaPoints(M);
203  bigintmat F = facets(Delta);
204
205  list slopes;
206  number slope;
207  for (i=1; i<=nrows(F); i++)
208  {
209    if (F[i,3]!=0)
210    {
211      slope = number(F[i,2])/number(F[i,3]);
212      slopes[size(slopes)+1] = slope;
213    }
214  }
215  return (slopes);
216}
217example
218{ "EXAMPLE:"; echo = 2;
219  ring r = (0,t),x,dp;
220  setUniformizingParameter(t);
221  poly g = tx2+x+1;
222  newtonPolygonNegSlopes(g);
223  // poor man's polynomials over p-adic numbers:
224  ring s = 0,x,dp;
225  setUniformizingParameter(3);
226  poly g = x2+9x+1;
227  newtonPolygonNegSlopes(g);
228}
229
230
231static proc listDot(list VV, intvec expv)
232{
233  number ld;
234  for (int i=1; i<=size(expv); i++)
235  {
236    if (expv[i]>0)
237    {
238      if (typeof(VV[i])=="number")
239      {
240        ld = ld + VV[i]*expv[i];
241      }
242      else
243      {
244        ERROR("listDot: unkown or invalid relevant entry");
245      }
246    }
247  }
248  return (ld);
249}
250
251
252static proc lexSmallestVariableIndex(poly g)
253{
254  intvec alpha = leadexp(g);
255  for (int i=1; i<=nvars(basering); i++)
256  {
257    if (alpha[i]>0)
258    {
259      return (i);
260    }
261  }
262  return (-1);
263}
264
265static proc numbersLessThanAsRationalNumbers(number a, number b)
266{
267  int a_num = int(numerator(a));
268  int a_denom = int(denominator(a));
269  int b_num = int(numerator(b));
270  int b_denom = int(denominator(b));
271  return (a_num*b_denom-b_num*a_denom<0);
272}
273
274static proc expectedValuation(poly h, list VV)
275{
276  number vh = val(leadcoef(h[1]))+listDot(VV,leadexp(h[1]));
277  number vhh;
278  int amb = 0;
279  for (int i=2; i<=size(h); i++)
280  {
281    vhh = val(leadcoef(h[i]))+listDot(VV,leadexp(h[i]));
282    if (vh == vhh)
283    {
284      amb = 1;
285    }
286    if (numbersLessThanAsRationalNumbers(vhh,vh))
287    {
288      vh = vhh;
289      amb = 0;
290    }
291  }
292  if (amb > 0)
293  {
294    "h:"; h;
295    "VV:"; VV;
296    ERROR("expectedValuation: valuation ambiguous");
297  }
298
299  return (vh);
300}
301
302static proc expectedNewtonPolygonNegSlopes(poly g, list VV, list #)
303{
304  int k = size(VV);
305  for (int i=1; i<=size(VV); i++)
306  {
307    if (typeof(VV[i])=="poly")
308    {
309      g = subst(g,var(i),VV[i]);
310    }
311  }
312
313  if (size(#)>0)
314  {
315    if (#[1]==1 && size(g)==2)
316    {
317      if (deg(g[1])==1 && deg(g[2])==0)
318      {
319        poly root = -leadcoef(g[2])/leadcoef(g[1]);
320        return (list(root));
321      }
322    }
323  }
324
325  int n = lexSmallestVariableIndex(g);
326  matrix G = coef(g,var(n));
327
328  k = ncols(G);
329  intmat M[k+1][3];
330  number vh;
331  for (i=1; i<=k; i++)
332  {
333    vh = expectedValuation(G[2,i],VV);
334    M[i,1] = int(denominator(vh));
335    M[i,2] = deg(G[1,i]);
336    M[i,3] = int(numerator(vh));
337  }
338  M[k+1,1] = 0;
339  M[k+1,2] = 0;
340  M[k+1,3] = 1;
341
342  cone Delta = coneViaPoints(M);
343  bigintmat F = facets(Delta);
344
345  list slopes;
346  number slope;
347  for (i=1; i<=nrows(F); i++)
348  {
349    if (F[i,3]!=0)
350    {
351      slope = number(F[i,2])/number(F[i,3]);
352      slopes[size(slopes)+1] = slope;
353    }
354  }
355  return (slopes);
356}
357
358
359static proc randomValuedNumber(list #)
360{
361  int n = 99;
362  if (size(#)>0)
363  {
364    n = #[1];
365  }
366  int v = random(0,n);
367
368  number c = random(1,99)*uniformizingParameter^v;
369  c = c+random(0,99)*uniformizingParameter^(v+random(0,n div 3));
370  c = c+random(0,99)*uniformizingParameter^(v+random(0,n div 3));
371
372  return (c);
373}
374
375
376static proc listOfNumbersToMatrix(list L)
377{
378  int k = size(L);
379  matrix Lmat[1][k];
380  for (int i=1; i<=k; i++)
381  {
382    if (typeof(L[i]) == "number")
383    {
384      Lmat[1,i] = L[i];
385    }
386    if (typeof(L[i]) == "poly")
387    {
388      Lmat[1,i] = val(number(L[i]));
389    }
390  }
391
392  return (Lmat);
393}
394
395proc cccMatrixToPositiveIntvec(matrix L)
396"USAGE:   cccMatrixToPositiveIntvec(M); M matrix
397RETURN:  intvec, strictly positive equivalent as weight vector to row vector in M
398ASSUME:  constant coefficient case only,
399         will scale weight vector and add vectors of ones to it
400EXAMPLE: example cccMatrixToPositiveIntvec; shows an example
401"
402{
403  int k = ncols(L);
404  intvec denoms;
405  intvec enums;
406  int denoms_gcd = int(denominator(number(L[1,1])));
407  int denoms_prod = int(denominator(number(L[1,1])));
408  for (int i=2; i<=k; i++)
409  {
410    denoms_gcd = gcd(denoms_gcd,int(denominator(number(L[1,i]))));
411    denoms_prod = denoms_prod*int(denominator(number(L[1,i])));
412  }
413  int denoms_lcm = denoms_prod div denoms_gcd;
414
415  intvec w;
416  for (i=1; i<=k; i++)
417  {
418    w[i] = int(number(L[1,i])*denoms_lcm);
419  }
420  return (addOneVectorUntilPositive(w));
421}
422example
423{ "EXAMPLE:"; echo = 2;
424  ring r = (0,t),(p01,p02,p12,p03,p13,p23,p04,p14,p24,p34),dp;
425  number uniformizingParameter = t;
426  export(uniformizingParameter);
427  ideal I =
428    p23*p14-p13*p24+p12*p34,
429    p23*p04-p03*p24+p02*p34,
430    p13*p04-p03*p14+p01*p34,
431    p12*p04-p02*p14+p01*p24,
432    p12*p03-p02*p13+p01*p23;
433  system("--random",1337);
434  matrix p = tropicalPointNewton(I);
435  print(p);
436  intvec w = cccMatrixToPositiveIntvec(p);
437  print(w);
438  def s = switchRingsAndComputeInitialIdeal(I,w);
439  kill uniformizingParameter;
440}
441
442static proc oneVector(int n)
443{
444  intvec w;
445  for (int i=1; i<=n; i++)
446  {
447    w[i] = 1;
448  }
449  return (w);
450}
451
452static proc addOneVectorUntilPositive(intvec w)
453{
454  int n = size(w);
455  int w_min = w[1];
456  for (int i=2; i<=n; i++)
457  {
458    if (w[i]<w_min)
459    {
460      w_min = w[i];
461    }
462  }
463  return (w-(w_min-1)*oneVector(n));
464}
465
466proc tropicalPointNewton(ideal I, list #)
467"USAGE:   tropicalPointsLasVegas(I); I ideal
468RETURN:  matrix, a matrix containing a tropical point as row vector
469ASSUME:  uniformizingParameter is set and I monomial free
470NOTE:    if printlevel sufficiently high will print intermediate data and timings
471         returns error if randomly chosen hyperplanes are not generic
472EXAMPLE: example tropicalPointNewton; shows an example
473"
474{
475  if (defined(uniformizingParameter)==0)
476  {
477    ERROR("tropicalPointNewton: uniformizingParameter not set");
478  }
479  int ttotal = timer;
480  ring origin = basering;
481
482  int tindepSet = timer;
483  ideal stdI = std(I);
484  attrib(stdI,"isSB",1);
485  intvec u = indepSet(stdI);
486  tindepSet = timer - tindepSet;
487  dbprint("maximal independent set: "+string(u));
488
489  int n = nvars(origin);
490  string vars;
491  string subststr = "ideal J = subst(I,";
492  list VV;
493  number c;
494  for (int i=1; i<=n; i++)
495  {
496    if (u[i]==0)
497    {
498      vars = vars+varstr(i)+",";
499    }
500    else
501    {
502      c = randomValuedNumber(#);
503      VV[i]=number(val(c));
504      dbprint("substituting "+string(var(i))+" with number of valuation "+string(VV[i]));
505      subststr = subststr+varstr(i)+","+string(c)+",";
506    }
507  }
508  vars = vars[1..size(vars)-1];
509  subststr = string(subststr[1..size(subststr)-1])+");";
510  execute(subststr);
511
512  execute("ring s = ("+charstr(origin)+"),("+vars+"),lp;");
513  ideal J = imap(origin,J);
514  int tstdlp = timer;
515  dbprint("computing triangular decomposition (picking first factor)");
516  J = satstd(J);
517  attrib(J,"isSB",1);
518  if (dim(J)>0)
519  {
520    dbprint("not zero-dimensional, please restart algorithm");
521    bigintmat O[1][1];
522    return (O);
523  }
524  if (!isTriangularSet(J))
525  {
526    J = triangL(J)[1];
527  }
528  tstdlp = timer - tstdlp;
529
530  dbprint("starting analysis of Newton polygons");
531  setring origin;
532  J = imap(s,J);
533
534  int tnewton = timer;
535  list Vlocal = newtonPolygonNegSlopes(J[1],1);
536  int l = lexSmallestVariableIndex(J[1]);
537  VV[l] = Vlocal[1];
538  dbprint("possible valuations for "+string(var(l))+" (picking first): "+string(Vlocal));
539
540  for (i=2; i<=size(J); i++)
541  {
542    Vlocal = expectedNewtonPolygonNegSlopes(J[i],VV,1);
543    l = lexSmallestVariableIndex(J[i]);
544    dbprint("possible valuations for "+string(var(l))+" (picking first): "+string(Vlocal));
545    VV[l] = Vlocal[1];
546  }
547  tnewton = timer-tnewton;
548  matrix w = -listOfNumbersToMatrix(VV);
549
550  ttotal = timer-ttotal;
551  dbprint("time used total: "+string(ttotal));
552  dbprint("computing independent set: "+string(tindepSet));
553  dbprint("computing triangular decomposition: "+string(tstdlp));
554  dbprint("analyzing newton polygons: "+string(tnewton));
555
556  return (w);
557}
558example
559{ "EXAMPLE:"; echo = 2;
560  ring r = (0,t),(p01,p02,p12,p03,p13,p23,p04,p14,p24,p34),dp;
561  number uniformizingParameter = t;
562  export(uniformizingParameter);
563  ideal I =
564    p23*p14-p13*p24+p12*p34,
565    p23*p04-p03*p24+p02*p34,
566    p13*p04-p03*p14+p01*p34,
567    p12*p04-p02*p14+p01*p24,
568    p12*p03-p02*p13+p01*p23;
569  system("--random",1337);
570  printlevel = 3;
571  matrix p = tropicalPointNewton(I);
572  print(p);
573  intvec w = cccMatrixToPositiveIntvec(p);
574  print(w);
575  def s = switchRingsAndComputeInitialIdeal(I,w);
576  kill uniformizingParameter;
577}
578
579
580static proc sumOfLeadExpV(poly f)
581{
582  intvec expvSum;
583  expvSum = leadexp(f[1]);
584  for (int i=2; i<=size(f); i++)
585  {
586    expvSum = expvSum + leadexp(f[i]);
587  }
588  return (expvSum);
589}
590
591
592static proc isTriangularSet(ideal F)
593{
594  // checks whether F has the right amount of elements
595  int n = size(F);
596  if (n!=nvars(basering))
597  {
598    return (0);
599  }
600  int i,j;
601  poly fi;
602  intvec expvSum;
603  for (i=1; i<=n; i++)
604  {
605    fi = F[i];
606    expvSum = sumOfLeadExpV(fi);
607    // checks whether fi has monomial containing x_i
608    if (expvSum[n-i+1]==0)
609    {
610      return (0);
611    }
612    // checks whether fi has no monomial containing x_j, j>i
613    for (j=i+1; j<=n; j++)
614    {
615      if (expvSum[n-j+1]>0)
616      {
617        return (0);
618      }
619    }
620  }
621  return (1);
622}
623
624
625proc tropicalVarietyNewton(ideal I)
626"USAGE:   tropicalVarietyNewton(I); I ideal
627RETURN:  matrix, a matrix containing all elements of the tropical variety
628ASSUME:  uniformizingParameter is set, I monomial free and zero-dimensional
629EXAMPLE: example tropicalVarietyNewton; shows an example
630"
631{
632  if (defined(uniformizingParameter)==0)
633  {
634    ERROR("tropicalVarietyNewton: uniformizingParameter not set");
635  }
636  if(ordstr(basering)[1,2]!="lp")
637  {
638    ERROR("tropicalVarietyNewton: order must be lp");
639  }
640
641  if (!isTriangularSet(I))
642  {
643    I = satstd(I);
644    attrib(I,"isSB",1);
645    if (!isTriangularSet(I))
646    {
647      list triangI = triangL(I);
648    }
649    else
650    {
651      list triangI = I;
652    }
653  }
654  else
655  {
656    list triangI = I;
657  }
658  int i,ii,j,jj,k,l;
659  list Vlocal;
660  ideal J;
661  list VV;
662  for (ii=1; ii<=size(triangI); ii++)
663  {
664    J = triangI[ii];
665    Vlocal = newtonPolygonNegSlopes(J[1],1);
666    for (jj=1; jj<=size(Vlocal); jj++)
667    {
668      list Vjj;
669      l = lexSmallestVariableIndex(J[1]);
670      Vjj[l] = Vlocal[jj];
671      VV[jj] = Vjj;
672      kill Vjj;
673    }
674
675    for (i=2; i<=size(J); i++)
676    {
677      k = size(VV);
678      for (jj=1; jj<=k; jj++)
679      {
680        list Vlocal(jj) = expectedNewtonPolygonNegSlopes(J[i],VV[jj],1);
681      }
682      l = lexSmallestVariableIndex(J[i]);
683      for (jj=1; jj<=k; jj++)
684      {
685        for (j=1; j<=size(Vlocal(jj)); j++)
686        {
687          if (j==1)
688          {
689            VV[jj][l] = Vlocal(jj)[1];
690          }
691          else
692          {
693            VV[size(VV)+1] = VV[jj];
694            VV[size(VV)][l] = Vlocal(jj)[j];
695          }
696        }
697      }
698      for (jj=1; jj<=k; jj++)
699      {
700        kill Vlocal(jj);
701      }
702    }
703  }
704  list TI;
705  list VVjj;
706  for (jj=1; jj<=size(VV); jj++)
707  {
708    VVjj = VV[jj];
709    matrix Vjj[1][size(VVjj)] = VVjj[1..size(VVjj)];
710    TI[jj] = -Vjj;
711    kill Vjj;
712  }
713
714  return (TI);
715}
716example
717{ "EXAMPLE:"; echo = 2;
718  ring r = (0,t),(z,y,x),lp;
719  number uniformizingParameter = t;
720  export(uniformizingParameter);
721  ideal I = tx2+x+1,txy2+xy+1,xyz+1;
722  list TI = tropicalVarietyNewton(I);
723  for (int i=1; i<=size(TI); i++)
724  { print(TI[i]); }
725  kill uniformizingParameter;
726}
727
728
729proc switchRingsAndComputeInitialIdeal(ideal I, intvec w)
730"USAGE:   switchRingsAndComputeInitialIdeal(I,w); I ideal, w intvec
731RETURN:  ring, a ring containing the initial ideal with respect to w
732ASSUME:  constant coefficient case and w strictly positive integer
733NOTE:    if printlevel sufficiently high will print timing
734EXAMPLE: example switchRingsAndComputeInitialIdeal; shows an example
735"
736{
737  def origin = basering;
738  intvec wOne = oneVector(nvars(origin));
739  execute("ring rInitialIdeal = ("+charstr(origin)+"),("+varstr(origin)+"),(a(wOne),wp(w));");
740  ideal I = imap(origin,I);
741  int tinI = timer;
742  option(redSB);
743  ideal stdI = satstd(I);
744  ideal inI = initial(stdI,w);
745  tinI = timer-tinI;
746  dbprint("time used computing initial ideal: "+string(tinI));
747
748  export(I);
749  export(stdI);
750  export(inI);
751  return (rInitialIdeal);
752}
753example
754{ "EXAMPLE:"; echo = 2;
755  ring r = (0,t),(p01,p02,p12,p03,p13,p23,p04,p14,p24,p34),dp;
756  number uniformizingParameter = t;
757  export(uniformizingParameter);
758  ideal I =
759    p23*p14-p13*p24+p12*p34,
760    p23*p04-p03*p24+p02*p34,
761    p13*p04-p03*p14+p01*p34,
762    p12*p04-p02*p14+p01*p24,
763    p12*p03-p02*p13+p01*p23;
764  system("--random",1337);
765  printlevel = 3;
766  matrix p = tropicalPointNewton(I);
767  print(p);
768  intvec w = cccMatrixToPositiveIntvec(p);
769  print(w);
770  def s = switchRingsAndComputeInitialIdeal(I,w);
771  kill uniformizingParameter;
772}
773
774static proc pivotIndices(matrix H)
775{
776  intvec p;
777  p[ncols(H)]=0;
778  int pp;
779  int i,j;
780  for (i=1; i<=ncols(H); i++)
781  {
782    for (j=nrows(H); j>=0; j--)
783    {
784      if (H[j,i]!=0)
785      {
786        if (j>pp)
787        {
788          p[i] = 1;
789          pp = j;
790        }
791        break;
792      }
793    }
794  }
795  return (p);
796}
797
798
799static proc varstrIntvec(intvec p)
800{
801  string s;
802  for (int i=1; i<=size(p); i++)
803  {
804    s = s+varstr(p[i])+",";
805  }
806  s = s[1..size(s)-1];
807  return (s);
808}
809
810
811static proc substRing(int i, string orderstring)
812{
813  int n = nvars(basering);
814  if (i==1)
815  {
816    intvec p = 2..n;
817  }
818  else
819  {
820    if (i==n)
821    {
822      intvec p = 1..n-1;
823    }
824    else
825    {
826      intvec p = 1..i-1,i+1..n;
827    }
828  }
829  execute("ring ssub = (0,t),("+varstrIntvec(p)+"),"+orderstring+";");
830  return (ssub);
831}
832
833
834static proc extendTropNewton(matrix Ti, int i, number toAdd, intvec toFill)
835{
836  // extend TI by one, inserting toAdd in position i
837  int n = ncols(Ti)+1;
838  matrix Tii[1][n];
839  Tii[1,i] = toAdd;
840  if (i==1)
841  {
842    Tii[1,2..n] = Ti[1,1..n-1];
843  }
844  else
845  {
846    if (i==n)
847    {
848      Tii[1,1..n-1] = Ti[1,1..n-1];
849    }
850    else
851    {
852      Tii[1,1..i-1] = Ti[1,1..i-1];
853      Tii[1,i+1..n] = Ti[1,i..n-1];
854    }
855  }
856  // extend Tii, inserting 0 in the 1 positions of toFill
857  n = n+sum(toFill);
858  matrix Tiii[1][n];
859  int TiiCounter=1;
860  for (int j=1; j<=n; j++)
861  {
862    if (toFill[j]==1)
863    {
864      Tiii[1,j]=0;
865    }
866    else
867    {
868      Tiii[1,j]=Tii[1,TiiCounter];
869      TiiCounter++;
870    }
871  }
872  return (Tiii);
873}
874
875
876static proc mergeTropNewton(list T, list Ti, int i, number toAdd, intvec toFill)
877{
878  int ii,j;
879  matrix Tii;
880  for (ii=1; ii<=size(Ti); ii++)
881  {
882    Tii=extendTropNewton(Ti[ii],i,toAdd,toFill);
883    for (j=1; j<=size(T); j++)
884    {
885      if (T[j]==Tii)
886      {
887        break;
888      }
889    }
890    if (j>size(T))
891    {
892      T[size(T)+1] = Tii;
893    }
894  }
895  return (T);
896}
897
898
899proc tropicalLinkNewton(ideal inI)
900"USAGE:   tropicalLinkNewton(inI); inI ideal
901RETURN:  matrix, a matrix containing generators of all rays of the tropical variety
902ASSUME:  constant coefficient case, inI is monomial free,
903         its tropical variety has codimension one lineality space and
904         is a polyhedral fan
905NOTE:    if printlevel sufficiently high will print intermediate results
906EXAMPLE: example tropicalLinkNewton; shows an example
907"
908{
909  ring origin=basering;
910  dbprint("reducing to one-dimensional fan");
911  cone C0 = homogeneitySpace(inI);
912  intmat HH = intmat(generatorsOfLinealitySpace(C0));
913  matrix H = gauss_nf(HH);
914
915  intvec p = pivotIndices(H);
916  string resVars;
917  ideal resImage;
918  for (int i=1; i<=nvars(basering); i++)
919  {
920    if (p[i]==1)
921    {
922      resImage[i]=1;
923    }
924    else
925    {
926      resImage[i]=var(i);
927      resVars = resVars+varstr(i)+",";
928    }
929  }
930  resVars = resVars[1..size(resVars)-1];
931
932  execute("ring srestr = ("+string(char(origin))+",t),("+resVars+"),dp;");
933  number uniformizingParameter = t;
934  map resMap = origin,imap(origin,resImage);
935  ideal inI = resMap(inI);
936  inI = satstd(inI);
937
938  dbprint("intersecting with pairs of affine hyperplanes");
939  ideal substImagePos;
940  ideal substImageNeg;
941  int n = nvars(srestr);
942  list T;
943  int j;
944  for (i=1; i<=n; i++)
945  {
946    setring srestr;
947    substImagePos = maxideal(1);
948    substImagePos[i] = t;
949    substImageNeg = maxideal(1);
950    substImageNeg[i] = t^(-1);
951
952    def ssubstDP = substRing(i,"dp");
953    setring ssubstDP;
954    ideal substImagePos = imap(srestr,substImagePos);
955    map substMapPos = srestr,substImagePos;
956    ideal inIPos = substMapPos(inI);
957    if (dim(std(inIPos))<0)
958    {
959      setring srestr;
960      kill ssubstDP;
961      dbprint(string(i)+": empty");
962      i++;
963      continue;
964    }
965    dbprint(string(i)+": non-empty, computing tropical variety");
966
967    setring srestr;
968    kill ssubstDP;
969    def ssubstLP = substRing(i,"lp");
970    setring ssubstLP;
971    ideal substImagePos = imap(srestr,substImagePos);
972    map substMapPos = srestr,substImagePos;
973    ideal inIPos = std(substMapPos(inI));
974    number uniformizingParameter=t;
975    export(uniformizingParameter);
976    list Tpos = tropicalVarietyNewton(inIPos);
977    ideal substImageNeg = imap(srestr,substImageNeg);
978    map substMapNeg = srestr,substImageNeg;
979    ideal inINeg = std(substMapNeg(inI));
980    list Tneg = tropicalVarietyNewton(inINeg);
981    setring srestr;
982    list Tpos = imap(ssubstLP,Tpos);
983    list Tneg = imap(ssubstLP,Tneg);
984    kill ssubstLP;
985    T = mergeTropNewton(T,Tpos,i,number(1),p);
986    T = mergeTropNewton(T,Tneg,i,number(-1),p);
987    dbprint("total number of rays: "+string(size(T)));
988    kill Tpos;
989    kill Tneg;
990  }
991  setring origin;
992  return(imap(srestr,T));
993}
994example
995{ "EXAMPLE:"; echo = 2;
996  // a 10 valent facet in tropical Grass(3,7)
997  ring r = (0,t),
998           (p012,p013,p023,p123,p014,p024,p124,p034,p134,p234,
999            p015,p025,p125,p035,p135,p235,p045,p145,p245,p345,
1000            p016,p026,p126,p036,p136,p236,p046,p146,p246,p346,
1001            p056,p156,p256,p356,p456),
1002           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,
1003              5,7,7,5,7,4,4,4,4,4,4);
1004  number uniformizingParameter = t;
1005  export(uniformizingParameter);
1006  ideal inI =
1007    p345*p136+p134*p356,  p125*p045+p015*p245,  p124*p015-p014*p125,
1008    p135*p245-p125*p345,  p135*p045+p015*p345,  p124*p045+p014*p245,
1009    p024*p125-p012*p245,  p145*p236-p124*p356,  p124*p135-p123*p145,
1010    p024*p015+p012*p045,  p134*p026+p023*p146-p024*p136,
1011    p145*p036+p014*p356,  p014*p135-p013*p145,  p234*p145+p124*p345,
1012    p034*p145-p014*p345,  p024*p135-p012*p345,  p125*p035+p015*p235,
1013    p235*p045-p035*p245,  p234*p136-p134*p236,  p134*p036-p034*p136,
1014    p146*p356-p136*p456,  p135*p146-p134*p156,
1015    p135*p026+p023*p156+p012*p356,  p124*p035+p014*p235,
1016    p123*p025+p012*p235,  p013*p025-p012*p035,  p345*p146+p134*p456,
1017    p125*p036+p015*p236,  p345*p026-p023*p456+p024*p356,
1018    p123*p015-p013*p125,  p234*p025-p024*p235,  p034*p025-p024*p035,
1019    p234*p125+p123*p245,  p245*p036-p045*p236,  p123*p045+p013*p245,
1020    p034*p125-p013*p245,  p234*p015+p013*p245,  p245*p156+p125*p456,
1021    p034*p015+p013*p045,  p045*p156-p015*p456,  p135*p236-p123*p356,
1022    p235*p146-p134*p256,  p135*p036+p013*p356,  p124*p036+p014*p236,
1023    p123*p014-p013*p124,  p035*p146-p134*p056,  p145*p126+p124*p156,
1024    p234*p045-p034*p245,  p235*p026+p023*p256-p025*p236,
1025    p145*p016+p014*p156,  p035*p026+p023*p056-p025*p036,
1026    p345*p236+p234*p356,  p234*p135+p123*p345,  p345*p036+p034*p356,
1027    p034*p135-p013*p345,  p345*p156+p135*p456,  p124*p034+p014*p234,
1028    p145*p246-p124*p456,  p123*p024+p012*p234,  p145*p046+p014*p456,
1029    p013*p024-p012*p034,  p024*p156+p012*p456,  p125*p056+p015*p256,
1030    p245*p056-p045*p256,  p236*p146-p136*p246,  p134*p126+p123*p146,
1031    p136*p046-p036*p146,  p235*p036-p035*p236,  p134*p016+p013*p146,
1032    p123*p035+p013*p235,  p235*p156-p135*p256,
1033    p123*p026-p023*p126+p012*p236,  p135*p056-p035*p156,
1034    p023*p016-p013*p026+p012*p036,  p124*p056+p014*p256,
1035    p234*p146-p134*p246,  p025*p126-p012*p256,  p134*p046-p034*p146,
1036    p025*p016+p012*p056,  p234*p035-p034*p235,  p345*p256+p235*p456,
1037    p234*p026+p023*p246-p024*p236,  p345*p056+p035*p456,
1038    p034*p026+p023*p046-p024*p036,  p125*p016-p015*p126,
1039    p025*p246-p024*p256,  p025*p046-p024*p056,  p245*p126-p125*p246,
1040    p125*p046+p015*p246,  p045*p126+p015*p246,  p245*p016-p015*p246,
1041    p045*p016-p015*p046,  p123*p036+p013*p236,  p236*p156+p126*p356,
1042    p135*p126+p123*p156,  p036*p156-p016*p356,  p135*p016+p013*p156,
1043    p124*p016-p014*p126,  p235*p056-p035*p256,  p245*p046-p045*p246,
1044    p234*p036-p034*p236,  p123*p034+p013*p234,  p246*p356-p236*p456,
1045    p234*p156-p123*p456,  p135*p246-p123*p456,  p345*p126-p123*p456,
1046    p046*p356-p036*p456,  p034*p156+p013*p456,  p135*p046+p013*p456,
1047    p345*p016-p013*p456,  p124*p046+p014*p246,  p024*p126-p012*p246,
1048    p024*p016+p012*p046,  p345*p246+p234*p456,  p345*p046+p034*p456,
1049    p235*p126+p123*p256,  p236*p056-p036*p256,  p123*p056+p013*p256,
1050    p035*p126-p013*p256,  p235*p016+p013*p256,  p035*p016+p013*p056,
1051    p235*p246-p234*p256,  p234*p056-p034*p256,  p035*p246-p034*p256,
1052    p235*p046-p034*p256,  p035*p046-p034*p056,  p126*p036+p016*p236,
1053    p123*p016-p013*p126,  p234*p126+p123*p246,  p236*p046-p036*p246,
1054    p123*p046+p013*p246,  p034*p126-p013*p246,  p234*p016+p013*p246,
1055    p246*p156+p126*p456,  p034*p016+p013*p046,  p046*p156-p016*p456,
1056    p234*p046-p034*p246,  p126*p056+p016*p256,  p246*p056-p046*p256,
1057    p126*p046+p016*p246,  p024*p235*p145+p124*p025*p345,
1058    p024*p035*p145-p014*p025*p345,  p123*p145*p245-p124*p125*p345,
1059    p013*p145*p245-p014*p125*p345,  p013*p045*p145+p014*p015*p345,
1060    p024*p235*p136-p134*p025*p236,  p123*p245*p136+p134*p125*p236,
1061    p013*p245*p136+p134*p015*p236,  p034*p245*p136-p134*p045*p236,
1062    p134*p156*p356-p135*p136*p456,  p123*p145*p146-p124*p134*p156,
1063    p013*p145*p146-p014*p134*p156,  p013*p145*p026+p023*p014*p156+p012*p014*p356,
1064    p124*p025*p156+p012*p145*p256,  p012*p145*p056-p014*p025*p156,
1065    p024*p145*p256-p124*p025*p456,  p024*p145*p056+p014*p025*p456,
1066    p034*p235*p136-p134*p035*p236,  p134*p256*p356-p235*p136*p456,
1067    p134*p056*p356-p035*p136*p456,  p025*p036*p146-p024*p136*p056,
1068    p013*p125*p026-p023*p015*p126+p012*p015*p236,
1069    p123*p245*p146+p134*p125*p246,  p013*p245*p146+p134*p015*p246,
1070    p013*p245*p026-p023*p015*p246-p012*p045*p236,
1071    p013*p045*p026-p023*p015*p046-p012*p045*p036,
1072    p034*p245*p146-p134*p045*p246,  p013*p124*p026-p023*p014*p126+p012*p014*p236,
1073    p013*p145*p056-p014*p035*p156,  p024*p256*p356-p025*p236*p456,
1074    p024*p056*p356-p025*p036*p456,  p234*p256*p356-p235*p236*p456,
1075    p034*p256*p356-p035*p236*p456,  p034*p056*p356-p035*p036*p456,
1076    p012*p235*p145*p245+p124*p025*p125*p345,
1077    p012*p035*p145*p245-p014*p025*p125*p345,
1078    p012*p035*p045*p145+p014*p015*p025*p345,
1079    p012*p235*p245*p136-p134*p025*p125*p236,
1080    p012*p035*p245*p136+p134*p015*p025*p236,
1081    p024*p035*p245*p136-p134*p025*p045*p236,
1082    p014*p025*p125*p156+p012*p015*p145*p256,
1083    p012*p145*p245*p256-p124*p025*p125*p456,
1084    p012*p045*p145*p256+p014*p025*p125*p456,
1085    p012*p245*p256*p356-p025*p125*p236*p456,
1086    p012*p045*p256*p356+p015*p025*p236*p456,
1087    p012*p045*p056*p356+p015*p025*p036*p456,
1088    p123*p245*p256*p356+p125*p235*p236*p456,
1089    p013*p245*p256*p356+p015*p235*p236*p456,
1090    p013*p045*p256*p356+p015*p035*p236*p456,
1091    p013*p045*p056*p356+p015*p035*p036*p456;
1092  system("--random",1337);
1093  printlevel = 3;
1094  list TinI = tropicalLinkNewton(inI);
1095  for (int i=1; i<=size(TinI); i++)
1096  { print(TinI[i]); }
1097}
1098
1099
1100// disabled routines to check characteristic-freeness of tropical points
1101
1102// proc saturateWithRespectToVariable(ideal I, int k)
1103// {
1104//   ASSUME(1,k>=1);
1105//   ASSUME(1,k<=nvars(basering));
1106
1107//   def origin = basering;
1108//   int n = nvars(basering);
1109//   intvec weightVector = ringlist(origin)[3][1][2];
1110
1111//   string newVars;
1112//   for (int i=1; i<k; i++)
1113//   {
1114//     newVars = newVars+string(var(i))+",";
1115//   }
1116//   for (i=k+1; i<=n; i++)
1117//   {
1118//     newVars = newVars+string(var(i))+",";
1119//   }
1120//   newVars = newVars+string(var(k));
1121//   execute("ring ringForSaturation = ("+charstr(origin)+"),("+newVars+"),dp;");
1122
1123//   ideal I = satstd(imap(origin,I));
1124//   if (I==-1)
1125//   {
1126//     return (-1);
1127//   }
1128//   I = simplify(I,2+4+32);
1129
1130//   setring origin;
1131//   I = imap(ringForSaturation,I);
1132//   return (I);
1133// }
1134
1135// proc stepwiseSaturation(ideal I)
1136// {
1137//   if (I!=1)
1138//   {
1139//     list variablesToBeSaturated;
1140//     int l = nvars(basering);
1141//     for (int i=1; i<=l; i++)
1142//     { variablesToBeSaturated[i]=l-i+1; }
1143
1144//     while (size(variablesToBeSaturated)>0)
1145//     {
1146//       dbprint("variablesToBeSaturated: "+string(variablesToBeSaturated));
1147//       I = saturateWithRespectToVariable(I,variablesToBeSaturated[1]);
1148//       variablesToBeSaturated = delete(variablesToBeSaturated,1);
1149//       if ((I==1) || (I==-1))
1150//       {
1151//         break;
1152//       }
1153//     }
1154//   }
1155
1156//   return (I);
1157// }
1158
1159
1160// proc checkForContainmentInTropicalVariety(ideal I, intvec w, int charInt)
1161// {
1162//   def origin = basering;
1163//   intvec wOne = oneVector(nvars(origin));
1164//   execute("ring rInitialIdeal = ("+string(charInt)+"),("+varstr(origin)+"),(a(wOne),wp(w));");
1165
1166//   ideal I = imap(origin,I);
1167//   int tinI = timer;
1168//   option(redSB);
1169//   ideal stdI = satstd(I);
1170//   attrib(stdI,"isSB",1);
1171//   ideal inI = initial(stdI,w);
1172//   tinI = timer-tinI;
1173//   dbprint("time used computing initial ideal: "+string(tinI));
1174
1175//   int tsat = timer;
1176//   ideal satinI = stepwiseSaturation(inI);
1177//   tsat = timer-tsat;
1178//   dbprint("time used computing saturation: "+string(tsat));
1179
1180
1181//   export(I);
1182//   export(stdI);
1183//   export(inI);
1184//   export(satinI);
1185//   return (rInitialIdeal);
1186// }
Note: See TracBrowser for help on using the repository browser.