1 | /////////////////////////////////////////////////////////////////////////////// |
---|
2 | version="version tropical.lib 4.0.3.3 Nov_2016 "; // $Id$ |
---|
3 | category="Tropical Geometry"; |
---|
4 | info=" |
---|
5 | LIBRARY: tropicalNewton.lib Computations in Tropical Geometry |
---|
6 | using Newton Polygon methods |
---|
7 | |
---|
8 | AUTHORS: Tommy Hofman, email: thofmann@mathematik.uni.kl.de |
---|
9 | Yue Ren, email: reny@cs.bgu.ac.il |
---|
10 | |
---|
11 | OVERVIEW: |
---|
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 | |
---|
18 | REFERENCES: 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 | |
---|
22 | SEE ALSO: tropicalVariety, tropical_lib |
---|
23 | |
---|
24 | KEYWORDS: tropical geometry; tropical varieties; Newton polygons |
---|
25 | |
---|
26 | PROCEDURES: |
---|
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 | |
---|
42 | LIB "gfan.lib"; |
---|
43 | LIB "triang.lib"; |
---|
44 | LIB "linalg.lib"; |
---|
45 | |
---|
46 | /////////////////////////////////////////////////////////////////////////////// |
---|
47 | |
---|
48 | proc setUniformizingParameter(number p) |
---|
49 | "USAGE: setUniformizingParameter(p); p number |
---|
50 | RETURN: none, sets the uniformizing parameter as p |
---|
51 | ASSUME: char(K)==0 and p prime number or |
---|
52 | trdeg(K)>0 and p transcendental variable or |
---|
53 | p==0 |
---|
54 | EXAMPLE: 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 | } |
---|
86 | example |
---|
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 | |
---|
100 | proc val(number c) |
---|
101 | "USAGE: val(c); c number |
---|
102 | RETURN: int, the valuation of a element in the ground field |
---|
103 | ASSUME: uniformizingParameter is set and c!=0 |
---|
104 | EXAMPLE: 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 | } |
---|
153 | example |
---|
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 | |
---|
167 | proc newtonPolygonNegSlopes(poly g, list #) |
---|
168 | "USAGE: newtonPolygonNegSlopes(g,b); g poly, b int |
---|
169 | RETURN: list, the negative slopes of the Newton Polygon of g |
---|
170 | if b==1, computes root (type poly) instead if (easily) possible |
---|
171 | ASSUME: uniformizingParameter is set and g univariate |
---|
172 | EXAMPLE: 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 | } |
---|
217 | example |
---|
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 | |
---|
231 | static 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 | |
---|
252 | static 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 | |
---|
265 | static 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 | |
---|
274 | static 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 | |
---|
302 | static 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 | |
---|
359 | static 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 | |
---|
376 | static 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 | |
---|
395 | proc cccMatrixToPositiveIntvec(matrix L) |
---|
396 | "USAGE: cccMatrixToPositiveIntvec(M); M matrix |
---|
397 | RETURN: intvec, strictly positive equivalent as weight vector to row vector in M |
---|
398 | ASSUME: constant coefficient case only, |
---|
399 | will scale weight vector and add vectors of ones to it |
---|
400 | EXAMPLE: 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 | } |
---|
422 | example |
---|
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 | |
---|
442 | static 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 | |
---|
452 | static 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 | |
---|
466 | proc tropicalPointNewton(ideal I, list #) |
---|
467 | "USAGE: tropicalPointsLasVegas(I); I ideal |
---|
468 | RETURN: matrix, a matrix containing a tropical point as row vector |
---|
469 | ASSUME: uniformizingParameter is set and I monomial free |
---|
470 | NOTE: if printlevel sufficiently high will print intermediate data and timings |
---|
471 | returns error if randomly chosen hyperplanes are not generic |
---|
472 | EXAMPLE: 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 | } |
---|
558 | example |
---|
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 | |
---|
580 | static 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 | |
---|
592 | static 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 | |
---|
625 | proc tropicalVarietyNewton(ideal I) |
---|
626 | "USAGE: tropicalVarietyNewton(I); I ideal |
---|
627 | RETURN: matrix, a matrix containing all elements of the tropical variety |
---|
628 | ASSUME: uniformizingParameter is set, I monomial free and zero-dimensional |
---|
629 | EXAMPLE: 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 | } |
---|
716 | example |
---|
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 | |
---|
729 | proc switchRingsAndComputeInitialIdeal(ideal I, intvec w) |
---|
730 | "USAGE: switchRingsAndComputeInitialIdeal(I,w); I ideal, w intvec |
---|
731 | RETURN: ring, a ring containing the initial ideal with respect to w |
---|
732 | ASSUME: constant coefficient case and w strictly positive integer |
---|
733 | NOTE: if printlevel sufficiently high will print timing |
---|
734 | EXAMPLE: 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 | } |
---|
753 | example |
---|
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 | |
---|
774 | static 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 | |
---|
799 | static 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 | |
---|
811 | static 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 | |
---|
834 | static 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 | |
---|
876 | static 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 | |
---|
899 | proc tropicalLinkNewton(ideal inI) |
---|
900 | "USAGE: tropicalLinkNewton(inI); inI ideal |
---|
901 | RETURN: matrix, a matrix containing generators of all rays of the tropical variety |
---|
902 | ASSUME: constant coefficient case, inI is monomial free, |
---|
903 | its tropical variety has codimension one lineality space and |
---|
904 | is a polyhedral fan |
---|
905 | NOTE: if printlevel sufficiently high will print intermediate results |
---|
906 | EXAMPLE: 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 | } |
---|
994 | example |
---|
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 | // } |
---|