Changeset 3360fb in git
- Timestamp:
- Apr 14, 2009, 2:00:15 PM (14 years ago)
- Branches:
- (u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
- Children:
- 334c21f9e573112667d39f4cd9dcad9c7b19a50c
- Parents:
- 2c3a5db391a25c28449b80c5f27e1dd77ed2c30f
- Location:
- Singular/LIB
- Files:
-
- 18 edited
-
algebra.lib (modified) (3 diffs)
-
bfun.lib (modified) (28 diffs)
-
crypto.lib (modified) (2 diffs)
-
decodegb.lib (modified) (2 diffs)
-
dmod.lib (modified) (31 diffs)
-
elim.lib (modified) (6 diffs)
-
freegb.lib (modified) (5 diffs)
-
jacobson.lib (modified) (7 diffs)
-
nctools.lib (modified) (2 diffs)
-
polymake.lib (modified) (94 diffs)
-
presolve.lib (modified) (16 diffs)
-
random.lib (modified) (3 diffs)
-
ratgb.lib (modified) (5 diffs)
-
sing4ti2.lib (modified) (12 diffs)
-
teachstd.lib (modified) (2 diffs)
-
triang.lib (modified) (2 diffs)
-
tropical.lib (modified) (324 diffs)
-
zeroset.lib (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/algebra.lib
r2c3a5d r3360fb 3 3 //new proc nonZeroEntry(id), used to fix a bug in proc finitenessTest 4 4 /////////////////////////////////////////////////////////////////////////////// 5 version="$Id: algebra.lib,v 1.2 1 2009-04-07 16:18:05 seelischExp $";5 version="$Id: algebra.lib,v 1.22 2009-04-14 12:00:13 Singular Exp $"; 6 6 category="Commutative Algebra"; 7 7 info=" … … 254 254 NOTE: the proc algebra_containment tests the same using a different 255 255 algorithm, which is often faster 256 if l[1] == 0 then l[2] may contain more than one relation h(y(0),y(1),...,y(k)), 256 if l[1] == 0 then l[2] may contain more than one relation h(y(0),y(1),...,y(k)), 257 257 separated by comma 258 258 EXAMPLE: example inSubring; shows an example … … 300 300 } 301 301 if (deg(ker[i]/y(0))>0) 302 { if( l[2] != "" ){ l[2] = l[2] + ","; } 302 { if( l[2] != "" ){ l[2] = l[2] + ","; } 303 303 l[2] = l[2] + string(ker[i]); 304 304 } -
Singular/LIB/bfun.lib
r2c3a5d r3360fb 1 1 ////////////////////////////////////////////////////////////////////////////// 2 version="$Id: bfun.lib,v 1. 7 2009-04-08 16:51:07 seelischExp $";2 version="$Id: bfun.lib,v 1.8 2009-04-14 12:00:13 Singular Exp $"; 3 3 category="Noncommutative"; 4 4 info=" … … 10 10 @* one is interested in the global b-function (also known as Bernstein-Sato 11 11 @* polynomial) b(s) in K[s], defined to be the non-zero monic polynomial of minimal 12 @* degree, satisfying a functional identity L * F^{s+1} = b(s) F^s, 12 @* degree, satisfying a functional identity L * F^{s+1} = b(s) F^s, 13 13 @* for some operator L in D[s] (* stands for the action of differential operator) 14 14 @* By D one denotes the n-th Weyl algebra … … 21 21 @* There is a constructive definition of a b-function of a holonomic ideal I in D 22 22 @* (that is, an ideal I in a Weyl algebra D, such that D/I is holonomic module) 23 @* with respect to the given weight vector w: For a poly p in D, its initial 24 @* form w.r.t. (-w,w) is defined as the sum of all terms of p which have 23 @* with respect to the given weight vector w: For a poly p in D, its initial 24 @* form w.r.t. (-w,w) is defined as the sum of all terms of p which have 25 25 @* maximal weighted total degree where the weight of x_i is -w_i and the weight 26 @* of d_i is w_i. Let J be the initial ideal of I w.r.t. (-w,w), i.e. the 26 @* of d_i is w_i. Let J be the initial ideal of I w.r.t. (-w,w), i.e. the 27 27 @* K-vector space generated by all initial forms w.r.t (-w,w) of elements of I. 28 28 @* Put s = w_1 x_1 d_1 + ... + w_n x_n d_n. Then the monic generator b_w(s) of 29 @* the intersection of J with the PID K[s] is called the b-function of I w.r.t. w. 29 @* the intersection of J with the PID K[s] is called the b-function of I w.r.t. w. 30 30 @* Unlike Bernstein-Sato polynomial, general b-function with respect to 31 31 @* arbitrary weights need not have rational roots at all. However, b-function 32 32 @* of a holonomic ideal is known to be non-zero as well. 33 @* 34 @* References: [SST] Saito, Sturmfels, Takayama: Groebner Deformations of 33 @* 34 @* References: [SST] Saito, Sturmfels, Takayama: Groebner Deformations of 35 35 @* Hypergeometric Differential Equations (2000), 36 36 @* Noro: An Efficient Modular Algorithm for Computing the Global b-function, … … 47 47 pIntersect(f,I[,s]); intersection of ideal with subalgebra K[f] for a poly f 48 48 pIntersectSyz(f,I[,p,s,t]); intersection of ideal with subalgebra K[f] with syz-solver 49 linReduce(f,I[,s]); reduce a poly by linear reductions wrt ideal 49 linReduce(f,I[,s]); reduce a poly by linear reductions wrt ideal 50 50 linReduceIdeal(I,[s,t]); interreduce generators of ideal by linear reductions 51 51 linSyzSolve(I[,s]); compute a linear dependency of elements of ideal I … … 53 53 AUXILIARY PROCEDURES: 54 54 55 allPositive(v); checks whether all entries of an intvec are positive 55 allPositive(v); checks whether all entries of an intvec are positive 56 56 scalarProd(v,w); computes the standard scalar product of two intvecs 57 57 vec2poly(v[,i]); constructs an univariate poly with given coefficients … … 97 97 EXAMPLE: example gradedWeyl; shows examples 98 98 NOTE: u[i] is the weight of x(i), v[i] the weight of D(i). 99 @* u+v has to be a non-negative intvec. 99 @* u+v has to be a non-negative intvec. 100 100 " 101 101 { … … 107 107 { 108 108 ERROR("weight vectors have wrong dimension"); 109 } 109 } 110 110 intvec uv,gr; 111 111 uv = u+v; … … 133 133 for (i=1; i<=n; i++) 134 134 { 135 if (gr[i] == 1) 135 if (gr[i] == 1) 136 136 { 137 137 l2[n+i] = "xi("+string(i)+")"; … … 358 358 for (i=1; i<=ncols(I); i++) { M[i] = gen(i); } 359 359 } 360 } 360 } 361 361 dbprint(ppl-1,"// initially sorted ideal:", I); 362 362 if (remembercoeffs <> 0) { dbprint(ppl-1,"// used permutations:", M); } … … 411 411 lmJ = insertGenerator(lmJ,lm,j); 412 412 ordJ = insertGenerator(ordJ,poly(ordlm),j); 413 if (remembercoeffs <> 0) 413 if (remembercoeffs <> 0) 414 414 { 415 415 v = M[i]; … … 434 434 { 435 435 if (ordJ[i] == ord(J[j][k])) 436 { 436 { 437 437 if (lm == normalize(J[j][k])) 438 438 { … … 469 469 RETURN: poly or list, linear reductum (over field) of f by elements from I 470 470 PURPOSE: reduce a poly only by linear reductions (no monomial multiplications) 471 NOTE: If s<>0, a list consisting of the reduced poly and the coefficient 471 NOTE: If s<>0, a list consisting of the reduced poly and the coefficient 472 472 @* vector of the used reductions is returned, otherwise (and by default) 473 473 @* only reduced poly is returned. 474 474 @* If t<>0 (and by default) all monomials are reduced (if possible), 475 475 @* otherwise, only leading monomials are reduced. 476 @* If u<>0 (and by default), the ideal is linearly pre-reduced, i.e. 476 @* If u<>0 (and by default), the ideal is linearly pre-reduced, i.e. 477 477 @* instead of the given ideal, the output of @code{linReduceIdeal} is used. 478 478 @* If u is set to 0 and the given ideal does not equal the output of … … 605 605 f = x3+y2+x2+y+x; 606 606 I = x3-y3, y3-x2,x3-y2,x2-y,y2-x; 607 list l = linReduce(f,I,1); 607 list l = linReduce(f,I,1); 608 608 l; 609 609 module M = I; … … 739 739 proc pIntersect (poly s, ideal I, list #) 740 740 "USAGE: pIntersect(f, I [,s]); f a poly, I an ideal, s an optional int 741 RETURN: vector, coefficient vector of the monic polynomial 741 RETURN: vector, coefficient vector of the monic polynomial 742 742 PURPOSE: compute the intersection of ideal I with the subalgebra K[f] 743 743 ASSUME: I is given as Groebner basis, basering is not a qring. … … 810 810 print("// Try a bound of at least " + string(degbound)); 811 811 return(vector(0)); 812 } 812 } 813 813 dbprint(ppl,"// lower bound for the degree of the insection is " +string(degbound)); 814 814 if (degbound == 0) // lm(s) does not appear in lm(I) … … 894 894 v = v + m[j,1]*gen(j); 895 895 } 896 setring save; 896 setring save; 897 897 v = imap(@R,v); 898 898 kill @R; … … 918 918 NOTE: If the intersection is zero, this procedure might not terminate. 919 919 @* If p>0 is given, this proc computes the generator of the intersection in 920 @* char p first and then only searches for a generator of the obtained 921 @* degree in the basering. Otherwise, it searches for all degrees by 920 @* char p first and then only searches for a generator of the obtained 921 @* degree in the basering. Otherwise, it searches for all degrees by 922 922 @* computing syzygies. 923 923 @* If s<>0, @code{std} is used for Groebner basis computations in char 0, 924 924 @* otherwise, and by default, @code{slimgb} is used. 925 @* If t<>0 and by default, @code{std} is used for Groebner basis 925 @* If t<>0 and by default, @code{std} is used for Groebner basis 926 926 @* computations in char >0, otherwise, @code{slimgb} is used. 927 927 DISPLAY: If printlevel=1, progress debug messages will be printed, … … 1131 1131 } 1132 1132 dbprint(ppl,"// pIntersectSyz finished"); 1133 if (solveincharp) { short = shortSave; } 1133 if (solveincharp) { short = shortSave; } 1134 1134 return(v); 1135 1135 } … … 1272 1272 { 1273 1273 return(list(ideal(0),intvec(0))); 1274 } 1274 } 1275 1275 if (inorann == 0) // bfct using initial ideal 1276 1276 { … … 1364 1364 RETURN: list of ideal and intvec 1365 1365 PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) 1366 @* for the hypersurface defined by f. 1366 @* for the hypersurface defined by f. 1367 1367 ASSUME: The basering is commutative and of characteristic 0. 1368 1368 BACKGROUND: In this proc, the initial Malgrange ideal is computed according to 1369 1369 @* the algorithm by Masayuki Noro and then a system of linear equations is 1370 1370 @* solved by linear reductions. 1371 NOTE: In the output list, the ideal contains all the roots 1371 NOTE: In the output list, the ideal contains all the roots 1372 1372 @* and the intvec their multiplicities. 1373 1373 @* If s<>0, @code{std} is used for GB computations, 1374 @* otherwise, and by default, @code{slimgb} is used. 1374 @* otherwise, and by default, @code{slimgb} is used. 1375 1375 @* If t<>0, a matrix ordering is used for Groebner basis computations, 1376 1376 @* otherwise, and by default, a block ordering is used. … … 1435 1435 @* the algorithm by Masayuki Noro and then a system of linear equations is 1436 1436 @* solved by computing syzygies. 1437 NOTE: In the output list, the ideal contains all the roots and the intvec 1437 NOTE: In the output list, the ideal contains all the roots and the intvec 1438 1438 @* their multiplicities. 1439 1439 @* If r<>0, @code{std} is used for GB computations in characteristic 0, … … 1520 1520 ASSUME: The basering is the n-th Weyl algebra in characteristic 0 and for all 1521 1521 @* 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the 1522 @* sequence of variables is given by x(1),...,x(n),D(1),...,D(n), 1522 @* sequence of variables is given by x(1),...,x(n),D(1),...,D(n), 1523 1523 @* where D(i) is the differential operator belonging to x(i). 1524 1524 @* Further we assume that I is holonomic. … … 1533 1533 @* L[3] is 1 (for nonzero constant) or 0 (for zero b-function). 1534 1534 @* If s<>0, @code{std} is used for GB computations in characteristic 0, 1535 @* otherwise, and by default, @code{slimgb} is used. 1535 @* otherwise, and by default, @code{slimgb} is used. 1536 1536 @* If t<>0, a matrix ordering is used for GB computations, otherwise, 1537 1537 @* and by default, a block ordering is used. … … 1622 1622 @* their multiplicities. 1623 1623 @* If s<>0, @code{std} is used for the GB computation, otherwise, 1624 @* and by default, @code{slimgb} is used. 1624 @* and by default, @code{slimgb} is used. 1625 1625 @* If t<>0, a matrix ordering is used for GB computations, 1626 1626 @* otherwise, and by default, a block ordering is used. … … 1634 1634 def save = basering; 1635 1635 int n = nvars(save); 1636 if (char(save) <> 0) 1636 if (char(save) <> 0) 1637 1637 { 1638 1638 ERROR("characteristic of basering has to be 0"); … … 1671 1671 // create names for vars 1672 1672 list Lvar; 1673 Lvar[1] = safeVarName("t"); 1673 Lvar[1] = safeVarName("t"); 1674 1674 Lvar[2] = safeVarName("s"); 1675 1675 Lvar[n+3] = safeVarName("D"+Lvar[1]); … … 1684 1684 intvec @a = 1:N; @a[2] = 2; 1685 1685 intvec @a2 = @a; @a2[2] = 0; @a2[2*n+4] = 0; 1686 list Lord; 1686 list Lord; 1687 1687 Lord[1] = list("a",@a); Lord[2] = list("a",@a2); 1688 1688 if (methodord == 0) // default: block ordering … … 1782 1782 @* their multiplicities. 1783 1783 @* If r<>0, @code{std} is used for GB computations, 1784 @* otherwise, and by default, @code{slimgb} is used. 1784 @* otherwise, and by default, @code{slimgb} is used. 1785 1785 DISPLAY: If printlevel=1, progress debug messages will be printed, 1786 1786 @* if printlevel>=2, all the debug messages will be printed. -
Singular/LIB/crypto.lib
r2c3a5d r3360fb 1 1 //GP, last modified 28.6.06 2 2 /////////////////////////////////////////////////////////////////////////////// 3 version="$Id: crypto.lib,v 1. 7 2009-04-06 12:39:02 seelischExp $";3 version="$Id: crypto.lib,v 1.8 2009-04-14 12:00:14 Singular Exp $"; 4 4 category="Teaching"; 5 5 info=" … … 1575 1575 "USAGE: generateG(a,b,m); 1576 1576 RETURN: m-th division polynomial 1577 NOTE: generate the so-called division polynomials, i.e., the recursively defined 1577 NOTE: generate the so-called division polynomials, i.e., the recursively defined 1578 1578 polynomials p_m=generateG(a,b,m) in Z[x, y] such that, for a point (x:y:1) on the 1579 1579 elliptic curve defined by y^2=x^3+a*x+b over Z/N the point@* -
Singular/LIB/decodegb.lib
r2c3a5d r3360fb 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: decodegb.lib,v 1.1 3 2009-04-10 15:18:18Singular Exp $";2 version="$Id: decodegb.lib,v 1.14 2009-04-14 12:00:14 Singular Exp $"; 3 3 category="Coding theory"; 4 4 info=" … … 1995 1995 "USAGE: decodeRandomFL(redun,p,e,n,t,ncodes,ntrials,minpol); 1996 1996 @format 1997 - n is length of codes generated, 1997 - n is length of codes generated, 1998 1998 - redun = redundancy of codes generated, 1999 1999 - p is the characteristic, -
Singular/LIB/dmod.lib
r2c3a5d r3360fb 1 1 ////////////////////////////////////////////////////////////////////////////// 2 version="$Id: dmod.lib,v 1. 39 2009-04-09 12:04:41 seelischExp $";2 version="$Id: dmod.lib,v 1.40 2009-04-14 12:00:14 Singular Exp $"; 3 3 category="Noncommutative"; 4 4 info=" … … 7 7 @* Jorge Martin Morales, jorge@unizar.es 8 8 9 THEORY: Let K be a field of characteristic 0. Given a polynomial ring 9 THEORY: Let K be a field of characteristic 0. Given a polynomial ring 10 10 @* R = K[x_1,...,x_n] and a polynomial F in R, 11 @* one is interested in the R[1/F]-module of rank one, generated by 11 @* one is interested in the R[1/F]-module of rank one, generated by 12 12 @* the symbol F^s for a symbolic discrete variable s. 13 13 @* In fact, the module R[1/F]*F^s has a structure of a D(R)[s]-module, where D(R) … … 19 19 @* One is interested in the following data: 20 20 @* - Ann F^s = I = I(F^s) in D(R)[s], denoted by LD in the output 21 @* - global Bernstein polynomial in K[s], denoted by bs, 21 @* - global Bernstein polynomial in K[s], denoted by bs, 22 22 @* - its minimal integer root s0, the list of all roots of bs, which are known 23 23 @* to be rational, with their multiplicities, which is denoted by BS 24 @* - Ann F^s0 = I(F^s0) in D(R), denoted by LD0 in the output 24 @* - Ann F^s0 = I(F^s0) in D(R), denoted by LD0 in the output 25 25 @* (LD0 is a holonomic ideal in D(R)) 26 26 @* - Ann^(1) F^s in D(R)[s], denoted by LD1 (logarithmic derivations) … … 28 28 @* PS*F^(s+1) = bs*F^s holds in K[x_1,...,x_n,1/F]*F^s. 29 29 30 REFERENCES: 30 REFERENCES: 31 31 @* We provide the following implementations of algorithms: 32 @* (OT) the classical Ann F^s algorithm from Oaku and Takayama (Journal of 32 @* (OT) the classical Ann F^s algorithm from Oaku and Takayama (Journal of 33 33 @* Pure and Applied Math., 1999), 34 34 @* (LOT) Levandovskyy's modification of the Oaku-Takayama algorithm (ISSAC 2007) … … 36 36 @* l'ideal de Bernstein associe a des polynomes, preprint, 2002) 37 37 @* (LM08) V. Levandovskyy and J. Martin-Morales, ISSAC 2008 38 @* (C) Countinho, A Primer of Algebraic D-Modules, 39 @* (SST) Saito, Sturmfels, Takayama 'Groebner Deformations of Hypergeometric 38 @* (C) Countinho, A Primer of Algebraic D-Modules, 39 @* (SST) Saito, Sturmfels, Takayama 'Groebner Deformations of Hypergeometric 40 40 @* Differential Equations', Springer, 2000 41 41 … … 143 143 "USAGE: annfs(f [,S,eng]); f a poly, S a string, eng an optional int 144 144 RETURN: ring 145 PURPOSE: compute the D-module structure of basering[1/f]*f^s with the algorithm 145 PURPOSE: compute the D-module structure of basering[1/f]*f^s with the algorithm 146 146 @* given in S and with the Groebner basis engine given in ''eng'' 147 147 NOTE: activate the output ring with the @code{setring} command. … … 267 267 "USAGE: Sannfs(f [,S,eng]); f a poly, S a string, eng an optional int 268 268 RETURN: ring 269 PURPOSE: compute the D-module structure of basering[f^s] with the algorithm 269 PURPOSE: compute the D-module structure of basering[f^s] with the algorithm 270 270 @* given in S and with the Groebner basis engine given in eng 271 271 NOTE: activate the output ring with the @code{setring} command. … … 388 388 PURPOSE: compute the D-module structure of basering[1/f]*f^s 389 389 NOTE: activate the output ring with the @code{setring} command. 390 @* In the output ring D[s], the ideal LD1 is generated by the elements 390 @* In the output ring D[s], the ideal LD1 is generated by the elements 391 391 @* in Ann F^s in D[s], coming from logarithmic derivations. 392 392 @* If eng <>0, @code{std} is used for Groebner basis computations, … … 536 536 "USAGE: ALTannfsBM(f [,eng]); f a poly, eng an optional int 537 537 RETURN: ring 538 PURPOSE: compute the annihilator ideal of f^s in D[s], where D is the Weyl 538 PURPOSE: compute the annihilator ideal of f^s in D[s], where D is the Weyl 539 539 @* algebra, according to the algorithm by Briancon and Maisonobe 540 540 NOTE: activate the output ring with the @code{setring} command. In this ring, … … 728 728 "USAGE: bernsteinBM(f [,eng]); f a poly, eng an optional int 729 729 RETURN: list (of roots of the Bernstein polynomial b and their multiplicies) 730 PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface, 730 PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface, 731 731 @* defined by f, according to the algorithm by Briancon and Maisonobe 732 732 NOTE: If eng <>0, @code{std} is used for Groebner basis computations, … … 1251 1251 "USAGE: annfs2(I, F [,eng]); I an ideal, F a poly, eng an optional int 1252 1252 RETURN: ring 1253 PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, 1253 PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, 1254 1254 @* based on the output of Sannfs-like procedure 1255 1255 @* annfs2 uses shorter expressions in the variable s (the idea of Noro). … … 1412 1412 "USAGE: annfsRB(I, F [,eng]); I an ideal, F a poly, eng an optional int 1413 1413 RETURN: ring 1414 PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, 1414 PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, 1415 1415 @* based on the output of Sannfs like procedure 1416 1416 NOTE: activate the output ring with the @code{setring} command. In this ring, … … 1601 1601 "USAGE: operatorBM(f [,eng]); f a poly, eng an optional int 1602 1602 RETURN: ring 1603 PURPOSE: compute the B-operator and other relevant data for Ann F^s, 1603 PURPOSE: compute the B-operator and other relevant data for Ann F^s, 1604 1604 @* using e.g. algorithm by Briancon and Maisonobe for Ann F^s and BS. 1605 1605 NOTE: activate the output ring with the @code{setring} command. In the output ring D[s] … … 2238 2238 "USAGE: annfsBMI(F [,eng]); F an ideal, eng an optional int 2239 2239 RETURN: ring 2240 PURPOSE: compute the D-module structure of basering[1/f]*f^s where 2240 PURPOSE: compute the D-module structure of basering[1/f]*f^s where 2241 2241 @* f = F[1]*..*F[P], according to the algorithm by Briancon and Maisonobe. 2242 2242 NOTE: activate the output ring with the @code{setring} command. In this ring, … … 2567 2567 "USAGE: annfsOT(f [,eng]); f a poly, eng an optional int 2568 2568 RETURN: ring 2569 PURPOSE: compute the D-module structure of basering[1/f]*f^s, 2569 PURPOSE: compute the D-module structure of basering[1/f]*f^s, 2570 2570 @* according to the algorithm by Oaku and Takayama 2571 2571 NOTE: activate the output ring with the @code{setring} command. In this ring, … … 2949 2949 "USAGE: SannfsOT(f [,eng]); f a poly, eng an optional int 2950 2950 RETURN: ring 2951 PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the 2951 PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the 2952 2952 @* 1st step of the algorithm by Oaku and Takayama in the ring D[s] 2953 2953 NOTE: activate the output ring with the @code{setring} command. 2954 @* In the output ring D[s], the ideal LD (which is NOT a Groebner basis) 2954 @* In the output ring D[s], the ideal LD (which is NOT a Groebner basis) 2955 2955 @* is the needed D-module structure. 2956 2956 @* If eng <>0, @code{std} is used for Groebner basis computations, … … 3233 3233 "USAGE: SannfsBM(f [,eng]); f a poly, eng an optional int 3234 3234 RETURN: ring 3235 PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the 3235 PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the 3236 3236 @* 1st step of the algorithm by Briancon and Maisonobe in the ring D[s]. 3237 3237 NOTE: activate the output ring with the @code{setring} command. 3238 @* In the output ring D[s], the ideal LD (which is NOT a Groebner basis) is 3238 @* In the output ring D[s], the ideal LD (which is NOT a Groebner basis) is 3239 3239 @* the needed D-module structure. 3240 3240 @* If eng <>0, @code{std} is used for Groebner basis computations, … … 3442 3442 PURPOSE: compute Ann f^s and Groebner basis of Ann f^s+f in D[s] 3443 3443 NOTE: activate the output ring with the @code{setring} command. 3444 @* This procedure, unlike SannfsBM, returns a ring with the degrevlex 3444 @* This procedure, unlike SannfsBM, returns a ring with the degrevlex 3445 3445 @* ordering in all variables. 3446 3446 @* In the ring D[s], the ideal LD is the ideal needed (which is returned as a Groebner basis). … … 3648 3648 PURPOSE: compute Ann f^s and Groebner basis of Ann f^s+f in D[s] 3649 3649 NOTE: activate the output ring with the @code{setring} command. 3650 @* This procedure, unlike SannfsBM, returns a ring with the degrevlex 3650 @* This procedure, unlike SannfsBM, returns a ring with the degrevlex 3651 3651 @* ordering in all variables. 3652 3652 @* In the ring D[s], the ideal LD (which IS a Groebner basis) is the needed ideal. … … 3657 3657 " 3658 3658 { 3659 // DEBUG INFO: ordering on the output ring = dp, 3659 // DEBUG INFO: ordering on the output ring = dp, 3660 3660 // use std(K,F); for reusing the std property of K 3661 3661 … … 3869 3869 "USAGE: SannfsLOT(f [,eng]); f a poly, eng an optional int 3870 3870 RETURN: ring 3871 PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the 3871 PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the 3872 3872 @* Levandovskyy's modification of the algorithm by Oaku and Takayama in D[s] 3873 3873 NOTE: activate the output ring with the @code{setring} command. 3874 @* In the ring D[s], the ideal LD (which is NOT a Groebner basis) is 3874 @* In the ring D[s], the ideal LD (which is NOT a Groebner basis) is 3875 3875 @* the needed D-module structure. 3876 3876 @* If eng <>0, @code{std} is used for Groebner basis computations, … … 4326 4326 "USAGE: annfsLOT(F [,eng]); F a poly, eng an optional int 4327 4327 RETURN: ring 4328 PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to 4328 PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to 4329 4329 @* the Levandovskyy's modification of the algorithm by Oaku and Takayama 4330 4330 NOTE: activate the output ring with the @code{setring} command. In this ring, … … 4371 4371 "USAGE: annfs0(I, F [,eng]); I an ideal, F a poly, eng an optional int 4372 4372 RETURN: ring 4373 PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, based 4373 PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, based 4374 4374 @* on the output of Sannfs-like procedure 4375 4375 NOTE: activate the output ring with the @code{setring} command. In this ring, … … 4860 4860 "USAGE: annfspecial(I,F,mir,n); I an ideal, F a poly, int mir, number n 4861 4861 RETURN: ideal 4862 PURPOSE: compute the annihilator ideal of F^n in the Weyl Algebra 4862 PURPOSE: compute the annihilator ideal of F^n in the Weyl Algebra 4863 4863 @* for the given rational number n 4864 4864 ASSUME: The basering is D[s] and contains 's' explicitly as a variable, … … 4866 4866 @* the integer 'mir' is the minimal integer root of the BS polynomial of F, 4867 4867 @* and the number n is rational. 4868 NOTE: We compute the real annihilator for any rational value of n (both 4869 @* generic and exceptional). The implementation goes along the lines of 4870 @* the Algorithm 5.3.15 from Saito-Sturmfels-Takayama. 4868 NOTE: We compute the real annihilator for any rational value of n (both 4869 @* generic and exceptional). The implementation goes along the lines of 4870 @* the Algorithm 5.3.15 from Saito-Sturmfels-Takayama. 4871 4871 DISPLAY: If printlevel=1, progress debug messages will be printed, 4872 4872 @* if printlevel>=2, all the debug messages will be printed. … … 5192 5192 RETURN: int 5193 5193 ASSUME: Basering is a commutative ring, alpha is a rational number. 5194 PURPOSE: check whether a rational number alpha is a root of the global 5194 PURPOSE: check whether a rational number alpha is a root of the global 5195 5195 @* Bernstein-Sato polynomial of f and compute its multiplicity, 5196 5196 @* with the algorithm given by S and with the Groebner basis engine given by eng. 5197 NOTE: The annihilator of f^s in D[s] is needed and hence it is computed with the 5197 NOTE: The annihilator of f^s in D[s] is needed and hence it is computed with the 5198 5198 @* algorithm by Briancon and Maisonobe. The value of a string S can be 5199 5199 @* 'alg1' (default) - for the algorithm 1 of [LM08] … … 5442 5442 proc checkRoot2 (ideal I, poly F, number a, list #) 5443 5443 "USAGE: checkRoot2(I,f,a [,eng]); I an ideal, f a poly, alpha a number, eng an optional int 5444 ASSUME: I is the annihilator of f^s in D[s], basering is D[s], 5444 ASSUME: I is the annihilator of f^s in D[s], basering is D[s], 5445 5445 @* that is basering and I are the output os Sannfs-like procedure, 5446 5446 @* f is a polynomial in K[_x] and alpha is a rational number. 5447 RETURN: int, the multiplicity of -alpha as a root of the BS polynomial of f. 5447 RETURN: int, the multiplicity of -alpha as a root of the BS polynomial of f. 5448 5448 PURPOSE: check whether a rational number alpha is a root of the global Bernstein- 5449 5449 @* Sato polynomial of f and compute its multiplicity from the known Ann F^s in D[s] … … 5457 5457 { 5458 5458 5459 5459 5460 5460 // to check: alpha is rational (has char 0 check inside) 5461 5461 if (!isRational(a)) … … 5592 5592 ASSUME: checkFactor is called from the basering, created by Sannfs-like proc, 5593 5593 @* that is, from the Weyl algebra in x1,..,xN,d1,..,dN tensored with K[s]. 5594 @* The ideal I is the annihilator of f^s in D[s], that is the ideal, computed 5594 @* The ideal I is the annihilator of f^s in D[s], that is the ideal, computed 5595 5595 @* by Sannfs-like procedure (usually called LD there). 5596 5596 @* Moreover, f is a polynomial in K[x1,..,xN] and qs is a polynomial in K[s]. 5597 5597 RETURN: int, 1 if qs is a factor of the global Bernstein polynomial of f and 0 otherwise 5598 PURPOSE: check whether a univariate polynomial qs is a factor of the 5598 PURPOSE: check whether a univariate polynomial qs is a factor of the 5599 5599 @* Bernstein-Sato polynomial of f without explicit knowledge of the latter. 5600 5600 NOTE: If eng <>0, @code{std} is used for Groebner basis computations, … … 5693 5693 "USAGE: indAR(L,n); list L, int n 5694 5694 RETURN: list 5695 PURPOSE: computes arrangement inductively, using L and 5695 PURPOSE: computes arrangement inductively, using L and 5696 5696 @* var(n) as the next variable 5697 5697 ASSUME: L has a structure of an arrangement … … 5741 5741 5742 5742 proc isRational(number n) 5743 "USAGE: isRational(n); n number 5743 "USAGE: isRational(n); n number 5744 5744 RETURN: int 5745 5745 PURPOSE: determine whether n is a rational number, -
Singular/LIB/elim.lib
r2c3a5d r3360fb 1 // $Id: elim.lib,v 1. 29 2009-04-07 16:18:05 seelischExp $1 // $Id: elim.lib,v 1.30 2009-04-14 12:00:14 Singular Exp $ 2 2 // (GMG, modified 22.06.96) 3 3 // GMG, last modified 30.10.08: new procedure elimRing; … … 10 10 // and can now choose as method slimgb or std. 11 11 /////////////////////////////////////////////////////////////////////////////// 12 version="$Id: elim.lib,v 1. 29 2009-04-07 16:18:05 seelischExp $";12 version="$Id: elim.lib,v 1.30 2009-04-14 12:00:14 Singular Exp $"; 13 13 category="Commutative Algebra"; 14 14 info=" … … 223 223 //------------------ set resp. compute ring weights ---------------------- 224 224 int ii; 225 intvec @w; //to store weights of all variables 226 @w[nvarBR] = 0; 227 @w = @w + 1; //initialize @w as 1..1 225 intvec @w=1:nvarBR; //to store weights of all variables 228 226 string str = "a"; //default for specifying elimination ordering 229 227 if (size(#) == 0) //default values … … 250 248 @w = #[1]; //take the given weights 251 249 str = #[2]; //string for specifying elimination ordering 252 253 250 } 254 251 if ( typeof(#[1]) == "string" and typeof(#[2]) == "intvec" ) … … 280 277 { 281 278 int local = (var(ii) < 1); 282 279 } 283 280 } 284 281 else … … 376 373 377 374 //define now the first a-block of the form a(w1,0..0) 378 intvec @v; 379 @v[nvarBR] = 0; 375 intvec @v=0:nvarBR; 380 376 @v = @v+w1; 381 377 B3[1] = list("a", @v); -
Singular/LIB/freegb.lib
r2c3a5d r3360fb 1 1 ////////////////////////////////////////////////////////////////////////////// 2 version="$Id: freegb.lib,v 1.2 1 2009-04-09 12:04:41 seelischExp $";2 version="$Id: freegb.lib,v 1.22 2009-04-14 12:00:14 Singular Exp $"; 3 3 category="Noncommutative"; 4 4 info=" … … 60 60 PURPOSE: sets attributes for a letterplace ring: 61 61 @* 'isLetterplaceRing' = true, 'uptodeg' = d, 'lV' = b, where 62 @* 'uptodeg' stands for the degree bound, 62 @* 'uptodeg' stands for the degree bound, 63 63 @* 'lV' for the number of variables in the block 0. 64 64 NOTE: Activate the resulting ring by using @code{setring} … … 69 69 ERROR("uptodeg and lV do not agree on the basering!"); 70 70 } 71 71 72 72 // Set letterplace-specific attributes for the output ring! 73 73 attrib(R, "uptodeg", uptodeg); 74 attrib(R, "lV", lV); 75 attrib(R, "isLetterplaceRing", 1); 76 return (R); 74 attrib(R, "lV", lV); 75 attrib(R, "isLetterplaceRing", 1); 76 return (R); 77 77 } 78 78 example … … 365 365 "USAGE: isVar(p); poly p 366 366 RETURN: int 367 PURPOSE: check, whether leading monomial of p is a power of a single variable 367 PURPOSE: check, whether leading monomial of p is a power of a single variable 368 368 @* from the basering. Returns the exponent or 0 if p is multivariate. 369 369 EXAMPLE: example isVar; shows examples … … 2148 2148 def R = makeLetterplaceRing(d); 2149 2149 setring R; 2150 int uptodeg = d; 2151 int lV = 2; 2150 int uptodeg = d; 2151 int lV = 2; 2152 2152 def R = setLetterplaceAttributes(r,uptodeg,2); // supply R with letterplace structure 2153 2153 setring R; -
Singular/LIB/jacobson.lib
r2c3a5d r3360fb 1 1 ////////////////////////////////////////////////////////////////////////////// 2 version="$Id: jacobson.lib,v 1.1 0 2009-03-05 20:36:11 levandovExp $";2 version="$Id: jacobson.lib,v 1.11 2009-04-14 12:00:14 Singular Exp $"; 3 3 category="System and Control Theory"; 4 4 info=" … … 13 13 @* with respect to the mult. closed set S = K[x] without 0. 14 14 @* Thus, we treat basering as principal ideal ring with d a polynomial 15 @* variable and x an invertible one. 15 @* variable and x an invertible one. 16 16 @* Note, that in computations no division by x will actually happen. 17 17 @* … … 22 22 REFERENCES: 23 23 @* (1) N. Jacobson, 'The theory of rings', AMS, 1943. 24 @* (2) Manuel Avelino Insua Hermo, 'Varias perspectives sobre las bases de Groebner : 25 @* Forma normal de Smith, Algorithme de Berlekamp y algebras de Leibniz'. 24 @* (2) Manuel Avelino Insua Hermo, 'Varias perspectives sobre las bases de Groebner : 25 @* Forma normal de Smith, Algorithme de Berlekamp y algebras de Leibniz'. 26 26 @* PhD thesis, Universidad de Santiago de Compostela, 2005. 27 27 … … 65 65 L = LL; 66 66 } 67 67 68 68 69 69 // divide units out … … 74 74 int i; int nsize = Min(intvec(nrM,ncM)); 75 75 poly p; number n; intvec lexp; 76 76 77 77 for(i=1; i<=nsize; i++) 78 78 { … … 123 123 RETURN: matrix or list of matrices, depending on arguments 124 124 ASSUME: Basering is a commutative polynomial ring in one variable 125 PURPOSE: compute the Smith Normal Form of M with (optionally) transformation matrices 125 PURPOSE: compute the Smith Normal Form of M with (optionally) transformation matrices 126 126 THEORY: Groebner bases are used for the Smith form like in (2). 127 127 NOTE: By default, just the Smith normal form of M is returned. … … 129 129 @* where U*M*V = D and the diagonal field entries of D are not normalized. 130 130 @* The normalization of the latter can be done with the 'divideUnits' procedure. 131 @* U and V above are square unimodular (invertible) matrices. 131 @* U and V above are square unimodular (invertible) matrices. 132 132 @* Note, that the procedure works for a rectangular matrix M. 133 133 @* -
Singular/LIB/nctools.lib
r2c3a5d r3360fb 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: nctools.lib,v 1.4 8 2009-04-14 08:25:19Singular Exp $";2 version="$Id: nctools.lib,v 1.49 2009-04-14 12:00:14 Singular Exp $"; 3 3 category="Noncommutative"; 4 4 info=" … … 901 901 { 902 902 print("Warning: Since the current ordering is not global there might be problems computing twostd(Q)!"); 903 "Q:"; 903 "Q:"; 904 904 @Q; 905 905 } -
Singular/LIB/polymake.lib
r2c3a5d r3360fb 1 version="$Id: polymake.lib,v 1.1 4 2009-04-07 16:18:06 seelischExp $";1 version="$Id: polymake.lib,v 1.15 2009-04-14 12:00:14 Singular Exp $"; 2 2 category="Tropical Geometry"; 3 3 info=" 4 LIBRARY: polymake.lib Computations with polytopes and fans, 4 LIBRARY: polymake.lib Computations with polytopes and fans, 5 5 interface to polymake and TOPCOM 6 6 AUTHOR: Thomas Markwig, email: keilen@mathematik.uni-kl.de … … 8 8 WARNING: 9 9 Most procedures will not work unless polymake or topcom is installed and 10 if so, they will only work with the operating system LINUX! 10 if so, they will only work with the operating system LINUX! 11 11 For more detailed information see the following note or consult the 12 12 help string of the procedures. 13 13 14 14 NOTE: 15 Even though this is a Singular library for computing polytopes and fans 16 such as the Newton polytope or the Groebner fan of a polynomial, most of 17 the hard computations are NOT done by Singular but by the program 15 Even though this is a Singular library for computing polytopes and fans 16 such as the Newton polytope or the Groebner fan of a polynomial, most of 17 the hard computations are NOT done by Singular but by the program 18 18 @* - polymake by Ewgenij Gawrilow, TU Berlin and Michael Joswig, TU Darmstadt 19 @* (see http://www.math.tu-berlin.de/polymake/), 19 @* (see http://www.math.tu-berlin.de/polymake/), 20 20 @* respectively (only in the procedure triangularions) by the program 21 21 @* - topcom by Joerg Rambau, Universitaet Bayreuth (see http://www.uni-bayreuth.de/ 22 22 departments/wirtschaftsmathematik/rambau/TOPCOM); 23 23 @* This library should rather be seen as an interface which allows to use a 24 (very limited) number of options which polymake respectively topcom offers 25 to compute with polytopes and fans and to make the results available in 24 (very limited) number of options which polymake respectively topcom offers 25 to compute with polytopes and fans and to make the results available in 26 26 Singular for further computations; 27 27 moreover, the user familiar with Singular does not have to learn the syntax 28 of polymake or topcom, if the options offered here are sufficient for his 28 of polymake or topcom, if the options offered here are sufficient for his 29 29 purposes. 30 @* Note, though, that the procedures concerned with planar polygons are 30 @* Note, though, that the procedures concerned with planar polygons are 31 31 independent of both, polymake and topcom. 32 32 … … 38 38 groebnerFan computes the Groebner fan of a polynomial 39 39 intmatToPolymake transforms an integer matrix into polymake format 40 polymakeToIntmat transforms polymake output into an integer matrix 40 polymakeToIntmat transforms polymake output into an integer matrix 41 41 42 42 PROCEDURES USING TOPCOM: … … 51 51 splitPolygon splits a marked polygon into vertices, facets, interior points 52 52 eta computes the eta-vector of a triangulation 53 findOrientedBoundary computes the boundary of a convex hull 53 findOrientedBoundary computes the boundary of a convex hull 54 54 cyclePoints computes lattice points connected to some lattice point 55 55 latticeArea computes the lattice area of a polygon … … 62 62 63 63 KEYWORDS: polytope; fan; secondary fan; secondary polytope; polymake; 64 Newton polytope; Groebner fan 64 Newton polytope; Groebner fan 65 65 "; 66 66 … … 95 95 proc polymakePolytope (intmat polytope,list #) 96 96 "USAGE: polymakePolytope(polytope[,#]); polytope list, # string 97 ASSUME: each row of polytope gives the coordinates of a lattice point of a 98 polytope with their affine coordinates as given by the output of 97 ASSUME: each row of polytope gives the coordinates of a lattice point of a 98 polytope with their affine coordinates as given by the output of 99 99 secondaryPolytope 100 PURPOSE: the procedure calls polymake to compute the vertices of the polytope 100 PURPOSE: the procedure calls polymake to compute the vertices of the polytope 101 101 as well as its dimension and information on its facets 102 102 RETURN: list L with four entries 103 103 @* L[1] : an integer matrix whose rows are the coordinates of vertices 104 of the polytope 104 of the polytope 105 105 @* L[2] : the dimension of the polytope 106 106 @* L[3] : a list whose i-th entry explains to which vertices the 107 ith vertex of the Newton polytope is connected 108 -- i.e. L[3][i] is an integer vector and an entry k in 109 there means that the vertex L[1][i] is connected to the 107 ith vertex of the Newton polytope is connected 108 -- i.e. L[3][i] is an integer vector and an entry k in 109 there means that the vertex L[1][i] is connected to the 110 110 vertex L[1][k] 111 @* L[4] : an integer matrix whose rows mulitplied by 112 (1,var(1),...,var(nvar)) give a linear system of equations 111 @* L[4] : an integer matrix whose rows mulitplied by 112 (1,var(1),...,var(nvar)) give a linear system of equations 113 113 describing the affine hull of the polytope, 114 114 i.e. the smallest affine space containing the polytope 115 NOTE: - for its computations the procedure calls the program polymake by 116 Ewgenij Gawrilow, TU Berlin and Michael Joswig, TU Darmstadt; 117 it therefore is necessary that this program is installed in order 115 NOTE: - for its computations the procedure calls the program polymake by 116 Ewgenij Gawrilow, TU Berlin and Michael Joswig, TU Darmstadt; 117 it therefore is necessary that this program is installed in order 118 118 to use this procedure; 119 119 see http://www.math.tu-berlin.de/polymake/ 120 @* - note that in the vertex edge graph we have changed the polymake 121 convention which starts indexing its vertices by zero while we start 120 @* - note that in the vertex edge graph we have changed the polymake 121 convention which starts indexing its vertices by zero while we start 122 122 with one ! 123 @* - the procedure creates the file /tmp/polytope.polymake which contains 124 the polytope in polymake format; if you wish to use this for further 125 computations with polymake, you have to use the procedure 123 @* - the procedure creates the file /tmp/polytope.polymake which contains 124 the polytope in polymake format; if you wish to use this for further 125 computations with polymake, you have to use the procedure 126 126 polymakeKeepTmpFiles in before 127 127 @* - moreover, the procedure creates the file /tmp/polytope.output which 128 128 it deletes again before ending 129 129 @* - it is possible to provide an optional second argument a string 130 which then will be used instead of 'polytope' in the name of the 130 which then will be used instead of 'polytope' in the name of the 131 131 polymake output file 132 132 EXAMPLE: example polymakePolytope; shows an example" … … 200 200 } 201 201 } 202 } 202 } 203 203 newveg=newveg[1,size(newveg)-1]; 204 204 execute("list nveg="+newveg+";"); … … 235 235 "EXAMPLE:"; 236 236 echo=2; 237 // the lattice points of the unit square in the plane 237 // the lattice points of the unit square in the plane 238 238 list points=intvec(0,0),intvec(0,1),intvec(1,0),intvec(1,1); 239 239 // the secondary polytope of this lattice point configuration is computed … … 261 261 of the Newton polytope of f 262 262 @* L[2] : the dimension of the Newton polytope of f 263 @* L[3] : a list whose ith entry explains to which vertices the 264 ith vertex of the Newton polytope is connected 265 -- i.e. L[3][i] is an integer vector and an entry k in 263 @* L[3] : a list whose ith entry explains to which vertices the 264 ith vertex of the Newton polytope is connected 265 -- i.e. L[3][i] is an integer vector and an entry k in 266 266 there means that the vertex L[1][i] is 267 267 connected to the vertex L[1][k] 268 @* L[4] : an integer matrix whose rows mulitplied by 269 (1,var(1),...,var(nvar)) give a linear system of equations 268 @* L[4] : an integer matrix whose rows mulitplied by 269 (1,var(1),...,var(nvar)) give a linear system of equations 270 270 describing the affine hull of the Newton polytope, i.e. the 271 271 smallest affine space containing the Newton polytope 272 NOTE: - if we replace the first column of L[4] by zeros, i.e. if we move 273 the affine hull to the origin, then we get the equations for the 274 orthogonal comploment of the linearity space of the normal fan dual 272 NOTE: - if we replace the first column of L[4] by zeros, i.e. if we move 273 the affine hull to the origin, then we get the equations for the 274 orthogonal comploment of the linearity space of the normal fan dual 275 275 to the Newton polytope, i.e. we get the EQUATIONS that 276 276 we need as input for polymake when computing the normal fan … … 278 278 TU Berlin and Michael Joswig, so it only works if polymake is installed; 279 279 see http://www.math.tu-berlin.de/polymake/ 280 @* - the procedure creates the file /tmp/newtonPolytope.polymake which 281 contains the polytope in polymake format and which can be used for 280 @* - the procedure creates the file /tmp/newtonPolytope.polymake which 281 contains the polytope in polymake format and which can be used for 282 282 further computations with polymake 283 @* - moreover, the procedure creates the file /tmp/newtonPolytope.output 283 @* - moreover, the procedure creates the file /tmp/newtonPolytope.output 284 284 and deletes it again before ending 285 @* - it is possible to give as an optional second argument a string 286 which then will be used instead of 'newtonPolytope' in the name of 285 @* - it is possible to give as an optional second argument a string 286 which then will be used instead of 'newtonPolytope' in the name of 287 287 the polymake output file 288 288 EXAMPLE: example newtonPolytope; shows an example" 289 289 { 290 290 int i,j; 291 // compute the list of exponent vectors of the polynomial, 291 // compute the list of exponent vectors of the polynomial, 292 292 // which are the lattice points 293 293 // whose convex hull is the Newton polytope of f … … 320 320 np[2]; 321 321 // np[3] contains information how the vertices are connected to each other, 322 // e.g. the first vertex (number 0) is connected to the second, third and 322 // e.g. the first vertex (number 0) is connected to the second, third and 323 323 // fourth vertex 324 324 np[3][1]; … … 330 330 np[1]; 331 331 // its dimension is 332 np[2]; 333 // the Newton polytope is contained in the affine space given 332 np[2]; 333 // the Newton polytope is contained in the affine space given 334 334 // by the equations 335 335 np[4]*M; … … 340 340 proc newtonPolytopeLP (poly f) 341 341 "USAGE: newtonPolytopeLP(f); f poly 342 RETURN: list, the exponent vectors of the monomials occuring in f, 342 RETURN: list, the exponent vectors of the monomials occuring in f, 343 343 i.e. the lattice points of the Newton polytope of f 344 344 EXAMPLE: example normalFan; shows an example" … … 368 368 proc normalFan (intmat vertices,intmat affinehull,list graph,int er,list #) 369 369 "USAGE: normalFan (vert,aff,graph,rays,[,#]); vert,aff intmat, graph list, rays int, # string 370 ASSUME: - vert is an integer matrix whose rows are the coordinate of 371 the vertices of a convex lattice polytope; 370 ASSUME: - vert is an integer matrix whose rows are the coordinate of 371 the vertices of a convex lattice polytope; 372 372 @* - aff describes the affine hull of this polytope, i.e. 373 the smallest affine space containing it, in the following sense: 374 denote by n the number of columns of vert, then multiply aff by 375 (1,x(1),...,x(n)) and set the resulting terms to zero in order to 373 the smallest affine space containing it, in the following sense: 374 denote by n the number of columns of vert, then multiply aff by 375 (1,x(1),...,x(n)) and set the resulting terms to zero in order to 376 376 get the equations for the affine hull; 377 @* - the ith entry of graph is an integer vector describing to which 378 vertices the ith vertex is connected, i.e. a k as entry means that 377 @* - the ith entry of graph is an integer vector describing to which 378 vertices the ith vertex is connected, i.e. a k as entry means that 379 379 the vertex vert[i] is connected to vert[k]; 380 @* - the integer rays is either one (if the extreme rays should be 380 @* - the integer rays is either one (if the extreme rays should be 381 381 computed) or zero (otherwise) 382 RETURN: list, the ith entry of L[1] contains information about the cone in the 383 normal fan dual to the ith vertex of the polytope 384 @* L[1][i][1] = integer matrix representing the inequalities which 382 RETURN: list, the ith entry of L[1] contains information about the cone in the 383 normal fan dual to the ith vertex of the polytope 384 @* L[1][i][1] = integer matrix representing the inequalities which 385 385 describe the cone dual to the ith vertex 386 @* L[1][i][2] = a list which contains the inequalities represented 387 by L[i][1] as a list of strings, where we use the 386 @* L[1][i][2] = a list which contains the inequalities represented 387 by L[i][1] as a list of strings, where we use the 388 388 variables x(1),...,x(n) 389 389 @* L[1][i][3] = only present if 'er' is set to 1; in that case it is 390 an interger matrix whose rows are the extreme rays 390 an interger matrix whose rows are the extreme rays 391 391 of the cone 392 @* L[2] = is an integer matrix whose rows span the linearity space 393 of the fan, i.e. the linear space which is contained in 392 @* L[2] = is an integer matrix whose rows span the linearity space 393 of the fan, i.e. the linear space which is contained in 394 394 each cone 395 395 NOTE: - the procedure calls for its computation polymake by Ewgenij Gawrilow, 396 TU Berlin and Michael Joswig, so it only works if polymake is 396 TU Berlin and Michael Joswig, so it only works if polymake is 397 397 installed; 398 398 see http://www.math.tu-berlin.de/polymake/ 399 @* - in the optional argument # it is possible to hand over other names 399 @* - in the optional argument # it is possible to hand over other names 400 400 for the variables to be used -- be careful, the format must be correct 401 401 which is not tested, e.g. if you want the variable names to be … … 405 405 list ineq; // stores the inequalities of the cones 406 406 int i,j,k; 407 // we work over the following ring 407 // we work over the following ring 408 408 execute("ring ineqring=0,x(1.."+string(ncols(vertices))+"),lp;"); 409 409 string greatersign=">"; … … 430 430 for (i=1;i<=nrows(vertices);i++) 431 431 { 432 // first we produce for each vertex in the polytope 432 // first we produce for each vertex in the polytope 433 433 // the inequalities describing the dual cone in the normal fan 434 list pp; // contain strings representing the inequalities 434 list pp; // contain strings representing the inequalities 435 435 // describing the normal cone 436 intmat ie[size(graph[i])][ncols(vertices)]; // contains the inequalities 436 intmat ie[size(graph[i])][ncols(vertices)]; // contains the inequalities 437 437 // as rows 438 // consider all the vertices to which the ith vertex in the 438 // consider all the vertices to which the ith vertex in the 439 439 // polytope is connected by an edge 440 440 for (j=1;j<=size(graph[i]);j++) 441 441 { 442 442 // produce the vector ie_j pointing from the jth vertex to the ith vertex; 443 // this will be the jth inequality for the cone in the normal fan dual to 443 // this will be the jth inequality for the cone in the normal fan dual to 444 444 // the ith vertex 445 445 ie[j,1..ncols(vertices)]=vertices[i,1..ncols(vertices)]-vertices[graph[i][j],1..ncols(vertices)]; … … 448 448 p=(VAR*EXP)[1,1]; 449 449 pl,pr=0,0; 450 // separate the terms with positive coefficients in p from 450 // separate the terms with positive coefficients in p from 451 451 // those with negative coefficients 452 452 for (k=1;k<=size(p);k++) … … 461 461 } 462 462 } 463 // build the string which represents the jth inequality 463 // build the string which represents the jth inequality 464 464 // for the cone dual to the ith vertex 465 // as polynomial inequality of type string, and store this 465 // as polynomial inequality of type string, and store this 466 466 // in the list pp as jth entry 467 467 pp[j]=string(pl)+" "+greatersign+" "+string(pr); … … 496 496 // create the file ineq.output 497 497 write(":w /tmp/ineq.output",""); 498 int dimension; // keeps the dimension of the intersection the 498 int dimension; // keeps the dimension of the intersection the 499 499 // bad cones with the u11tobeseencone 500 500 for (i=1;i<=size(ineq);i++) 501 501 { 502 i,". Cone of ",nrows(vertices); // indicate how many 502 i,". Cone of ",nrows(vertices); // indicate how many 503 503 // vertices have been dealt with 504 504 ungleichungen=intmatToPolymake(ineq[i][1],"rays"); … … 525 525 } 526 526 // get the linearity space 527 return(list(ineq,linearity)); 527 return(list(ineq,linearity)); 528 528 } 529 529 example … … 554 554 proc groebnerFan (poly f,list #) 555 555 "USAGE: groebnerFan(f[,#]); f poly, # string 556 RETURN: list, the ith entry of L[1] contains information about the ith cone 557 in the Groebner fan dual to the ith vertex in the Newton 556 RETURN: list, the ith entry of L[1] contains information about the ith cone 557 in the Groebner fan dual to the ith vertex in the Newton 558 558 polytope of the f 559 @* L[1][i][1] = integer matrix representing the inequalities 560 which describe the cone 561 @* L[1][i][2] = a list which contains the inequalities represented 559 @* L[1][i][1] = integer matrix representing the inequalities 560 which describe the cone 561 @* L[1][i][2] = a list which contains the inequalities represented 562 562 by L[1][i][1] as a list of strings 563 @* L[1][i][3] = an interger matrix whose rows are the extreme rays 563 @* L[1][i][3] = an interger matrix whose rows are the extreme rays 564 564 of the cone 565 @* L[2] = is an integer matrix whose rows span the linearity space 566 of the fan, i.e. the linear space which is contained 567 in each cone 568 @* L[3] = the Newton polytope of f in the format of the procedure 565 @* L[2] = is an integer matrix whose rows span the linearity space 566 of the fan, i.e. the linear space which is contained 567 in each cone 568 @* L[3] = the Newton polytope of f in the format of the procedure 569 569 newtonPolytope 570 @* L[4] = integer matrix where each row represents the exponet 570 @* L[4] = integer matrix where each row represents the exponet 571 571 vector of one monomial occuring in the input polynomial 572 572 NOTE: - if you have already computed the Newton polytope of f then you might want 573 to use the procedure normalFan instead in order to avoid doing costly 573 to use the procedure normalFan instead in order to avoid doing costly 574 574 computation twice 575 575 @* - the procedure calls for its computation polymake by Ewgenij Gawrilow, 576 576 TU Berlin and Michael Joswig, so it only works if polymake is installed; 577 577 see http://www.math.tu-berlin.de/polymake/ 578 @* - the procedure creates the file /tmp/newtonPolytope.polymake which 579 contains the Newton polytope of f in polymake format and which can 578 @* - the procedure creates the file /tmp/newtonPolytope.polymake which 579 contains the Newton polytope of f in polymake format and which can 580 580 be used for further computations with polymake 581 @* - it is possible to give as an optional second argument as string which 582 then will be used instead of 'newtonPolytope' in the name of the 581 @* - it is possible to give as an optional second argument as string which 582 then will be used instead of 'newtonPolytope' in the name of the 583 583 polymake output file 584 584 EXAMPLE: example groebnerFan; shows an example" 585 585 { 586 586 int i,j; 587 // compute the list of exponent vectors of the polynomial, which are 587 // compute the list of exponent vectors of the polynomial, which are 588 588 // the lattice points whose convex hull is the Newton polytope of f 589 589 intmat exponents[size(f)][nvars(basering)]; … … 649 649 proc intmatToPolymake (intmat M,string art) 650 650 "USAGE: intmatToPolymake(M,art); M intmat, art string 651 ASSUME: - M is an integer matrix which should be transformed into polymake 651 ASSUME: - M is an integer matrix which should be transformed into polymake 652 652 format; 653 653 @* - art is one of the following strings: 654 654 @* + 'rays' : indicating that a first column of 0's should be added 655 @* + 'points' : indicating that a first column of 1's should be added 656 RETURN: string, the matrix is transformed in a string and a first column has 655 @* + 'points' : indicating that a first column of 1's should be added 656 RETURN: string, the matrix is transformed in a string and a first column has 657 657 been added 658 658 EXAMPLE: example intmatToPolymake; shows an example" … … 662 662 string anf="0 "; 663 663 } 664 else 664 else 665 665 { 666 666 string anf="1 "; … … 697 697 proc polymakeToIntmat (string pm,string art) 698 698 "USAGE: polymakeToIntmat(pm,art); pm, art string 699 ASSUME: pm is the result of calling polymake with one 'argument' like 700 VERTICES, AFFINE_HULL, etc., so that the first row of the string is 701 the name of the corresponding 'argument' and the further rows contain 699 ASSUME: pm is the result of calling polymake with one 'argument' like 700 VERTICES, AFFINE_HULL, etc., so that the first row of the string is 701 the name of the corresponding 'argument' and the further rows contain 702 702 the result which consists of vectors either over the integers 703 703 or over the rationals … … 705 705 from the second row, where each row has been multiplied with the 706 706 lowest common multiple of the denominators of its entries as if 707 it is an integer matrix; moreover, if art=='affine', then 708 the first column is omitted since we only want affine 707 it is an integer matrix; moreover, if art=='affine', then 708 the first column is omitted since we only want affine 709 709 coordinates 710 710 EXAMPLE: example polymakeToIntmat; shows an example" … … 718 718 pm=stringdelete(pm,1); 719 719 } 720 pm=stringdelete(pm,1); 721 // find out how many entries each vector has, namely one more 720 pm=stringdelete(pm,1); 721 // find out how many entries each vector has, namely one more 722 722 // than 'spaces' in a row 723 723 int i=1; … … 734 734 // if we want to have affine coordinates 735 735 if (art=="affine") 736 { 736 { 737 737 s--; // then there is one column less 738 738 // and the entry of the first column (in the first row) has to be removed … … 743 743 pm=stringdelete(pm,1); 744 744 } 745 // we add two line breaks at the end in order to have this as 745 // we add two line breaks at the end in order to have this as 746 746 // a stopping criterion 747 747 pm=pm+zeilenumbruch+zeilenumbruch; … … 761 761 z++; 762 762 pm[i]=","; 763 // if we want to have affine coordinates, 763 // if we want to have affine coordinates, 764 764 // then we have to delete the first entry in each row 765 765 if (art=="affine") … … 775 775 if (pm[i]==" ") 776 776 { 777 pm[i]=","; 777 pm[i]=","; 778 778 } 779 779 } … … 784 784 pm=stringdelete(pm,size(pm)); 785 785 } 786 // since the matrix could be over the rationals, 786 // since the matrix could be over the rationals, 787 787 // we need a ring with rational coefficients 788 ring zwischering=0,x,lp; 788 ring zwischering=0,x,lp; 789 789 // create the matrix with the elements of pm as entries 790 790 execute("matrix mm["+string(z)+"]["+string(s)+"]="+pm+";"); … … 835 835 proc triangulations (list polygon) 836 836 "USAGE: triangulations(polygon); list polygon 837 ASSUME: polygon is a list of integer vectors of the same size representing 838 the affine coordinates of the lattice points 837 ASSUME: polygon is a list of integer vectors of the same size representing 838 the affine coordinates of the lattice points 839 839 PURPOSE: the procedure considers the marked polytope given as the convex hull of 840 840 the lattice points and with these lattice points as markings; it then 841 computes all possible triangulations of this marked polytope 841 computes all possible triangulations of this marked polytope 842 842 RETURN: list, each entry corresponds to one triangulation and the ith entry is 843 843 itself a list of integer vectors of size three, where each integer … … 846 846 NOTE:- the procedure calls for its computations the program points2triangs 847 847 from the program topcom by Joerg Rambau, Universitaet Bayreuth; it 848 therefore is necessary that this program is installed in order to use 848 therefore is necessary that this program is installed in order to use 849 849 this procedure; see 850 850 @* http://www.uni-bayreuth.de/departments/wirtschaftsmathematik/rambau/TOPCOM 851 @* - the procedure creates the files /tmp/triangulationsinput and 851 @* - the procedure creates the files /tmp/triangulationsinput and 852 852 /tmp/triangulationsoutput; 853 the former is used as input for points2triangs and the latter is its 854 output containing the triangulations of corresponding to points in the 855 format of points2triangs; if you wish to use this for further 856 computations with topcom, you have to use the procedure 853 the former is used as input for points2triangs and the latter is its 854 output containing the triangulations of corresponding to points in the 855 format of points2triangs; if you wish to use this for further 856 computations with topcom, you have to use the procedure 857 857 polymakeKeepTmpFiles in before 858 @* - note that an integer i in an integer vector representing a triangle 859 refers to the ith lattice point, i.e. polygon[i]; this convention is 860 different from TOPCOM's convention, where i would refer to the i-1st 858 @* - note that an integer i in an integer vector representing a triangle 859 refers to the ith lattice point, i.e. polygon[i]; this convention is 860 different from TOPCOM's convention, where i would refer to the i-1st 861 861 lattice point 862 862 EXAMPLE: example triangulations; shows an example" 863 863 { 864 864 int i,j; 865 // prepare the input for points2triangs by writing the input polygon in the 865 // prepare the input for points2triangs by writing the input polygon in the 866 866 // necessary format 867 867 string spi="["; … … 885 885 system("sh","cd /tmp; rm -f triangulationsinput; rm -f triangulationsoutput"); 886 886 } 887 // preprocessing of p2t if points2triangs is version >= 0.15 887 // preprocessing of p2t if points2triangs is version >= 0.15 888 888 // brings p2t to the format of version 0.14 889 889 string np2t; // takes the triangulations in Singular format … … 907 907 } 908 908 else 909 { 909 { 910 910 np2t=np2t+p2t[i]; 911 911 } … … 919 919 { 920 920 if (np2t[size(np2t)]!=";") 921 { 921 { 922 922 np2t=np2t+p2t[size(p2t)-1]+p2t[size(p2t)]; 923 923 } … … 941 941 np2t=np2t+"))"; 942 942 i++; 943 } 943 } 944 944 else 945 945 { … … 979 979 "EXAMPLE:"; 980 980 echo=2; 981 // the lattice points of the unit square in the plane 981 // the lattice points of the unit square in the plane 982 982 list polygon=intvec(0,0),intvec(0,1),intvec(1,0),intvec(1,1); 983 983 // the triangulations of this lattice point configuration are computed … … 1025 1025 int i,j,k,l; 1026 1026 intmat N[2][2]; // is used to compute areas of triangles 1027 intvec vertex; // stores a point in the secondary polytope as 1027 intvec vertex; // stores a point in the secondary polytope as 1028 1028 // intermediate result 1029 1029 int eintrag; 1030 1030 int halt; 1031 intmat secpoly[size(triangs)][size(polygon)]; // stores all lattice points 1031 intmat secpoly[size(triangs)][size(polygon)]; // stores all lattice points 1032 1032 // of the secondary polytope 1033 // consider each triangulation and compute the corresponding point 1033 // consider each triangulation and compute the corresponding point 1034 1034 // in the secondary polytope 1035 1035 for (i=1;i<=size(triangs);i++) 1036 1036 { 1037 // for each triangulation we have to compute the coordinates 1037 // for each triangulation we have to compute the coordinates 1038 1038 // corresponding to each marked point 1039 1039 for (j=1;j<=size(polygon);j++) 1040 1040 { 1041 1041 eintrag=0; 1042 // for each marked point we have to consider all triangles in the 1042 // for each marked point we have to consider all triangles in the 1043 1043 // triangulation which involve this particular point 1044 1044 for (k=1;k<=size(triangs[i]);k++) … … 1062 1062 secpoly[i,1..size(polygon)]=vertex; 1063 1063 } 1064 return(list(secpoly,triangs)); 1064 return(list(secpoly,triangs)); 1065 1065 } 1066 1066 example … … 1084 1084 proc secondaryFan (list polygon,list #) 1085 1085 "USAGE: secondaryFan(polygon[,#]); list polygon, list # 1086 ASSUME: - polygon is a list of integer vectors of the same size representing 1086 ASSUME: - polygon is a list of integer vectors of the same size representing 1087 1087 the affine coordinates of lattice points 1088 @* - if the triangulations of the corresponding polygon have already been 1089 computed with the procedure triangulations then these can be given 1090 as a second (optional) argument in order to avoid doing this 1088 @* - if the triangulations of the corresponding polygon have already been 1089 computed with the procedure triangulations then these can be given 1090 as a second (optional) argument in order to avoid doing this 1091 1091 computation again 1092 1092 PURPOSE: the procedure considers the marked polytope given as the convex hull of 1093 1093 the lattice points and with these lattice points as markings; it then 1094 computes the lattice points of the secondary polytope given by this 1094 computes the lattice points of the secondary polytope given by this 1095 1095 marked polytope which correspond to the triangulations computed by 1096 1096 the procedure triangulations 1097 RETURN: list, the ith entry of L[1] contains information about the ith cone in 1098 the secondary fan of the polygon, i.e. the cone dual to the 1097 RETURN: list, the ith entry of L[1] contains information about the ith cone in 1098 the secondary fan of the polygon, i.e. the cone dual to the 1099 1099 ith vertex of the secondary polytope 1100 @* L[1][i][1] = integer matrix representing the inequalities which 1100 @* L[1][i][1] = integer matrix representing the inequalities which 1101 1101 describe the cone dual to the ith vertex 1102 @* L[1][i][2] = a list which contains the inequalities represented 1102 @* L[1][i][2] = a list which contains the inequalities represented 1103 1103 by L[1][i][1] as a list of strings, where we use the 1104 1104 variables x(1),...,x(n) 1105 1105 @* L[1][i][3] = only present if 'er' is set to 1; in that case it is 1106 an interger matrix whose rows are the extreme rays 1106 an interger matrix whose rows are the extreme rays 1107 1107 of the cone 1108 @* L[2] = is an integer matrix whose rows span the linearity space 1109 of the fan, i.e. the linear space which is contained in 1108 @* L[2] = is an integer matrix whose rows span the linearity space 1109 of the fan, i.e. the linear space which is contained in 1110 1110 each cone 1111 @* L[3] = the secondary polytope in the format of the procedure 1111 @* L[3] = the secondary polytope in the format of the procedure 1112 1112 polymakePolytope 1113 @* L[4] = the list of triangulations corresponding to the vertices 1113 @* L[4] = the list of triangulations corresponding to the vertices 1114 1114 of the secondary polytope 1115 1115 NOTE:- the procedure calls for its computation polymake by Ewgenij Gawrilow, 1116 1116 TU Berlin and Michael Joswig, so it only works if polymake is installed; 1117 1117 see http://www.math.tu-berlin.de/polymake/ 1118 @* - in the optional argument # it is possible to hand over other names for 1118 @* - in the optional argument # it is possible to hand over other names for 1119 1119 the variables to be used -- be careful, the format must be correct 1120 1120 which is not tested, e.g. if you want the variable names to be 1121 1121 u00,u10,u01,u11 then you must hand over the string 'u11,u10,u01,u11' 1122 @* - if the triangluations are not handed over as optional argument the 1123 procedure calls for its computation of these triangulations the program 1124 points2triangs from the program topcom by Joerg Rambau, Universitaet 1125 Bayreuth; it therefore is necessary that this program is installed in 1122 @* - if the triangluations are not handed over as optional argument the 1123 procedure calls for its computation of these triangulations the program 1124 points2triangs from the program topcom by Joerg Rambau, Universitaet 1125 Bayreuth; it therefore is necessary that this program is installed in 1126 1126 order to use this procedure; see 1127 1127 @* http://www.uni-bayreuth.de/departments/wirtschaftsmathematik/rambau/TOPCOM … … 1172 1172 proc cycleLength (list boundary,intvec interior) 1173 1173 "USAGE: cycleLength(boundary,interior); list boundary, intvec interior 1174 ASSUME: boundary is a list of integer vectors describing a cycle in some 1175 convex lattice polygon around the lattice point interior ordered 1174 ASSUME: boundary is a list of integer vectors describing a cycle in some 1175 convex lattice polygon around the lattice point interior ordered 1176 1176 clock wise 1177 1177 RETURN: string, the cycle length of the corresponding cycle in the dual … … 1180 1180 { 1181 1181 int j; 1182 // create a ring whose variables are indexed by the points in 1182 // create a ring whose variables are indexed by the points in 1183 1183 // boundary resp. by interior 1184 1184 string rst="ring cyclering=0,(u"+string(interior[1])+string(interior[2]); … … 1218 1218 // interior is a lattice point in the interior of this lattice polygon 1219 1219 intvec interior=1,1; 1220 // compute the general cycle length of a cycle of the corresponding cycle 1220 // compute the general cycle length of a cycle of the corresponding cycle 1221 1221 // in the dual tropical curve, note that (0,1) and (2,1) do not contribute 1222 1222 cycleLength(boundary,interior); … … 1227 1227 proc splitPolygon (list markings) 1228 1228 "USAGE: splitPolygon (markings); markings list 1229 ASSUME: markings is a list of integer vectors representing lattice points in 1230 the plane which we consider as the marked points of the convex lattice 1229 ASSUME: markings is a list of integer vectors representing lattice points in 1230 the plane which we consider as the marked points of the convex lattice 1231 1231 polytope spanned by them 1232 PURPOSE: split the marked points in the vertices, the points on the facets 1232 PURPOSE: split the marked points in the vertices, the points on the facets 1233 1233 which are not vertices, and the interior points 1234 1234 RETURN: list, L consisting of three lists … … 1236 1236 @* L[1][i][1] = intvec, the coordinates of the ith vertex 1237 1237 @* L[1][i][2] = int, the position of L[1][i][1] in markings 1238 @* L[2][i] : represents the lattice points on the facet of the 1239 polygon with endpoints L[1][i] and L[1][i+1] 1238 @* L[2][i] : represents the lattice points on the facet of the 1239 polygon with endpoints L[1][i] and L[1][i+1] 1240 1240 (i considered modulo size(L[1])) 1241 @* L[2][i][j][1] = intvec, the coordinates of the jth 1241 @* L[2][i][j][1] = intvec, the coordinates of the jth 1242 1242 lattice point on that facet 1243 @* L[2][i][j][2] = int, the position of L[2][i][j][1] 1243 @* L[2][i][j][2] = int, the position of L[2][i][j][1] 1244 1244 in markings 1245 @* L[3] : represents the interior lattice points of the polygon 1245 @* L[3] : represents the interior lattice points of the polygon 1246 1246 @* L[3][i][1] = intvec, coordinates of ith interior point 1247 1247 @* L[3][i][2] = int, the position of L[3][i][1] in markings … … 1254 1254 vert[1]=pb[2]; 1255 1255 int i,j,k; // indices 1256 list boundary; // stores the points on the facets of the 1256 list boundary; // stores the points on the facets of the 1257 1257 // polygon which are not vertices 1258 // append to the boundary points as well as to the vertices 1258 // append to the boundary points as well as to the vertices 1259 1259 // the first vertex a second time 1260 1260 pb[1]=pb[1]+list(pb[1][1]); … … 1281 1281 // store the information on the boundary in vert[2] 1282 1282 vert[2]=boundary; 1283 // find the remaining points in the input which are not on 1283 // find the remaining points in the input which are not on 1284 1284 // the boundary by checking 1285 1285 // for each point in markings if it is contained in pb[1] … … 1298 1298 // store the interior points in vert[3] 1299 1299 vert[3]=interior; 1300 // add to each point in vert the index which it gets from 1300 // add to each point in vert the index which it gets from 1301 1301 // its position in the input markings; 1302 1302 // do so for ver[1] … … 1332 1332 } 1333 1333 vert[3][i]=list(vert[3][i],j); 1334 } 1334 } 1335 1335 return(vert); 1336 1336 } … … 1339 1339 "EXAMPLE:"; 1340 1340 echo=2; 1341 // the lattice polygon spanned by the points (0,0), (3,0) and (0,3) 1341 // the lattice polygon spanned by the points (0,0), (3,0) and (0,3) 1342 1342 // with all integer points as markings 1343 1343 list polygon=intvec(1,1),intvec(3,0),intvec(2,0),intvec(1,0), … … 1359 1359 proc eta (list triang,list polygon) 1360 1360 "USAGE: eta(triang,polygon); triang, polygon list 1361 ASSUME: polygon has the format of the output of splitPolygon, i.e. it is a 1362 list with three entries describing a convex lattice polygon in the 1361 ASSUME: polygon has the format of the output of splitPolygon, i.e. it is a 1362 list with three entries describing a convex lattice polygon in the 1363 1363 following way: 1364 @* polygon[1] : is a list of lists; for each i the entry polygon[1][i][1] 1365 is a lattice point which is a vertex of the lattice 1364 @* polygon[1] : is a list of lists; for each i the entry polygon[1][i][1] 1365 is a lattice point which is a vertex of the lattice 1366 1366 polygon, and polygon[1][i][2] is an integer assigned to 1367 1367 this lattice point as identifying index 1368 @* polygon[2] : is a list of lists; for each vertex of the polygon, 1369 i.e. for each entry in polygon[1], it contains a list 1370 polygon[2][i], which contains the lattice points on the 1371 facet with endpoints polygon[1][i] and polygon[1][i+1] 1368 @* polygon[2] : is a list of lists; for each vertex of the polygon, 1369 i.e. for each entry in polygon[1], it contains a list 1370 polygon[2][i], which contains the lattice points on the 1371 facet with endpoints polygon[1][i] and polygon[1][i+1] 1372 1372 - i considered mod size(polygon[1]); 1373 each such lattice point contributes an entry 1373 each such lattice point contributes an entry 1374 1374 polygon[2][i][j][1] which is an integer 1375 vector giving the coordinate of the lattice point and an 1375 vector giving the coordinate of the lattice point and an 1376 1376 entry polygon[2][i][j][2] which is the identifying index 1377 @* polygon[3] : is a list of lists, where each entry corresponds to a 1378 lattice point in the interior of the polygon, with 1377 @* polygon[3] : is a list of lists, where each entry corresponds to a 1378 lattice point in the interior of the polygon, with 1379 1379 polygon[3][j][1] being the coordinates of the point 1380 1380 and polygon[3][j][2] being the identifying index; 1381 @* triang is a list of integer vectors all of size three describing a 1382 triangulation of the polygon described by polygon; if an entry of 1381 @* triang is a list of integer vectors all of size three describing a 1382 triangulation of the polygon described by polygon; if an entry of 1383 1383 triang is the vector (i,j,k) then the triangle is built from the vertices 1384 1384 with indices i, j and k 1385 RETURN: intvec, the integer vector eta describing that vertex of the Newton 1386 polytope discriminant of the polygone whose dual cone in the 1387 Groebner fan contains the cone of the secondary fan of the 1385 RETURN: intvec, the integer vector eta describing that vertex of the Newton 1386 polytope discriminant of the polygone whose dual cone in the 1387 Groebner fan contains the cone of the secondary fan of the 1388 1388 polygon corresponding to the given triangulation 1389 NOTE: for a better description of eta see Gelfand, Kapranov, 1389 NOTE: for a better description of eta see Gelfand, Kapranov, 1390 1390 Zelevinski: Discriminants, Resultants and multidimensional Determinants. 1391 1391 Chapter 10. … … 1393 1393 { 1394 1394 int i,j,k,l,m,n; // index variables 1395 list ordpolygon; // stores the lattice points in the order 1395 list ordpolygon; // stores the lattice points in the order 1396 1396 // used in the triangulation 1397 1397 list triangarea; // stores the areas of the triangulations … … 1419 1419 for (i=1;i<=size(triang);i++) 1420 1420 { 1421 // Note that the ith lattice point in orderedpolygon has the 1421 // Note that the ith lattice point in orderedpolygon has the 1422 1422 // number i-1 in the triangulation! 1423 1423 N=ordpolygon[triang[i][1]]-ordpolygon[triang[i][2]],ordpolygon[triang[i][1]]-ordpolygon[triang[i][3]]; … … 1425 1425 } 1426 1426 intvec ETA; // stores the eta_ij 1427 int etaij; // stores the part of eta_ij during computations 1427 int etaij; // stores the part of eta_ij during computations 1428 1428 // which comes from triangle areas 1429 int seitenlaenge; // stores the part of eta_ij during computations 1429 int seitenlaenge; // stores the part of eta_ij during computations 1430 1430 // which comes from boundary facets 1431 1431 list seiten; // stores the lattice points on facets of the polygon 1432 1432 intvec v; // used to compute a facet length 1433 // 3) store first in seiten[i] all lattice points on the facet 1433 // 3) store first in seiten[i] all lattice points on the facet 1434 1434 // connecting the ith vertex, 1435 // i.e. polygon[1][i], with the i+1st vertex, i.e. polygon[1][i+1], 1435 // i.e. polygon[1][i], with the i+1st vertex, i.e. polygon[1][i+1], 1436 1436 // where we replace i+1 1437 1437 // 1 if i=size(polygon[1]); 1438 // then append the last entry of seiten once more at the very 1438 // then append the last entry of seiten once more at the very 1439 1439 // beginning of seiten, so 1440 1440 // that the index is shifted by one … … 1462 1462 if ((polygon[1][j][2]==triang[k][1]) or (polygon[1][j][2]==triang[k][2]) or (polygon[1][j][2]==triang[k][3])) 1463 1463 { 1464 // ... if so, add the area of the triangle to etaij 1464 // ... if so, add the area of the triangle to etaij 1465 1465 etaij=etaij+triangarea[k]; 1466 // then check if that triangle has a facet which is contained 1467 // in one of the 1466 // then check if that triangle has a facet which is contained 1467 // in one of the 1468 1468 // two facets of the polygon which are adjecent to the given vertex ... 1469 1469 // these two facets are seiten[j] and seiten[j+1] … … 1479 1479 if ((seiten[n][l][2]==triang[k][m]) and (seiten[n][l][2]!=polygon[1][j][2])) 1480 1480 { 1481 // if so, then compute the vector pointing from this 1481 // if so, then compute the vector pointing from this 1482 1482 // lattice point to the vertex 1483 1483 v=polygon[1][j][1]-seiten[n][l][1]; 1484 // and the lattice length of this vector has to be 1484 // and the lattice length of this vector has to be 1485 1485 // subtracted from etaij 1486 1486 etaij=etaij-abs(gcd(v[1],v[2])); … … 1494 1494 ETA[polygon[1][j][2]]=etaij; 1495 1495 } 1496 // 5) compute the eta_ij for all lattice points on the facets 1496 // 5) compute the eta_ij for all lattice points on the facets 1497 1497 // of the polygon which are not vertices, these are the 1498 1498 // lattice points in polygon[2][1] to polygon[2][size(polygon[1])] … … 1500 1500 { 1501 1501 for (j=1;j<=size(polygon[2][i]);j++) 1502 { 1502 { 1503 1503 // initialise etaij 1504 1504 etaij=0; … … 1511 1511 if ((polygon[2][i][j][2]==triang[k][1]) or (polygon[2][i][j][2]==triang[k][2]) or (polygon[2][i][j][2]==triang[k][3])) 1512 1512 { 1513 // ... if so, add the area of the triangle to etaij 1513 // ... if so, add the area of the triangle to etaij 1514 1514 etaij=etaij+triangarea[k]; 1515 // then check if that triangle has a facet which is contained in the 1515 // then check if that triangle has a facet which is contained in the 1516 1516 // facet of the polygon which contains the lattice point in question, 1517 1517 // this is the facet seiten[i+1]; … … 1521 1521 // ... and for each lattice point in the triangle ... 1522 1522 for (m=1;m<=size(triang[k]);m++) 1523 { 1523 { 1524 1524 // ... if they coincide and are not the vertex itself ... 1525 1525 if ((seiten[i+1][l][2]==triang[k][m]) and (seiten[i+1][l][2]!=polygon[2][i][j][2])) 1526 1526 { 1527 // if so, then compute the vector pointing from 1527 // if so, then compute the vector pointing from 1528 1528 // this lattice point to the vertex 1529 1529 v=polygon[2][i][j][1]-seiten[i+1][l][1]; 1530 // and the lattice length of this vector contributes 1530 // and the lattice length of this vector contributes 1531 1531 // to seitenlaenge 1532 1532 seitenlaenge=seitenlaenge+abs(gcd(v[1],v[2])); … … 1536 1536 } 1537 1537 } 1538 // if the lattice point was a vertex of any triangle 1538 // if the lattice point was a vertex of any triangle 1539 1539 // in the triangulation ... 1540 1540 if (etaij!=0) … … 1561 1561 if ((polygon[3][j][2]==triang[k][1]) or (polygon[3][j][2]==triang[k][2]) or (polygon[3][j][2]==triang[k][3])) 1562 1562 { 1563 // ... if so, add the area of the triangle to etaij 1563 // ... if so, add the area of the triangle to etaij 1564 1564 etaij=etaij+triangarea[k]; 1565 1565 } … … 1574 1574 "EXAMPLE:"; 1575 1575 echo=2; 1576 // the lattice polygon spanned by the points (0,0), (3,0) and (0,3) 1576 // the lattice polygon spanned by the points (0,0), (3,0) and (0,3) 1577 1577 // with all integer points as markings 1578 1578 list polygon=intvec(1,1),intvec(3,0),intvec(2,0),intvec(1,0), … … 1581 1581 // split the polygon in its vertices, its facets and its interior points 1582 1582 list sp=splitPolygon(polygon); 1583 // define a triangulation by connecting the only interior point 1583 // define a triangulation by connecting the only interior point 1584 1584 // with the vertices 1585 1585 list triang=intvec(1,2,5),intvec(1,5,10),intvec(1,5,10); … … 1587 1587 eta(triang,sp); 1588 1588 } 1589 1589 1590 1590 ///////////////////////////////////////////////////////////////////////////// 1591 1591 … … 1614 1614 } 1615 1615 // check is the polygon is only a line segment given by more than two points; 1616 // for this first compute sum of the absolute values of the determinants 1616 // for this first compute sum of the absolute values of the determinants 1617 1617 // of the matrices whose 1618 // rows are the vectors pointing from the first to the second point 1618 // rows are the vectors pointing from the first to the second point 1619 1619 // and from the 1620 // the first point to the ith point for i=3,...,size(polygon); 1620 // the first point to the ith point for i=3,...,size(polygon); 1621 1621 // if this sum is zero 1622 1622 // then the polygon is a line segment and we have to find its end points … … 1631 1631 intmat laenge[size(polygon)][size(polygon)]; 1632 1632 intvec mp; 1633 // for this collect first all vectors pointing from one lattice 1633 // for this collect first all vectors pointing from one lattice 1634 1634 // point to the next, 1635 1635 // compute their pairwise angles and their lengths 1636 1636 for (i=1;i<=size(polygon)-1;i++) 1637 { 1637 { 1638 1638 for (j=i+1;j<=size(polygon);j++) 1639 1639 { … … 1659 1659 polygon=sortlistbyintvec(polygon,abstand); 1660 1660 return(list(polygon,endpoints)); 1661 } 1661 } 1662 1662 /////////////////////////////////////////////////////////////// 1663 1663 list orderedvertices; // stores the vertices in an ordered way 1664 list minimisedorderedvertices; // stores the vertices in an ordered way; 1664 list minimisedorderedvertices; // stores the vertices in an ordered way; 1665 1665 // redundant ones removed 1666 list comparevertices; // stores vertices which should be compared to 1666 list comparevertices; // stores vertices which should be compared to 1667 1667 // the testvertex 1668 1668 orderedvertices[1]=polygon[1]; // set the starting vertex 1669 1669 minimisedorderedvertices[1]=polygon[1]; // set the starting vertex 1670 1670 intvec testvertex=polygon[1]; //vertex to which the others have to be compared 1671 intvec startvertex=polygon[1]; // keep the starting vertex to test, 1671 intvec startvertex=polygon[1]; // keep the starting vertex to test, 1672 1672 // when the end is reached 1673 1673 int endtest; // is set to one, when the end is reached 1674 int startvertexfound;// is 1, once for some testvertex a candidate 1675 // for the next vertex has been found 1674 int startvertexfound;// is 1, once for some testvertex a candidate 1675 // for the next vertex has been found 1676 1676 polygon=delete(polygon,1); // delete the testvertex 1677 1677 intvec v,w; 1678 1678 int l=1; // counts the vertices 1679 // the basic idea is that a vertex can be 1679 // the basic idea is that a vertex can be 1680 1680 // the next one on the boundary if all other vertices 1681 // lie to the right of the vector v pointing 1681 // lie to the right of the vector v pointing 1682 1682 // from the testvertex to this one; this can be tested 1683 // by checking if the determinant of the 2x2-matrix 1683 // by checking if the determinant of the 2x2-matrix 1684 1684 // with first column v and second column the vector w, 1685 // pointing from the testvertex to the new vertex, 1685 // pointing from the testvertex to the new vertex, 1686 1686 // is non-positive; if this is the case for all 1687 // new vertices, then the one in consideration is 1687 // new vertices, then the one in consideration is 1688 1688 // a possible choice for the next vertex on the boundary 1689 // and it is stored in naechste; we can then order 1689 // and it is stored in naechste; we can then order 1690 1690 // the candidates according to their distance from 1691 1691 // the testvertex; then they occur on the boundary in that order! … … 1699 1699 v=polygon[i]-testvertex; // points from the testvertex to the ith vertex 1700 1700 comparevertices=delete(polygon,i); // we needn't compare v to itself 1701 // we should compare v to the startvertex-testvertex; 1701 // we should compare v to the startvertex-testvertex; 1702 1702 // in the first calling of the loop 1703 // this is irrelevant since the difference will be zero; 1703 // this is irrelevant since the difference will be zero; 1704 1704 // however, later on it will 1705 // be vital, since we delete the vertices 1705 // be vital, since we delete the vertices 1706 1706 // which we have already tested from the list 1707 // of all vertices, and when all vertices 1707 // of all vertices, and when all vertices 1708 1708 // on the boundary have been found we would 1709 // therefore find a vertex in the interior 1709 // therefore find a vertex in the interior 1710 1710 // as candidate; but always testing against 1711 1711 // the starting vertex, this can not happen 1712 comparevertices[size(comparevertices)+1]=startvertex; 1712 comparevertices[size(comparevertices)+1]=startvertex; 1713 1713 for (j=1;(j<=size(comparevertices)) and (d<=0);j++) 1714 1714 { … … 1718 1718 d=det(D); 1719 1719 } 1720 if (d<=0) // if all determinants are non-positive, 1720 if (d<=0) // if all determinants are non-positive, 1721 1721 { // then the ith vertex is a candidate 1722 1722 naechste[k]=list(polygon[i],i,scalarproduct(v,v));// we store the vertex, … … 1726 1726 } 1727 1727 if (size(naechste)>0) // then a candidate for the next vertex has been found 1728 { 1728 { 1729 1729 startvertexfound=1; // at least once a candidate has been found 1730 naechste=sortlist(naechste,3); // we order the candidates according 1730 naechste=sortlist(naechste,3); // we order the candidates according 1731 1731 // to their distance from testvertex; 1732 for (j=1;j<=size(naechste);j++) // then we store them in this 1732 for (j=1;j<=size(naechste);j++) // then we store them in this 1733 1733 { // order in orderedvertices 1734 1734 l++; 1735 1735 orderedvertices[l]=naechste[j][1]; 1736 1736 } 1737 testvertex=naechste[size(naechste)][1]; // we store the last one as 1737 testvertex=naechste[size(naechste)][1]; // we store the last one as 1738 1738 // next testvertex; 1739 1739 // store the next corner of NSD 1740 minimisedorderedvertices[size(minimisedorderedvertices)+1]=testvertex; 1741 naechste=sortlist(naechste,2); // then we reorder the vertices 1740 minimisedorderedvertices[size(minimisedorderedvertices)+1]=testvertex; 1741 naechste=sortlist(naechste,2); // then we reorder the vertices 1742 1742 // according to their position 1743 1743 for (j=size(naechste);j>=1;j--) // and we delete them from the vertices … … 1746 1746 } 1747 1747 } 1748 else // that means either that the vertex was inside the polygon, 1749 { // or that we have reached the last vertex on the boundary 1748 else // that means either that the vertex was inside the polygon, 1749 { // or that we have reached the last vertex on the boundary 1750 1750 // of the polytope 1751 if (startvertexfound==0) // the vertex was in the interior; 1751 if (startvertexfound==0) // the vertex was in the interior; 1752 1752 { // we delete it and start all over again 1753 orderedvertices[1]=polygon[1]; 1754 minimisedorderedvertices[1]=polygon[1]; 1753 orderedvertices[1]=polygon[1]; 1754 minimisedorderedvertices[1]=polygon[1]; 1755 1755 testvertex=polygon[1]; 1756 1756 startvertex=polygon[1]; 1757 1757 polygon=delete(polygon,1); 1758 1758 } 1759 else // we have reached the last vertex on the boundary of 1759 else // we have reached the last vertex on the boundary of 1760 1760 { // the polytope and can stop 1761 1761 endtest=1; … … 1764 1764 kill naechste; 1765 1765 } 1766 // test if the first vertex in minimisedorderedvertices 1766 // test if the first vertex in minimisedorderedvertices 1767 1767 // is on the same line with the second and 1768 // the last, i.e. if we started our search in the 1768 // the last, i.e. if we started our search in the 1769 1769 // middle of a face; if so, delete it 1770 1770 v=minimisedorderedvertices[2]-minimisedorderedvertices[1]; … … 1775 1775 minimisedorderedvertices=delete(minimisedorderedvertices,1); 1776 1776 } 1777 // test if the first vertex in minimisedorderedvertices 1777 // test if the first vertex in minimisedorderedvertices 1778 1778 // is on the same line with the two 1779 // last ones, i.e. if we started our search at the end of a face; 1779 // last ones, i.e. if we started our search at the end of a face; 1780 1780 // if so, delete it 1781 1781 v=minimisedorderedvertices[size(minimisedorderedvertices)-1]-minimisedorderedvertices[1]; … … 1809 1809 proc cyclePoints (list triang,list points,int pt) 1810 1810 "USAGE: cyclePoints(triang,points,pt) triang,points list, pt int 1811 ASSUME: - points is a list of integer vectors describing the lattice 1811 ASSUME: - points is a list of integer vectors describing the lattice 1812 1812 points of a marked polygon; 1813 @* - triang is a list of integer vectors describing a triangulation 1814 of the marked polygon in the sense that an integer vector of 1815 the form (i,j,k) describes the triangle formed by polygon[i], 1813 @* - triang is a list of integer vectors describing a triangulation 1814 of the marked polygon in the sense that an integer vector of 1815 the form (i,j,k) describes the triangle formed by polygon[i], 1816 1816 polygon[j] and polygon[k]; 1817 @* - pt is an integer between 1 and size(points), singling out a 1817 @* - pt is an integer between 1 and size(points), singling out a 1818 1818 lattice point among the marked points 1819 PURPOSE: consider the convex lattice polygon, say P, spanned by all lattice 1820 points in points which in the triangulation triang are connected 1821 to the point points[pt]; the procedure computes all marked points 1819 PURPOSE: consider the convex lattice polygon, say P, spanned by all lattice 1820 points in points which in the triangulation triang are connected 1821 to the point points[pt]; the procedure computes all marked points 1822 1822 in points which lie on the boundary of that polygon, ordered 1823 1823 clockwise 1824 RETURN: list, of integer vectors which are the coordinates of the lattice 1825 points on the boundary of the above mentioned polygon P, if 1826 this polygon is not the empty set (that would be the case if 1827 points[pt] is not a vertex of any triangle in the 1824 RETURN: list, of integer vectors which are the coordinates of the lattice 1825 points on the boundary of the above mentioned polygon P, if 1826 this polygon is not the empty set (that would be the case if 1827 points[pt] is not a vertex of any triangle in the 1828 1828 triangulation); otherwise return the empty list 1829 1829 EXAMPLE: example cyclePoints; shows an example" 1830 1830 { 1831 1831 int i,j; // indices 1832 list v; // saves the indices of lattice points connected to the 1832 list v; // saves the indices of lattice points connected to the 1833 1833 // interior point in the triangulation 1834 1834 // save all points in triangulations containing pt in v … … 1866 1866 pts[i]=points[v[i]]; 1867 1867 } 1868 // consider the convex polytope spanned by the points in pts, 1868 // consider the convex polytope spanned by the points in pts, 1869 1869 // find the points on the 1870 1870 // boundary and order them clockwise … … 1875 1875 "EXAMPLE:"; 1876 1876 echo=2; 1877 // the lattice polygon spanned by the points (0,0), (3,0) and (0,3) 1877 // the lattice polygon spanned by the points (0,0), (3,0) and (0,3) 1878 1878 // with all integer points as markings 1879 1879 list points=intvec(1,1),intvec(3,0),intvec(2,0),intvec(1,0), 1880 1880 intvec(0,0),intvec(2,1),intvec(0,1),intvec(1,2), 1881 1881 intvec(0,2),intvec(0,3); 1882 // define a triangulation 1882 // define a triangulation 1883 1883 list triang=intvec(1,2,5),intvec(1,5,7),intvec(1,7,9),intvec(8,9,10), 1884 1884 intvec(1,8,9),intvec(1,2,8); … … 1892 1892 "USAGE: latticeArea(polygon); polygon list 1893 1893 ASSUME: polygon is a list of integer vectors in the plane 1894 RETURN: int, the lattice area of the convex hull of the lattice points in 1894 RETURN: int, the lattice area of the convex hull of the lattice points in 1895 1895 polygon, i.e. twice the Euclidean area 1896 1896 EXAMPLE: example polygonlatticeArea; shows an example" … … 1921 1921 proc picksFormula (list polygon) 1922 1922 "USAGE: picksFormula(polygon); polygon list 1923 ASSUME: polygon is a list of integer vectors in the plane and consider their 1924 convex hull C 1925 RETURN: list, L of three integersthe 1923 ASSUME: polygon is a list of integer vectors in the plane and consider their 1924 convex hull C 1925 RETURN: list, L of three integersthe 1926 1926 @* L[1] : the lattice area of C, i.e. twice the Euclidean area 1927 1927 @* L[2] : the number of lattice points on the boundary of C … … 1948 1948 bdpts=bdpts+abs(gcd(edge[1],edge[2])); 1949 1949 } 1950 // Pick's formula says that the lattice area A, the number g of interior 1950 // Pick's formula says that the lattice area A, the number g of interior 1951 1951 // points and 1952 1952 // the number b of boundary points are connected by the formula: A=b+2g-2 … … 1976 1976 "USAGE: ellipticNF(polygon); polygon list 1977 1977 ASSUME: polygon is a list of integer vectors in the plane such that their 1978 convex hull C has precisely one interior lattice point, i.e. C is the 1978 convex hull C has precisely one interior lattice point, i.e. C is the 1979 1979 Newton polygon of an elliptic curve 1980 PURPOSE: compute the normal form of the polygon with respect to the unimodular 1980 PURPOSE: compute the normal form of the polygon with respect to the unimodular 1981 1981 affine transformations T=A*x+v; there are sixteen different normal forms 1982 (see e.g. Bjorn Poonen, Fernando Rodriguez-Villegas: Lattice Polygons 1983 and the number 12. Amer. Math. Monthly 107 (2000), no. 3, 1982 (see e.g. Bjorn Poonen, Fernando Rodriguez-Villegas: Lattice Polygons 1983 and the number 12. Amer. Math. Monthly 107 (2000), no. 3, 1984 1984 238--250.) 1985 1985 RETURN: list, L such that 1986 @* L[1] : list whose entries are the vertices of the normal form of 1986 @* L[1] : list whose entries are the vertices of the normal form of 1987 1987 the polygon 1988 1988 @* L[2] : the matrix A of the unimodular transformation 1989 1989 @* L[3] : the translation vector v of the unimodular transformation 1990 @* L[4] : list such that the ith entry is the image of polygon[i] 1990 @* L[4] : list such that the ith entry is the image of polygon[i] 1991 1991 under the unimodular transformation T 1992 1992 EXAMPLE: example ellipticNF; shows an example" … … 2020 2020 intvec trans; // stores the vector by which we have to translate the polygon 2021 2021 intmat A[2][2]; // stores the matrix by which we have to transform the polygon 2022 matrix M[3][3]; // stores the projective coordinates of the points 2022 matrix M[3][3]; // stores the projective coordinates of the points 2023 2023 // which are to be transformed 2024 matrix N[3][3]; // stores the projective coordinates of the points to 2024 matrix N[3][3]; // stores the projective coordinates of the points to 2025 2025 // which M is to be transformed 2026 intmat T[3][3]; // stores the unimodular affine transformation in 2026 intmat T[3][3]; // stores the unimodular affine transformation in 2027 2027 // projective form 2028 2028 // add the second point of pg once again at the end 2029 2029 pg=insert(pg,pg[2],size(pg)); 2030 // if there is only one edge which has the maximal number of lattice points, 2030 // if there is only one edge which has the maximal number of lattice points, 2031 2031 // then M should be: 2032 2032 M=pg[max],1,pg[max+1],1,pg[max+2],1; … … 2118 2118 M=pg[max],1,pg[max+1],1,pg[max+2],1; 2119 2119 // the orientation of the polygon matters 2120 A=pg[max-1]-pg[max],pg[max+1]-pg[max]; 2120 A=pg[max-1]-pg[max],pg[max+1]-pg[max]; 2121 2121 if (det(A)==4) 2122 2122 { … … 2167 2167 { 2168 2168 max++; 2169 } 2169 } 2170 2170 M=pg[max],1,pg[max+1],1,pg[max+2],1; 2171 2171 N=0,1,1,1,2,1,2,1,1; … … 2230 2230 // the vertices of the normal form are 2231 2231 nf[1]; 2232 // it has been transformed by the unimodular affine transformation A*x+v 2232 // it has been transformed by the unimodular affine transformation A*x+v 2233 2233 // with matrix A 2234 2234 nf[2]; … … 2247 2247 "USAGE: ellipticNFDB(n[,#]); n int, # list 2248 2248 ASSUME: n is an integer between 1 and 16 2249 PURPOSE: this is a database storing the 16 normal forms of planar polygons with 2249 PURPOSE: this is a database storing the 16 normal forms of planar polygons with 2250 2250 precisely one interior point up to unimodular affine transformations 2251 @* (see e.g. Bjorn Poonen, Fernando Rodriguez-Villegas: Lattice Polygons 2251 @* (see e.g. Bjorn Poonen, Fernando Rodriguez-Villegas: Lattice Polygons 2252 2252 and the number 12. Amer. Math. Monthly 107 (2000), no. 3, 2253 2253 238--250.) 2254 2254 RETURN: list, L such that 2255 @* L[1] : list whose entries are the vertices of the nth normal form 2256 @* L[2] : list whose entries are all the lattice points of the 2257 nth normal form 2258 @* L[3] : only present if the optional parameter # is present, and 2259 then it is a polynomial in the variables (x,y) whose 2255 @* L[1] : list whose entries are the vertices of the nth normal form 2256 @* L[2] : list whose entries are all the lattice points of the 2257 nth normal form 2258 @* L[3] : only present if the optional parameter # is present, and 2259 then it is a polynomial in the variables (x,y) whose 2260 2260 Newton polygon is the nth normal form 2261 NOTE: the optional parameter is only allowed if the basering has the 2261 NOTE: the optional parameter is only allowed if the basering has the 2262 2262 variables x and y 2263 2263 EXAMPLE: example ellipticNFDB; shows an example" … … 2310 2310 proc polymakeKeepTmpFiles (int i) 2311 2311 "USAGE: polymakeKeepTmpFiles(int i); i int 2312 PURPOSE: some procedures create files in the directory /tmp which are used for 2312 PURPOSE: some procedures create files in the directory /tmp which are used for 2313 2313 computations with polymake respectively topcom; these will be removed 2314 when the corresponding procedure is left; however, it might be 2314 when the corresponding procedure is left; however, it might be 2315 2315 desireable to keep them for further computations with either polymake or 2316 2316 topcom; this can be achieved by this procedure; call the procedure as: … … 2355 2355 static proc scalarproduct (intvec w,intvec v) 2356 2356 "USAGE: scalarproduct(w,v); w,v intvec 2357 ASSUME: w and v are integer vectors of the same length 2357 ASSUME: w and v are integer vectors of the same length 2358 2358 RETURN: int, the scalarproduct of v and w 2359 2359 NOTE: the procedure is called by findOrientedBoundary" … … 2402 2402 { 2403 2403 int m=nrows(M); 2404 2404 2405 2405 } 2406 2406 else … … 2460 2460 { 2461 2461 return(""); 2462 2462 2463 2463 } 2464 2464 if (i==1) … … 2570 2570 k++; 2571 2571 } 2572 else 2572 else 2573 2573 { 2574 2574 stop=1; … … 2613 2613 k++; 2614 2614 } 2615 else 2615 else 2616 2616 { 2617 2617 stop=1; … … 2662 2662 static proc polygonToCoordinates (list points) 2663 2663 "USAGE: polygonToCoordinates(points); points list 2664 ASSUME: points is a list of integer vectors each of size two describing the 2665 marked points of a convex lattice polygon like the output of 2664 ASSUME: points is a list of integer vectors each of size two describing the 2665 marked points of a convex lattice polygon like the output of 2666 2666 polygonDB 2667 RETURN: list, the first entry is a string representing the coordinates 2667 RETURN: list, the first entry is a string representing the coordinates 2668 2668 corresponding to the latticpoints seperated by commata 2669 the second entry is a list where the ith entry is a string 2670 representing the coordinate of corresponding to the ith 2671 lattice point the third entry is the latex format of the 2669 the second entry is a list where the ith entry is a string 2670 representing the coordinate of corresponding to the ith 2671 lattice point the third entry is the latex format of the 2672 2672 first entry 2673 2673 NOTE: the procedure is called by fan" -
Singular/LIB/presolve.lib
r2c3a5d r3360fb 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: presolve.lib,v 1. 29 2009-04-06 09:17:01 seelischExp $";2 version="$Id: presolve.lib,v 1.30 2009-04-14 12:00:14 Singular Exp $"; 3 3 category="Symbolic-numerical solving"; 4 4 info=" … … 190 190 if ( size(#)!=0 ) { n=#[1]; } 191 191 ideal maxi,rest = maxideal(1),0; 192 if ( n < nvars(BAS) ) 193 { 194 rest = maxi[n+1..nvars(BAS)]; 192 if ( n < nvars(BAS) ) 193 { 194 rest = maxi[n+1..nvars(BAS)]; 195 195 } 196 196 attrib(rest,"isSB",1); … … 200 200 // which do not contain elements not to be eliminated 201 201 202 //ideal id = interred(i); 203 //## gmg, gendert 9/2008: interred sehr lange z.B. bei Leonard1 in normal, 202 //ideal id = interred(i); 203 //## gmg, gendert 9/2008: interred sehr lange z.B. bei Leonard1 in normal, 204 204 //daher interred ersetzt durch: std nur auf linearpart angewendet 205 205 //Ordnung muss global sein, sonst egal (da Lin affin linear) … … 207 207 //--------------- replace ordering by dp if it is not global ----------------- 208 208 if ( ord_test(BAS) <= 0 ) 209 { 210 intvec V; 209 { 210 intvec V; 211 211 V[n]=0; V=V+1; //weights for dp ordering 212 212 gnirlist[3] = list("dp",V), g32; … … 215 215 ideal i = imap(BAS,i); 216 216 } 217 217 218 218 list Lin = linearpart(i); 219 219 ideal lin = std(Lin[1]); //SB of ideal generated by polys of i … … 231 231 //------------- check for special case of unit ideal and return --------------- 232 232 int check; 233 if( lin[1] == 1 ) 234 { 235 check = 1; 236 } 237 else 233 if( lin[1] == 1 ) 234 { 235 check = 1; 236 } 237 else 238 238 { 239 239 for (ii=1; ii<=size(id); ii++ ) 240 240 { 241 241 if ( id[ii] == 1 ) 242 { 242 { 243 243 check = 1; break; 244 244 } … … 356 356 /* 357 357 Alte Version mit interred: 358 // Then go to ring newBAS with ordering c,dp(n) and create a matrix with 359 // size(k1) colums and 2 rows, such that if [f1,f2] is a column of M then f1+f2 360 // is one of the polys of lin containing a pure degree 1 part and f1 is this 361 // part interreduce this matrix (i.e. Gauss elimination on linear part, with 358 // Then go to ring newBAS with ordering c,dp(n) and create a matrix with 359 // size(k1) colums and 2 rows, such that if [f1,f2] is a column of M then f1+f2 360 // is one of the polys of lin containing a pure degree 1 part and f1 is this 361 // part interreduce this matrix (i.e. Gauss elimination on linear part, with 362 362 // rest transformed accordingly). 363 363 //Ist jetzt durch direkte Substitution gemacht (schneller!) … … 365 365 //ideal k12 = k1,k2; 366 366 //matrix M = matrix(k12,2,kk); //degree 1 part is now in row 1 367 //M = interred(M); 368 //### interred zu teuer, muss nicht sein. Wenn interred angewendet 367 //M = interred(M); 368 //### interred zu teuer, muss nicht sein. Wenn interred angewendet 369 369 //werden soll, vorher in Ring mit Ordnung (c,dp) wechseln! 370 370 //Abfrage: if( ordstr(BAS) != "c,dp("+string(n)+")" ) … … 378 378 -z ergibt ich auch i[2]-z*i[3] mit option(redThrough) 379 379 statt interred kann man hier auch NF(i,i[3])+i[3] verwenden 380 hier lifert elimpart(i) 2 Substitutionen (x,y) elimpart(interred(i)) 380 hier lifert elimpart(i) 2 Substitutionen (x,y) elimpart(interred(i)) 381 381 aber 3 (x,y,z) 382 382 Da interred oder NF aber die Laenge der polys vergroessern kann, nicht gemacht … … 397 397 //since lin1 != 0 there are candidates for substituting variables 398 398 399 lin2 = lin - lin1; //difference as matrix 399 lin2 = lin - lin1; //difference as matrix 400 400 // rest of lin, part of pure degree 1 substracted from each generator of lin 401 401 … … 410 410 } 411 411 } 412 //Now each !=0 generator of lin2 contains only constant terms or terms of 412 //Now each !=0 generator of lin2 contains only constant terms or terms of 413 413 //degree >= 2, hence lin 2 can never be used for further substitutions 414 414 //We have: lin = ideal(matrix(k1)+matrix(k2)), lin2 415 415 416 ideal kin = matrix(k1)+matrix(k2); 416 ideal kin = matrix(k1)+matrix(k2); 417 417 //kin = polys of lin which contained a pure degree 1 part. 418 418 kin = simplify(kin,2); … … 421 421 int count=1; 422 422 while ( count != 0 ) 423 { 423 { 424 424 count = 0; 425 425 for ( ii=1; ii<=n; ii++ ) //start direct substitution of var(ii) … … 432 432 { 433 433 //we look for the shortest candidate to substitute var(ii) 434 if ( cand == 0 ) 435 { 434 if ( cand == 0 ) 435 { 436 436 cand = kin[kk]; //candidate for substituting var(ii) 437 } 437 } 438 438 else 439 439 { 440 if ( size(kin[kk]) < size(cand) ) 441 { 442 cand = kin[kk]; 440 if ( size(kin[kk]) < size(cand) ) 441 { 442 cand = kin[kk]; 443 443 } 444 444 } 445 445 } 446 } 446 } 447 447 if ( cand != 0 ) 448 448 { … … 452 452 neva[ii] = 0; 453 453 sub = sub+kip; //poly defining substituion 454 //## gmg: gendert 08/2008, map durch subst ersetzt 454 //## gmg: gendert 08/2008, map durch subst ersetzt 455 455 //(viel schneller) 456 456 vip = var(ii) - kip; //poly to be substituted … … 465 465 } 466 466 } 467 467 468 468 lin = kin+lin; 469 469 470 470 for( ii=1; ii<=size(lin); ii++ ) 471 471 { … … 473 473 } 474 474 475 for( ii=1; ii<=n; ii++ ) 475 for( ii=1; ii<=n; ii++ ) 476 476 { 477 477 for( kk=1; kk<=size(eva); kk++ ) -
Singular/LIB/random.lib
r2c3a5d r3360fb 1 1 //(GMG/BM, last modified 22.06.96) 2 2 /////////////////////////////////////////////////////////////////////////////// 3 version="$Id: random.lib,v 1.1 8 2009-03-30 18:36:35 motsakExp $";3 version="$Id: random.lib,v 1.19 2009-04-14 12:00:14 Singular Exp $"; 4 4 category="General purpose"; 5 5 info=" … … 154 154 proc sparseHomogIdeal (int k, int u, list #) 155 155 "USAGE: sparseid(k,u[,o,p,b]); k,u,o,p,b integers 156 RETURN: ideal having k homogeneous generators, each of random degree in the 157 interval [u,o], p percent of terms in degree d are 0, the remaining 158 have random coefficients in the interval [1,b], (default: o=u, p=75, 156 RETURN: ideal having k homogeneous generators, each of random degree in the 157 interval [u,o], p percent of terms in degree d are 0, the remaining 158 have random coefficients in the interval [1,b], (default: o=u, p=75, 159 159 b=30000) 160 160 EXAMPLE: example sparseid; shows an example … … 172 172 { 173 173 id = maxideal(random(u, o)); // monomial basis of some degree 174 m = sparsemat(size(id),1,p,b); // random coefficients 174 m = sparsemat(size(id),1,p,b); // random coefficients 175 175 i[ii] = (matrix(id)*m)[1,1]; 176 176 } -
Singular/LIB/ratgb.lib
r2c3a5d r3360fb 1 1 ////////////////////////////////////////////////////////////////////////////// 2 version="$Id: ratgb.lib,v 1.1 4 2009-02-21 15:26:42 levandovExp $";2 version="$Id: ratgb.lib,v 1.15 2009-04-14 12:00:15 Singular Exp $"; 3 3 category="Noncommutative"; 4 4 info=" … … 174 174 va = L3[w][2]; 175 175 for(z=1;z<=nvars(save)-is;z++) 176 { 177 vb[z] = va[is+z]; 176 { 177 vb[z] = va[is+z]; 178 178 } 179 179 tmp3[1] = "a"; … … 372 372 setring A; 373 373 IAppel1; 374 def F1 = ratstd(IAppel1,2); 375 lead(pGBid); 374 def F1 = ratstd(IAppel1,2); 375 lead(pGBid); 376 376 setring F1; rGBid; 377 377 } … … 385 385 setring A; 386 386 IAppel2; 387 def F1 = ratstd(IAppel2,2); 388 lead(pGBid); 387 def F1 = ratstd(IAppel2,2); 388 lead(pGBid); 389 389 setring F1; rGBid; 390 390 } … … 398 398 setring A; 399 399 IAppel4; 400 def F1 = ratstd(IAppel4,2); 401 lead(pGBid); 400 def F1 = ratstd(IAppel4,2); 401 lead(pGBid); 402 402 setring F1; rGBid; 403 403 } -
Singular/LIB/sing4ti2.lib
r2c3a5d r3360fb 1 1 /////////////////////////////////////////////////////////////////// 2 version="$Id: sing4ti2.lib,v 1. 2 2009-04-07 16:18:06 seelischExp $";2 version="$Id: sing4ti2.lib,v 1.3 2009-04-14 12:00:15 Singular Exp $"; 3 3 category="Commutative Algebra"; 4 4 info=" … … 12 12 @* the returned result 13 13 14 REQUIRES: External programs 4ti2, sed and awk to be installed 14 REQUIRES: External programs 4ti2, sed and awk to be installed 15 15 16 16 PROCEDURES: … … 88 88 89 89 //---------------------------------------------------------------------- 90 // calling 4ti2 and converting output 90 // calling 4ti2 and converting output 91 91 // Singular's string is too clumsy for this, hence we first prepare 92 92 // using standard unix commands 93 93 //---------------------------------------------------------------------- 94 94 j=system("sh","markov sing4ti2"); 95 j=system("sh","awk \'BEGIN{ORS=\",\";}{print $0;}\' sing4ti2.mar | sed s/[\\\ \\\t\\\v\\\f]/,/g | sed s/,+/,/g|sed s/,,/,/g|sed s/,,/,/g > sing4ti2.converted"); 95 j=system("sh","awk \'BEGIN{ORS=\",\";}{print $0;}\' sing4ti2.mar | sed s/[\\\ \\\t\\\v\\\f]/,/g | sed s/,+/,/g|sed s/,,/,/g|sed s/,,/,/g > sing4ti2.converted"); 96 96 if(!defined(keepfiles)) 97 97 { 98 98 j=system("sh",("rm -f sing4ti2.mar sing4ti2."+fileending)); 99 } 99 } 100 100 //---------------------------------------------------------------------- 101 101 // reading output of 4ti2 … … 115 115 { 116 116 if(erglist[2+(i-1)*erglist[2]+j]>=0) 117 { 117 { 118 118 //--- positive exponents 119 119 temppol1=temppol1*(var(j)^erglist[2+(i-1)*erglist[2]+j]); … … 125 125 } 126 126 } 127 toric=toric,temppol1-temppol2; 127 toric=toric,temppol1-temppol2; 128 128 } 129 129 //--- get rid of leading entry 0; … … 210 210 211 211 //---------------------------------------------------------------------- 212 // calling 4ti2 and converting output 212 // calling 4ti2 and converting output 213 213 // Singular's string is too clumsy for this, hence we first prepare 214 214 // using standard unix commands 215 215 //---------------------------------------------------------------------- 216 216 j=system("sh","graver sing4ti2"); 217 j=system("sh","awk \'BEGIN{ORS=\",\";}{print $0;}\' sing4ti2.gra | sed s/[\\\ \\\t\\\v\\\f]/,/g | sed s/,+/,/g |sed s/,,/,/g|sed s/,,/,/g > sing4ti2.converted"); 217 j=system("sh","awk \'BEGIN{ORS=\",\";}{print $0;}\' sing4ti2.gra | sed s/[\\\ \\\t\\\v\\\f]/,/g | sed s/,+/,/g |sed s/,,/,/g|sed s/,,/,/g > sing4ti2.converted"); 218 218 if(!defined(keepfiles)) 219 219 { 220 220 j=system("sh",("rm -f sing4ti2.gra sing4ti2."+fileending)); 221 } 221 } 222 222 //---------------------------------------------------------------------- 223 223 // reading output of 4ti2 … … 237 237 { 238 238 if(erglist[2+(i-1)*erglist[2]+j]>=0) 239 { 239 { 240 240 //--- positive exponents 241 241 temppol1=temppol1*(var(j)^erglist[2+(i-1)*erglist[2]+j]); … … 247 247 } 248 248 } 249 toric=toric,temppol1-temppol2; 249 toric=toric,temppol1-temppol2; 250 250 } 251 251 //--- get rid of leading entry 0; … … 272 272 @* - number of variables of basering equals number of columns of A 273 273 @* (for ker(A)) resp. of rows of A (for Im(A)) 274 CREATE: temporary files sing4ti2.mat, sing4ti2.lat, sing4ti2.mar 274 CREATE: temporary files sing4ti2.mat, sing4ti2.lat, sing4ti2.mar 275 275 @* in the current directory (I/O files for communication with 4ti2) 276 276 NOTE: input rules for 4ti2 also apply to input to this procedure … … 330 330 331 331 //---------------------------------------------------------------------- 332 // calling 4ti2 and converting output 332 // calling 4ti2 and converting output 333 333 // Singular's string is too clumsy for this, hence we first prepare 334 334 // using standard unix commands 335 335 //---------------------------------------------------------------------- 336 336 j=system("sh","hilbert sing4ti2"); 337 j=system("sh","awk \'BEGIN{ORS=\",\";}{print $0;}\' sing4ti2.hil | sed s/[\\\ \\\t\\\v\\\f]/,/g | sed s/,+/,/g |sed s/,,/,/g|sed s/,,/,/g > sing4ti2.converted"); 337 j=system("sh","awk \'BEGIN{ORS=\",\";}{print $0;}\' sing4ti2.hil | sed s/[\\\ \\\t\\\v\\\f]/,/g | sed s/,+/,/g |sed s/,,/,/g|sed s/,,/,/g > sing4ti2.converted"); 338 338 if(!defined(keepfiles)) 339 339 { 340 340 j=system("sh",("rm -f sing4ti2.hil sing4ti2."+fileending)); 341 } 341 } 342 342 //---------------------------------------------------------------------- 343 343 // reading output of 4ti2 … … 357 357 { 358 358 if(erglist[2+(i-1)*erglist[2]+j]>=0) 359 { 359 { 360 360 //--- positive exponents 361 361 temppol1=temppol1*(var(j)^erglist[2+(i-1)*erglist[2]+j]); … … 367 367 } 368 368 } 369 toric=toric,temppol1-temppol2; 369 toric=toric,temppol1-temppol2; 370 370 } 371 371 //--- get rid of leading entry 0; -
Singular/LIB/teachstd.lib
r2c3a5d r3360fb 2 2 //GMG, last modified 28.9.01 3 3 /////////////////////////////////////////////////////////////////////////////// 4 version="$Id: teachstd.lib,v 1.1 1 2009-04-06 12:39:02 seelischExp $";4 version="$Id: teachstd.lib,v 1.12 2009-04-14 12:00:15 Singular Exp $"; 5 5 category="Teaching"; 6 6 info=" … … 422 422 if( size(#) > 0 ) 423 423 {// "size(#): ", size(#); "typeof(#[1]): ", typeof(#[1]); 424 424 425 425 if( typeof(#[1]) == "int" ) 426 426 {// "#[1] = int ", #[1]; 427 427 if( #[1] > 0 ) 428 { 428 { 429 429 return(0); 430 430 } -
Singular/LIB/triang.lib
r2c3a5d r3360fb 1 1 //last change: 13.02.2001 (Eric Westenberger) 2 2 ////////////////////////////////////////////////////////////////////////////// 3 version="$Id: triang.lib,v 1.1 3 2009-04-06 09:17:01 seelischExp $";3 version="$Id: triang.lib,v 1.14 2009-04-14 12:00:15 Singular Exp $"; 4 4 category="Symbolic-numerical solving"; 5 5 info=" … … 629 629 If i = 2, then each polynomial of the triangular systems 630 630 is factorized. 631 NOTE: Algorithm of Moeller (see: Moeller, H.M.: On decomposing systems of 632 polynomial equations with finitely many solutions, Appl. Algebra Eng. 631 NOTE: Algorithm of Moeller (see: Moeller, H.M.: On decomposing systems of 632 polynomial equations with finitely many solutions, Appl. Algebra Eng. 633 633 Commun. Comput. 4, 217 - 230, 1993). 634 634 EXAMPLE: example triangM; shows an example -
Singular/LIB/tropical.lib
r2c3a5d r3360fb 1 version="$Id: tropical.lib,v 1.1 4 2009-04-08 12:42:19 seelischExp $";1 version="$Id: tropical.lib,v 1.15 2009-04-14 12:00:15 Singular Exp $"; 2 2 category="Tropical Geometry"; 3 3 info=" … … 8 8 @* Thomas Markwig, email: keilen@mathematik.uni-kl.de 9 9 10 WARNING: 10 WARNING: 11 11 - tropicalLifting will only work with LINUX and if in addition gfan is installed. 12 @*- drawTropicalCurve and drawTropicalNewtonSubdivision will only display the 13 @* tropical curve with LINUX and if in addition latex and kghostview 12 @*- drawTropicalCurve and drawTropicalNewtonSubdivision will only display the 13 @* tropical curve with LINUX and if in addition latex and kghostview 14 14 @* are installed. 15 @*- For tropicalLifting in the definition of the basering the parameter t 15 @*- For tropicalLifting in the definition of the basering the parameter t 16 16 @* from the Puiseux series field C{{t}} must be defined as a variable, 17 17 @* while for all other procedures it must be defined as a parameter. 18 18 19 19 THEORY: 20 Fix some base field K and a bunch of lattice points v0,...,vm in the integer 21 lattice Z^n, then this defines a toric variety as the closure of (K*)^n in 22 the projective space P^m, where the torus is embedded via the map sending a 20 Fix some base field K and a bunch of lattice points v0,...,vm in the integer 21 lattice Z^n, then this defines a toric variety as the closure of (K*)^n in 22 the projective space P^m, where the torus is embedded via the map sending a 23 23 point x in (K*)^n to the point (x^v0,...,x^vm). 24 The generic hyperplane sections are just the images of the hypersurfaces 24 The generic hyperplane sections are just the images of the hypersurfaces 25 25 in (K*)^n defined by the polynomials f=a0*x^v0+...+am*x^vm=0. Some properties 26 of these hypersurfaces can be studied via tropicalisation. 27 28 For this we suppose that K=C{{t}} is the field of Puiseux series over the 26 of these hypersurfaces can be studied via tropicalisation. 27 28 For this we suppose that K=C{{t}} is the field of Puiseux series over the 29 29 field of complex numbers (or any other field with a valuation into the real 30 numbers). One associates to the hypersurface given by f=a0*x^v0+...+am*x^vm 30 numbers). One associates to the hypersurface given by f=a0*x^v0+...+am*x^vm 31 31 the tropical hypersurface defined by the tropicalisation 32 trop(f)=min{val(a0)+<v0,x>,...,val(am)+<vm,x>}. 32 trop(f)=min{val(a0)+<v0,x>,...,val(am)+<vm,x>}. 33 33 Here, <v,x> denotes the standard scalar product of the integer vector v in Z^n 34 34 with the vector x=(x1,...,xn) of variables, so that trop(f) is a piecewise 35 linear function on R^n. The corner locus of this function (i.e. the points 36 at which the minimum is attained a least twice) is the tropical hypersurface 37 defined by trop(f). 38 The theorem of Newton-Kapranov states that this tropical hypersurface is 39 the same as if one computes pointwise the valuation of the hypersurface 40 given by f. The analogue holds true if one replaces one equation f by an 41 ideal I. A constructive proof of the theorem is given by an adapted 42 version of the Newton-Puiseux algorithm. The hard part is to find a point 43 in the variety over C{{t}} which corresponds to a given point in the 35 linear function on R^n. The corner locus of this function (i.e. the points 36 at which the minimum is attained a least twice) is the tropical hypersurface 37 defined by trop(f). 38 The theorem of Newton-Kapranov states that this tropical hypersurface is 39 the same as if one computes pointwise the valuation of the hypersurface 40 given by f. The analogue holds true if one replaces one equation f by an 41 ideal I. A constructive proof of the theorem is given by an adapted 42 version of the Newton-Puiseux algorithm. The hard part is to find a point 43 in the variety over C{{t}} which corresponds to a given point in the 44 44 tropical variety. 45 45 46 It is the purpose of this library to provide basic means to deal with 47 tropical varieties. Of course we cannot represent the field of Puiseux 48 series over C in its full strength, however, in order to compute interesting 49 examples it will be sufficient to replace the complex numbers C by the 50 rational numbers Q and to replace Puiseux series in t by rational functions 51 in t, i.e. we replace C{{t}} by Q(t), or sometimes even by Q[t]. 46 It is the purpose of this library to provide basic means to deal with 47 tropical varieties. Of course we cannot represent the field of Puiseux 48 series over C in its full strength, however, in order to compute interesting 49 examples it will be sufficient to replace the complex numbers C by the 50 rational numbers Q and to replace Puiseux series in t by rational functions 51 in t, i.e. we replace C{{t}} by Q(t), or sometimes even by Q[t]. 52 52 Note, that this in particular forbids rational exponents for the t's. 53 53 54 Moreover, in Singular no negative exponents of monomials are allowed, so 55 that the integer vectors vi will have to have non-negative entries. 56 Shifting all exponents by a fixed integer vector does not change the 57 tropicalisation nor does it change the toric variety. Thus this does not 54 Moreover, in Singular no negative exponents of monomials are allowed, so 55 that the integer vectors vi will have to have non-negative entries. 56 Shifting all exponents by a fixed integer vector does not change the 57 tropicalisation nor does it change the toric variety. Thus this does not 58 58 cause any restriction. 59 If, however, for some reason you prefer to work with general vi, then you 60 have to pass right away to the tropicalisation of the equations, whereever 59 If, however, for some reason you prefer to work with general vi, then you 60 have to pass right away to the tropicalisation of the equations, whereever 61 61 this is allowed -- these are linear polynomials where the constant coefficient 62 corresponds to the valuation of the original coefficient and where 63 the non-constant coefficient correspond to the exponents of the monomials, 64 thus they may be rational numbers respectively negative numbers: 62 corresponds to the valuation of the original coefficient and where 63 the non-constant coefficient correspond to the exponents of the monomials, 64 thus they may be rational numbers respectively negative numbers: 65 65 e.g. if f=t^{1/2}*x^{-2}*y^3+2t*x*y+4 then trop(f)=min{1/2-2x+3y,1+x+y,0}. 66 66 67 67 The main tools provided in this library are as follows: 68 @* - tropicalLifting implements the constructive proof of the Theorem of 69 Newton-Kapranov and constructs a point in the variety 70 over C{{t}} corresponding to a given point in the 71 corresponding tropical variety associated to an 72 ideal I; the generators of I have to be in the 73 polynomial ring Q[t,x1,...,xn] considered as a 74 subring of C{{t}}[x1,...,xn]; a solution will be 75 constructed up to given order; note that several 68 @* - tropicalLifting implements the constructive proof of the Theorem of 69 Newton-Kapranov and constructs a point in the variety 70 over C{{t}} corresponding to a given point in the 71 corresponding tropical variety associated to an 72 ideal I; the generators of I have to be in the 73 polynomial ring Q[t,x1,...,xn] considered as a 74 subring of C{{t}}[x1,...,xn]; a solution will be 75 constructed up to given order; note that several 76 76 field extensions of Q might be necessary throughout 77 77 the intermediate computations; the procedures use 78 78 the external program gfan 79 @* - drawTropicalCurve visualises a tropical plane curve either given by a 80 polynomial in Q(t)[x,y] or by a list of linear 81 polynomials of the form ax+by+c with a,b in Z and c 79 @* - drawTropicalCurve visualises a tropical plane curve either given by a 80 polynomial in Q(t)[x,y] or by a list of linear 81 polynomials of the form ax+by+c with a,b in Z and c 82 82 in Q; latex must be installed on your computer 83 @* - tropicalJInvariant computes the tropical j-invaiant of a tropical 83 @* - tropicalJInvariant computes the tropical j-invaiant of a tropical 84 84 elliptic curve 85 85 @* - jInvariant computes the j-invariant of an elliptic curve 86 @* - weierstrassForm computes the Weierstrass form of an elliptic curve 86 @* - weierstrassForm computes the Weierstrass form of an elliptic curve 87 87 88 88 PROCEDURES FOR TROPICAL LIFTING: 89 tropicalLifting computes a point in the tropical variety 89 tropicalLifting computes a point in the tropical variety 90 90 displayTropicalLifting displays the output of tropicalLifting 91 91 … … 104 104 tropicalise computes the tropicalisation of a polynomial 105 105 tropicaliseSet computes the tropicalisation several polynomials 106 tInitialForm computes the tInitial form of a poly in Q[t,x_1,...,x_n] 106 tInitialForm computes the tInitial form of a poly in Q[t,x_1,...,x_n] 107 107 tInitialIdeal computes the tInitial ideal of an ideal in Q[t,x_1,...,x_n] 108 initialForm computes the initial form of poly in Q[x_1,...,x_n] 109 initialIdeal computes the initial ideal of an ideal in Q[x_1,...,x_n] 108 initialForm computes the initial form of poly in Q[x_1,...,x_n] 109 initialIdeal computes the initial ideal of an ideal in Q[x_1,...,x_n] 110 110 111 111 PROCEDURES FOR LATEX CONVERSION: 112 112 texNumber outputs the texcommand for the leading coefficient of poly 113 texPolynomial outputs the texcommand for the polynomial poly 113 texPolynomial outputs the texcommand for the polynomial poly 114 114 texMatrix outputs the texcommand for the matrix 115 115 texDrawBasic embeds output of texDrawTropical in a texdraw environment 116 116 texDrawTropical computes the texdraw commands for a tropical curve 117 texDrawNewtonSubdivision computes texdraw commands for a Newton subdivision 117 texDrawNewtonSubdivision computes texdraw commands for a Newton subdivision 118 118 texDrawTriangulation computes texdraw commands for a triangulation 119 119 … … 121 121 radicalMemberShip checks radical membership 122 122 tInitialFormPar computes the t-initial form of poly in Q(t)[x_1,...,x_n] 123 tInitialFormParMax same as tInitialFormPar, but uses maximum 123 tInitialFormParMax same as tInitialFormPar, but uses maximum 124 124 solveTInitialFormPar displays approximated solution of a 0-dim ideal 125 125 detropicalise computes the detropicalisation of a linear form 126 dualConic computes the dual of an affine plane conic 126 dualConic computes the dual of an affine plane conic 127 127 parameterSubstitute substitutes in the poly the parameter t by t^N 128 128 tropicalSubst makes certain substitutions in a tropical polynomial … … 155 155 /// - eliminatecomponents 156 156 /// - findzerosAndBasictransform 157 /// - ordermaximalideals 157 /// - ordermaximalideals 158 158 /// - verticesTropicalCurve 159 159 /// - bunchOfLines … … 204 204 205 205 /////////////////////////////////////////////////////////////////////////////// 206 /// Procedures concerned with tropical parametrisation 206 /// Procedures concerned with tropical parametrisation 207 207 /////////////////////////////////////////////////////////////////////////////// 208 208 209 209 proc tropicalLifting (ideal i,intvec w,int ordnung,list #) 210 210 "USAGE: tropicalLifting(i,w,ord[,opt]); i ideal, w intvec, ord int, opt string 211 ASSUME: - i is an ideal in Q[t,x_1,...,x_n], w=(w_0,w_1,...,w_n) 212 and (w_1/w_0,...,w_n/w_0) is in the tropical variety of i, 213 and ord is the order up to which a point in V(i) over Q{{t}} 214 lying over (w_1/w_0,...,w_n/w_0) shall be computed; 211 ASSUME: - i is an ideal in Q[t,x_1,...,x_n], w=(w_0,w_1,...,w_n) 212 and (w_1/w_0,...,w_n/w_0) is in the tropical variety of i, 213 and ord is the order up to which a point in V(i) over Q{{t}} 214 lying over (w_1/w_0,...,w_n/w_0) shall be computed; 215 215 w_0 may NOT be ZERO 216 @* - the basering should not have any parameters on its own 217 and it should have a global monomial ordering, 216 @* - the basering should not have any parameters on its own 217 and it should have a global monomial ordering, 218 218 e.g. ring r=0,(t,x(1..n)),dp; 219 @* - the first variable of the basering will be treated as the 219 @* - the first variable of the basering will be treated as the 220 220 parameter t in the Puiseux series field 221 @* - the optional parameter opt should be one or more strings among 221 @* - the optional parameter opt should be one or more strings among 222 222 the following: 223 223 @* 'isZeroDimensional' : the dimension i is zero (not to be checked); 224 @* 'isPrime' : the ideal is prime over Q(t)[x_1,...,x_n] 224 @* 'isPrime' : the ideal is prime over Q(t)[x_1,...,x_n] 225 225 (not to be checked); 226 @* 'isInTrop' : (w_1/w_0,...,w_n/w_0) is in the tropical 226 @* 'isInTrop' : (w_1/w_0,...,w_n/w_0) is in the tropical 227 227 variety (not to be checked); 228 @* 'oldGfan' : uses gfan version 0.2.1 or less 229 @* 'findAll' : find all solutions of a zero-dimensional 228 @* 'oldGfan' : uses gfan version 0.2.1 or less 229 @* 'findAll' : find all solutions of a zero-dimensional 230 230 ideal over (w_1/w_0,...,w_n/w_0) 231 231 @* 'noAbs' : do NOT use absolute primary decomposition … … 233 233 RETURN: IF THE OPTION 'findAll' WAS NOT SET THEN: 234 234 @* list, containing one lifting of the given point (w_1/w_0,...,w_n/w_0) 235 in the tropical variety of i to a point in V(i) over Puiseux 235 in the tropical variety of i to a point in V(i) over Puiseux 236 236 series field up to the first ord terms; more precisely: 237 237 @* IF THE OPTION 'noAbs' WAS NOT SET, THEN: … … 248 248 @* IF THE OPITON 'findAll' WAS SET, THEN: 249 249 @* list, containing ALL liftings of the given point ((w_1/w_0,...,w_n/w_0) 250 in the tropical variety of i to a point in V(i) over Puiseux 251 series field up to the first ord terms, if the ideal is 250 in the tropical variety of i to a point in V(i) over Puiseux 251 series field up to the first ord terms, if the ideal is 252 252 zero-dimensional over Q{{t}}; 253 more precisely, each entry of the list is a list l as computed 253 more precisely, each entry of the list is a list l as computed 254 254 if 'findAll' was NOT set 255 255 @* WE NOW DESCRIBE THE LIST ENTRIES IF 'findAll' WAS NOT SET: 256 @* - the ring l[1] contains an ideal LIFT, which contains 256 @* - the ring l[1] contains an ideal LIFT, which contains 257 257 a point in V(i) lying over w up to the first ord terms; 258 @* - and if the integer l[2] is N then t has to be replaced by t^1/N 258 @* - and if the integer l[2] is N then t has to be replaced by t^1/N 259 259 in the lift, or alternatively replace t by t^N in the defining ideal 260 @* - if the k+1st entry of l[3] is non-zero, then the kth component of 261 LIFT has to be multiplied t^(-l[3][k]/l[3][1]) AFTER substituting t 260 @* - if the k+1st entry of l[3] is non-zero, then the kth component of 261 LIFT has to be multiplied t^(-l[3][k]/l[3][1]) AFTER substituting t 262 262 by t^1/N 263 @* - unless the option 'noResubst' was set, the kth entry of list l[4] 263 @* - unless the option 'noResubst' was set, the kth entry of list l[4] 264 264 is a string which represents the kth generator of 265 the ideal i where the coordinates have been replaced by the result 265 the ideal i where the coordinates have been replaced by the result 266 266 of the lift; 267 the t-order of the kth entry should in principle be larger than the 267 the t-order of the kth entry should in principle be larger than the 268 268 t-degree of LIFT 269 @* - if the option 'noAbs' was set, then the string in l[5] defines 270 a maximal ideal in the field Q[X(1),...,X(k)], where X(1),...,X(k) 269 @* - if the option 'noAbs' was set, then the string in l[5] defines 270 a maximal ideal in the field Q[X(1),...,X(k)], where X(1),...,X(k) 271 271 are the parameters of the ring in l[1]; 272 the basefield of the ring in l[1] should be considered modulo this 272 the basefield of the ring in l[1] should be considered modulo this 273 273 ideal 274 REMARK: - it is best to use the procedure displayTropicalLifting to 274 REMARK: - it is best to use the procedure displayTropicalLifting to 275 275 display the result 276 276 @* - the option 'findAll' cannot be used if 'noAbs' is set 277 277 @* - if the parameter 'findAll' is set AND the ideal i is zero-dimensional 278 in Q{{t}}[x_1,...,x_n] then ALL points in V(i) lying over w are 278 in Q{{t}}[x_1,...,x_n] then ALL points in V(i) lying over w are 279 279 computed up to order ord; if the ideal is not-zero dimenisonal, then 280 280 only the points in the ideal after cutting down to dimension zero 281 281 will be computed 282 282 @* - the procedure requires that the program GFAN is installed on your 283 computer; if you have GFAN version less than 0.3.0 then you must 283 computer; if you have GFAN version less than 0.3.0 then you must 284 284 use the optional parameter 'oldGfan' 285 @* - the procedure requires the Singular procedure absPrimdecGTZ to be 285 @* - the procedure requires the Singular procedure absPrimdecGTZ to be 286 286 present in the package primdec.lib, unless the option 'noAbs' is set; 287 but even if absPrimdecGTZ is present it might be necessary to set 288 the option 'noAbs' in order to avoid the costly absolute primary 289 decomposition; the side effect is that the field extension which is 287 but even if absPrimdecGTZ is present it might be necessary to set 288 the option 'noAbs' in order to avoid the costly absolute primary 289 decomposition; the side effect is that the field extension which is 290 290 computed throughout the recursion might need more than one 291 291 parameter to be described 292 292 @* - since Q is infinite, the procedure finishes with probability one 293 @* - you can call the procedure with Z/pZ as base field instead of Q, 293 @* - you can call the procedure with Z/pZ as base field instead of Q, 294 294 but there are some problems you should be aware of: 295 @* + the Puiseux series field over the algebraic closure of Z/pZ is 296 NOT algebraicall closed, and thus there may not exist a point in 297 V(i) over the Puiseux series field with the desired valuation; 295 @* + the Puiseux series field over the algebraic closure of Z/pZ is 296 NOT algebraicall closed, and thus there may not exist a point in 297 V(i) over the Puiseux series field with the desired valuation; 298 298 so there is no chance that the procedure produced a sensible output 299 - e.g. if i=tx^p-tx-1 300 @* + if the dimension of i over Z/pZ(t) is not zero the process of 301 reduction to zero might not work if the characteristic is small 299 - e.g. if i=tx^p-tx-1 300 @* + if the dimension of i over Z/pZ(t) is not zero the process of 301 reduction to zero might not work if the characteristic is small 302 302 and you are unlucky 303 @* + the option 'noAbs' has to be used since absolute primary 303 @* + the option 'noAbs' has to be used since absolute primary 304 304 decomposition in Singular only works in characteristic zero 305 @* - the basefield should either be Q or Z/pZ for some prime p; 306 field extensions will be computed if necessary; if you need 307 parameters or field extensions from the beginning they should 308 rather be simulated as variables possibly adding their relations to 305 @* - the basefield should either be Q or Z/pZ for some prime p; 306 field extensions will be computed if necessary; if you need 307 parameters or field extensions from the beginning they should 308 rather be simulated as variables possibly adding their relations to 309 309 the ideal; the weights for the additional variables should be zero 310 310 EXAMPLE: example tropicalLifting; shows an example" … … 357 357 noabs=1; 358 358 } 359 // this option is not documented -- it prevents the execution of gfan and 360 // just asks for wneu to be inserted -- it can be used to check problems 361 // with the precedure without calling gfan, if wneu is know from previous 359 // this option is not documented -- it prevents the execution of gfan and 360 // just asks for wneu to be inserted -- it can be used to check problems 361 // with the precedure without calling gfan, if wneu is know from previous 362 362 // computations 363 363 if (#[j]=="noGfan") … … 370 370 } 371 371 } 372 // if the basering has characteristic not equal to zero, 372 // if the basering has characteristic not equal to zero, 373 373 // then absolute factorisation 374 374 // is not available, and thus we need the option noAbs … … 384 384 { 385 385 Error("The first coordinate of your input w must be NON-ZERO, since it is a DENOMINATOR!"); 386 } 386 } 387 387 // if w_0<0, then replace w by -w, so that the "denominator" w_0 is positive 388 388 if (w[1]<0) … … 391 391 } 392 392 intvec prew=w; // stores w for later reference 393 // for our computations, w[1] represents the weight of t and this 393 // for our computations, w[1] represents the weight of t and this 394 394 // should be -w_0 !!! 395 395 w[1]=-w[1]; … … 401 401 w[1]=-1; 402 402 } 403 // if some entry of w is positive, we have to make a transformation, 403 // if some entry of w is positive, we have to make a transformation, 404 404 // which moves it to something non-positive 405 405 for (j=2;j<=nvars(basering);j++) … … 427 427 { 428 428 variablen=variablen+var(j); 429 } 429 } 430 430 map GRUNDPHI=BASERING,t,variablen; 431 431 ideal i=GRUNDPHI(i); 432 // compute the initial ideal of i and test if w is in the tropical 433 // variety of i 432 // compute the initial ideal of i and test if w is in the tropical 433 // variety of i 434 434 // - the last entry 1 only means that t is the last variable in the ring 435 435 ideal ini=tInitialIdeal(i,w,1); 436 436 if (isintrop==0) // test if w is in trop(i) only if isInTrop has not been set 437 { 437 { 438 438 poly product=1; 439 439 for (j=1;j<=nvars(basering)-1;j++) … … 453 453 int dd=dim(i); 454 454 setring GRUNDRING; 455 // if the dimension is not zero, we cut the ideal down to dimension zero 455 // if the dimension is not zero, we cut the ideal down to dimension zero 456 456 // and compute the 457 457 // t-initial ideal of the new ideal at the same time 458 458 if(dd!=0) 459 459 { 460 // the procedurce cutdown computes a new ring, in which there lives a 460 // the procedurce cutdown computes a new ring, in which there lives a 461 461 // zero-dimensional 462 // ideal which has been computed by cutting down the input with 462 // ideal which has been computed by cutting down the input with 463 463 // generic linear forms 464 // of the type x_i1-p_1,...,x_id-p_d for some polynomials 465 // p_1,...,p_d not depending 466 // on the variables x_i1,...,x_id; that way we have reduced 464 // of the type x_i1-p_1,...,x_id-p_d for some polynomials 465 // p_1,...,p_d not depending 466 // on the variables x_i1,...,x_id; that way we have reduced 467 467 // the number of variables by dd !!! 468 // the new zero-dimensional ideal is called i, its t-initial 468 // the new zero-dimensional ideal is called i, its t-initial 469 469 // ideal (with respect to 470 // the new w=CUTDOWN[2]) is ini, and finally there is a list 471 // repl in the ring 470 // the new w=CUTDOWN[2]) is ini, and finally there is a list 471 // repl in the ring 472 472 // which contains at the polynomial p_j at position i_j and 473 473 //a zero otherwise; … … 492 492 list liftrings; // will contain the final result 493 493 // if the procedure is called without 'findAll' then it may happen, that no 494 // proper solution is found when dd>0; in that case we have 494 // proper solution is found when dd>0; in that case we have 495 495 // to start all over again; 496 496 // this is controlled by the while-loop … … 508 508 // compute the liftrings by resubstitution 509 509 kk=1; // counts the liftrings 510 int isgood; // test in the non-zerodimensional case 510 int isgood; // test in the non-zerodimensional case 511 511 // if the result has the correct valuation 512 512 for (jj=1;jj<=size(TP);jj++) 513 513 { 514 // the list TP contains as a first entry the ring over which the 515 // tropical parametrisation 514 // the list TP contains as a first entry the ring over which the 515 // tropical parametrisation 516 516 // of the (possibly cutdown ideal) i lives 517 517 def LIFTRING=TP[jj][1]; 518 // if the dimension of i originally was not zero, 518 // if the dimension of i originally was not zero, 519 519 // then we have to fill in the missing 520 520 // parts of the parametrisation 521 521 if (dd!=0) 522 522 { 523 // we need a ring where the parameters X_1,...,X_k 523 // we need a ring where the parameters X_1,...,X_k 524 524 // from LIFTRING are present, 525 525 // and where also the variables of CUTDOWNRING live 526 526 execute("ring REPLACEMENTRING=("+charstr(LIFTRING)+"),("+varstr(CUTDOWNRING)+"),dp;"); 527 list repl=imap(CUTDOWNRING,repl); // get the replacement rules 527 list repl=imap(CUTDOWNRING,repl); // get the replacement rules 528 528 // from CUTDOWNRING 529 ideal PARA=imap(LIFTRING,PARA); // get the zero-dim. parametrisatio 529 ideal PARA=imap(LIFTRING,PARA); // get the zero-dim. parametrisatio 530 530 // from LIFTRING 531 531 // compute the lift of the solution of the original ideal i … … 533 533 k=1; 534 534 // the lift has as many components as GRUNDRING has variables!=t 535 for (j=1;j<=nvars(GRUNDRING)-1;j++) 535 for (j=1;j<=nvars(GRUNDRING)-1;j++) 536 536 { 537 537 // if repl[j]=0, then the corresponding variable was not eliminated 538 if (repl[j]==0) 538 if (repl[j]==0) 539 539 { 540 LIFT[j]=PARA[k]; // thus the lift has been 540 LIFT[j]=PARA[k]; // thus the lift has been 541 541 // computed by tropicalparametrise 542 542 k++; // k checks how many entries of PARA have already been used … … 544 544 else // if repl[j]!=0, repl[j] contains replacement rule for the lift 545 545 { 546 LIFT[j]=repl[j]; // we still have to replace the vars 546 LIFT[j]=repl[j]; // we still have to replace the vars 547 547 // in repl[j] by the corresp. entries of PARA 548 548 // replace all variables!=t (from CUTDOWNRING) 549 for (l=1;l<=nvars(CUTDOWNRING)-1;l++) 549 for (l=1;l<=nvars(CUTDOWNRING)-1;l++) 550 550 { 551 551 // substitute the kth variable by PARA[k] 552 LIFT[j]=subst(LIFT[j],var(l),PARA[l]); 552 LIFT[j]=subst(LIFT[j],var(l),PARA[l]); 553 553 } 554 554 } 555 555 } 556 556 setring LIFTRING; 557 ideal LIFT=imap(REPLACEMENTRING,LIFT); 558 // test now if the LIFT has the correct valuation !!! 559 // note: it may happen, that when resubstituting PARA into 557 ideal LIFT=imap(REPLACEMENTRING,LIFT); 558 // test now if the LIFT has the correct valuation !!! 559 // note: it may happen, that when resubstituting PARA into 560 560 // the replacement rules 561 // there occured some unexpected cancellation; 561 // there occured some unexpected cancellation; 562 562 // we only know that for SOME 563 // solution of the zero-dimensional reduction NO 564 // canellation will occur, 565 // but for others this may very well happen; 563 // solution of the zero-dimensional reduction NO 564 // canellation will occur, 565 // but for others this may very well happen; 566 566 // this in particular means that 567 // we possibly MUST compute all zero-dimensional 567 // we possibly MUST compute all zero-dimensional 568 568 // solutions when cutting down! 569 569 intvec testw=precutdownw[1]; … … 589 589 kill PARA; 590 590 // only if LIFT has the right valuation we have to do something 591 if (isgood==1) 592 { 593 // it remains to reverse the original substitutions, 591 if (isgood==1) 592 { 593 // it remains to reverse the original substitutions, 594 594 // where appropriate !!! 595 // if some entry of the original w was positive, 595 // if some entry of the original w was positive, 596 596 // we replace the corresponding 597 597 // variable x_i by t^-w[i]*x_i, so we must now replace … … 610 610 */ 611 611 // if LIFTRING contains a parameter @a, change it to a 612 if ((noabs==0) and (defined(@a)==-1)) 612 if ((noabs==0) and (defined(@a)==-1)) 613 613 { 614 // pass first to a ring where a and @a 614 // pass first to a ring where a and @a 615 615 // are variables in order to use maps 616 616 poly mp=minpoly; … … 621 621 // replace @a by a in minpoly and in LIFT 622 622 map phi=INTERRING,t,a,a; 623 mp=phi(mp); 623 mp=phi(mp); 624 624 LIFT=phi(LIFT); 625 625 // pass now to a ring whithout @a and with a as parameter … … 628 628 ideal LIFT=imap(INTERRING,LIFT); 629 629 kill INTERRING; 630 } 630 } 631 631 // then export LIFT 632 export(LIFT); 632 export(LIFT); 633 633 // test the result by resubstitution 634 setring GRUNDRING; 634 setring GRUNDRING; 635 635 list resubst; 636 636 if (noresubst==0) … … 641 641 } 642 642 else 643 { 643 { 644 644 resubst=tropicalliftingresubstitute(substitute(i,t,t^(TP[jj][2])),list(LIFTRING),N*TP[jj][2]); 645 645 } 646 646 } 647 647 setring BASERING; 648 // Finally, if t has been replaced by t^N, then we have to change the 648 // Finally, if t has been replaced by t^N, then we have to change the 649 649 // third entry of TP by multiplying by N. 650 650 if (noabs==1) … … 662 662 kill LIFTRING; 663 663 } 664 // if dd!=0 and the procedure was called without the 664 // if dd!=0 and the procedure was called without the 665 665 // option findAll, then it might very well 666 // be the case that no solution is found, since 666 // be the case that no solution is found, since 667 667 // only one solution for the zero-dimensional 668 // reduction was computed and this one might have 668 // reduction was computed and this one might have 669 669 // had cancellations when resubstituting; 670 670 // if so we have to restart the process with the option findAll … … 674 674 "The procedure will be restarted with the option 'findAll'."; 675 675 "Go on by hitting RETURN!"; 676 findall=1; 676 findall=1; 677 677 noabs=0; 678 678 setring CUTDOWNRING; … … 680 680 "i";i; 681 681 "ini";tInitialIdeal(i,w,1); 682 682 683 683 /* 684 684 setring GRUNDRING; … … 698 698 } 699 699 } 700 // if internally the option findall was set, then return 700 // if internally the option findall was set, then return 701 701 // only the first solution 702 702 if (defined(hadproblems)!=0) … … 707 707 if (voice+printlevel>=2) 708 708 { 709 709 710 710 "The procedure has created a list of lists. The jth entry of this list 711 711 contains a ring, an integer and an intvec. 712 712 In this ring lives an ideal representing the wanted lifting, 713 713 if the integer is N then in the parametrisation t has to be replaced by t^1/N, 714 and if the ith component of the intvec is w[i] then the ith component in LIFT 714 and if the ith component of the intvec is w[i] then the ith component in LIFT 715 715 should be multiplied by t^-w[i]/N in order to get the parametrisation. 716 716 717 717 Suppose your list has the name L, then you can access the 1st ring via: 718 718 "; 719 719 if (findall==1) 720 720 { 721 "def LIFTRing=L[1][1]; setring LIFTRing; LIFT; 721 "def LIFTRing=L[1][1]; setring LIFTRing; LIFT; 722 722 "; 723 723 } 724 724 else 725 725 { 726 "def LIFTRing=L[1]; setring LIFTRing; LIFT; 726 "def LIFTRing=L[1]; setring LIFTRing; LIFT; 727 727 "; 728 } 728 } 729 729 } 730 730 if (findall==1) // if all solutions have been computed, return a list of lists … … 751 751 def LIFTRing=LIST[1]; 752 752 setring LIFTRing; 753 // LIFT contains the first 4 terms of a point in the variety of i 753 // LIFT contains the first 4 terms of a point in the variety of i 754 754 // over the Puiseux series field C{{t}} whose order is -w[1]/w[0]=1 755 755 LIFT; … … 779 779 // NOTE: since the last component of v is positive, the lifting 780 780 // must start with a negative power of t, which in Singular 781 // is not allowed for a variable. 781 // is not allowed for a variable. 782 782 def LIFTRing3=LIST[1]; 783 783 setring LIFTRing3; … … 834 834 string Kstring="Z/"+string(char(LIFTRing))+"Z"; 835 835 } 836 // this means that tropicalLifting was called with 836 // this means that tropicalLifting was called with 837 837 // absolute primary decomposition 838 if (size(troplift)==4) 839 { 838 if (size(troplift)==4) 839 { 840 840 setring LIFTRing; 841 841 "The lifting of the point in the tropical variety lives in the ring"; 842 842 if ((size(LIFTpar)==0) and (N==1)) 843 843 { 844 Kstring+"[[t]]"; 844 Kstring+"[[t]]"; 845 845 } 846 846 if ((size(LIFTpar)==0) and (N!=1)) 847 847 { 848 Kstring+"[[t^(1/"+string(N)+")]]"; 848 Kstring+"[[t^(1/"+string(N)+")]]"; 849 849 } 850 850 if ((size(LIFTpar)!=0) and (N!=1)) 851 { 852 Kstring+"["+LIFTpar+"]/"+string(minpoly)+"[[t^(1/"+string(N)+")]]"; 851 { 852 Kstring+"["+LIFTpar+"]/"+string(minpoly)+"[[t^(1/"+string(N)+")]]"; 853 853 } 854 854 if ((size(LIFTpar)!=0) and (N==1)) 855 { 856 Kstring+"["+LIFTpar+"]/"+string(minpoly)+"[[t]]"; 855 { 856 Kstring+"["+LIFTpar+"]/"+string(minpoly)+"[[t]]"; 857 857 } 858 858 } … … 871 871 } 872 872 if ((size(LIFTpar)!=0) and (N!=1)) 873 { 874 Kstring+"["+LIFTpar+"]/M[[t^(1/"+string(N)+")]]"; 873 { 874 Kstring+"["+LIFTpar+"]/M[[t^(1/"+string(N)+")]]"; 875 875 "where M is the maximal ideal"; 876 876 "M=<"+m+">"; 877 877 } 878 878 if ((size(LIFTpar)!=0) and (N==1)) 879 { 880 Kstring+"["+LIFTpar+"]/M[[t]]"; 879 { 880 Kstring+"["+LIFTpar+"]/M[[t]]"; 881 881 "where M is the maximal ideal"; 882 882 "M=<"+m+">"; 883 } 883 } 884 884 } 885 885 ""; … … 908 908 } 909 909 } 910 } 910 } 911 911 } 912 912 example … … 922 922 923 923 /////////////////////////////////////////////////////////////////////////////// 924 /// Procedures concerned with drawing a tropical curve or a Newton subdivision 924 /// Procedures concerned with drawing a tropical curve or a Newton subdivision 925 925 /////////////////////////////////////////////////////////////////////////////// 926 926 927 927 proc tropicalCurve (def tp,list #) 928 928 "USAGE: tropicalCurve(tp[,#]); tp list, # optional list 929 ASSUME: tp is list of linear polynomials of the form ax+by+c 930 with integers a, b and a rational number c representing 929 ASSUME: tp is list of linear polynomials of the form ax+by+c 930 with integers a, b and a rational number c representing 931 931 a tropical Laurent polynomial defining a tropical plane curve; 932 alternatively tp can be a polynomial in Q(t)[x,y] defining a 933 tropical plane curve via the valuation map; 934 the basering must have a global monomial ordering, 932 alternatively tp can be a polynomial in Q(t)[x,y] defining a 933 tropical plane curve via the valuation map; 934 the basering must have a global monomial ordering, 935 935 two variables and up to one parameter! 936 RETURN: list, each entry i=1,...,size(l)-1 corresponds to a vertex 936 RETURN: list, each entry i=1,...,size(l)-1 corresponds to a vertex 937 937 in the tropical plane curve defined by tp 938 938 l[i][1] = x-coordinate of the ith vertex 939 939 l[i][2] = y-coordinate of the ith vertex 940 l[i][3] = intmat, if j is an entry in the first row 940 l[i][3] = intmat, if j is an entry in the first row 941 941 of intmat then the ith vertex of 942 the tropical curve is connected to the 942 the tropical curve is connected to the 943 943 jth vertex with multiplicity given 944 944 by the corresponding entry in the second row 945 l[i][4] = list of lists, the first entry of a list is 946 a primitive integer vector defining the direction 947 of an unbounded edge emerging from the ith vertex 948 of the graph, the corresponding second entry in 945 l[i][4] = list of lists, the first entry of a list is 946 a primitive integer vector defining the direction 947 of an unbounded edge emerging from the ith vertex 948 of the graph, the corresponding second entry in 949 949 the list is the multiplicity of the unbounded edge 950 l[i][5] = a polynomial whose monomials mark the vertices 950 l[i][5] = a polynomial whose monomials mark the vertices 951 951 in the Newton polygon corresponding to the entries 952 in tp which take the common minimum at the ith 953 vertex -- if some coefficient a or b of the 954 linear polynomials in the input was negative, 952 in tp which take the common minimum at the ith 953 vertex -- if some coefficient a or b of the 954 linear polynomials in the input was negative, 955 955 then each monomial has to be shifted by 956 956 the values in l[size(l)][3] 957 l[size(l)][1] = list, the entries describe the boundary 957 l[size(l)][1] = list, the entries describe the boundary 958 958 points of the Newton subdivision 959 l[size(l)][2] = list, the entries are pairs of integer 959 l[size(l)][2] = list, the entries are pairs of integer 960 960 vectors defining an interior 961 961 edge of the Newton subdivision 962 l[size(l)][3] = intvec, the monmials occuring in l[i][5] 963 have to be shifted by this vector 964 in order to represent marked 962 l[size(l)][3] = intvec, the monmials occuring in l[i][5] 963 have to be shifted by this vector 964 in order to represent marked 965 965 vertices in the Newton polygon 966 NOTE: here the tropical polynomial is supposed to be the MINIMUM 967 of the linear forms in tp, unless the optional input #[1] 966 NOTE: here the tropical polynomial is supposed to be the MINIMUM 967 of the linear forms in tp, unless the optional input #[1] 968 968 is the string 'max' 969 969 EXAMPLE: example tropicalCurve; shows an example" … … 978 978 ERROR("The basering should have a global monomial ordering, e.g. ring r=(0,t),(x,y),dp;"); 979 979 } 980 // if you insert a single polynomial instead of an ideal 980 // if you insert a single polynomial instead of an ideal 981 981 // representing a tropicalised polynomial, 982 // then we compute first the tropicalisation of this polynomial 982 // then we compute first the tropicalisation of this polynomial 983 983 // -- this feature is not documented in the above help string 984 984 if (typeof(tp)=="poly") 985 985 { 986 // exclude the case that the basering has not precisely 986 // exclude the case that the basering has not precisely 987 987 // one parameter and two indeterminates 988 988 if ((npars(basering)!=1) or (nvars(basering)!=2)) 989 989 { 990 ERROR("The basering should have precisely one parameter and two indeterminates!"); 990 ERROR("The basering should have precisely one parameter and two indeterminates!"); 991 991 } 992 992 poly f=tp; … … 997 997 if (nvars(basering) != 2) 998 998 { 999 ERROR("The basering should have precisely two indeterminates!"); 1000 } 1001 // -1) Exclude the pathological case that the defining 999 ERROR("The basering should have precisely two indeterminates!"); 1000 } 1001 // -1) Exclude the pathological case that the defining 1002 1002 // tropical polynomial has only one term, 1003 1003 // so that the tropical variety is not defined. … … 1007 1007 intmat M[2][1]=0,0; 1008 1008 return(list(list(0,0,M,list(),detropicalise(tp[1])),list(list(leadexp(detropicalise(tp[1]))),list()))); 1009 } 1010 // 0) If the input was a list of linear polynomials, 1009 } 1010 // 0) If the input was a list of linear polynomials, 1011 1011 // then some coefficient of x or y can be negative, 1012 // i.e. the input corresponds to the tropical curve 1012 // i.e. the input corresponds to the tropical curve 1013 1013 // of a Laurent polynomial. In that case we should 1014 // add some ax+by, so that all coefficients are positive. 1014 // add some ax+by, so that all coefficients are positive. 1015 1015 // This does not change the tropical curve. 1016 // however, we have to save (a,b), since the Newton 1016 // however, we have to save (a,b), since the Newton 1017 1017 // polygone has to be shifted by (-a,-b). 1018 1018 poly aa,bb; // koeffizienten … … 1026 1026 { 1027 1027 bb=koeffizienten(tp[i],2); 1028 } 1028 } 1029 1029 } 1030 1030 if ((aa!=0) or (bb!=0)) … … 1035 1035 } 1036 1036 } 1037 // 1) compute the vertices of the tropical curve 1037 // 1) compute the vertices of the tropical curve 1038 1038 // defined by tp and the Newton subdivision 1039 1039 list vtp=verticesTropicalCurve(tp,#); 1040 // if vtp is empty, then the Newton polygone is just 1040 // if vtp is empty, then the Newton polygone is just 1041 1041 // a line segment and constitutes a bunch of lines 1042 1042 // which can be computed by bunchOfLines … … 1045 1045 return(bunchOfLines(tp)); 1046 1046 } 1047 // 2) store all vertices belonging to the ith part of the 1047 // 2) store all vertices belonging to the ith part of the 1048 1048 // Newton subdivision in the list vtp[i] as 4th entry, 1049 1049 // and store those, which are not corners of the ith subdivision polygon 1050 1050 // in vtp[i][6] 1051 poly nwt; 1051 poly nwt; 1052 1052 list boundaryNSD; // stores the boundary of a Newton subdivision 1053 intmat zwsp[2][1]; // used for intermediate storage 1053 intmat zwsp[2][1]; // used for intermediate storage 1054 1054 for (i=1;i<=size(vtp);i++) 1055 1055 { 1056 1056 k=1; 1057 nwt=vtp[i][3]; // the polynomial representing the 1057 nwt=vtp[i][3]; // the polynomial representing the 1058 1058 // ith part of the Newton subdivision 1059 // store the vertices of the ith part of the 1059 // store the vertices of the ith part of the 1060 1060 // Newton subdivision in the list newton 1061 list newton; 1061 list newton; 1062 1062 while (nwt!=0) 1063 1063 { … … 1066 1066 k++; 1067 1067 } 1068 boundaryNSD=findOrientedBoundary(newton);// a list of the vertices 1069 // of the Newton subdivision 1070 // as integer vectors (only those 1071 // on the boundary, and oriented 1068 boundaryNSD=findOrientedBoundary(newton);// a list of the vertices 1069 // of the Newton subdivision 1070 // as integer vectors (only those 1071 // on the boundary, and oriented 1072 1072 // clockwise) 1073 1073 vtp[i][4]=boundaryNSD[1]; 1074 1074 vtp[i][5]=boundaryNSD[2]; 1075 vtp[i][6]=zwsp; // the entries of the first row will denote to which 1076 // vertex the ith one is connected 1077 // and the entries of the second row will denote 1075 vtp[i][6]=zwsp; // the entries of the first row will denote to which 1076 // vertex the ith one is connected 1077 // and the entries of the second row will denote 1078 1078 //with which multiplicity 1079 1079 kill newton; // we kill the superflous list 1080 1080 } 1081 // 3) Next we build for each part of the Newton 1081 // 3) Next we build for each part of the Newton 1082 1082 // subdivision the list of all pairs of vertices on the 1083 1083 // boundary, which are involved, including those which are not corners … … 1096 1096 kill ipairs; 1097 1097 } 1098 // 4) Check for all pairs of verticies in the Newton diagram if they 1098 // 4) Check for all pairs of verticies in the Newton diagram if they 1099 1099 // occur in two different parts of the Newton subdivision 1100 int deleted; // if a pair occurs in two NSD, it can be removed 1100 int deleted; // if a pair occurs in two NSD, it can be removed 1101 1101 // from both - deleted is then set to 1 1102 list inneredges; // contains the list of all pairs contained in two NSD 1102 list inneredges; // contains the list of all pairs contained in two NSD 1103 1103 // - these are inner the edges of NSD 1104 1104 int ggt; 1105 1105 d=1; // counts the inner edges 1106 1106 for (i=1;i<=size(pairs)-1;i++) 1107 { 1107 { 1108 1108 for (j=i+1;j<=size(pairs);j++) 1109 1109 { … … 1112 1112 deleted=0; 1113 1113 for (l=size(pairs[j]);l>=1 and deleted==0;l--) 1114 { 1114 { 1115 1115 if (((pairs[i][k][1]==pairs[j][l][1]) and (pairs[i][k][2]==pairs[j][l][2])) or ((pairs[i][k][1]==pairs[j][l][2]) and (pairs[i][k][2]==pairs[j][l][1]))) 1116 1116 { … … 1118 1118 d++; 1119 1119 ggt=abs(gcd(pairs[i][k][1][1]-pairs[i][k][2][1],pairs[i][k][1][2]-pairs[i][k][2][2])); 1120 zwsp=j,ggt; // and it is recorded that the ith and jth 1120 zwsp=j,ggt; // and it is recorded that the ith and jth 1121 1121 // vertex should be connected with mult ggt 1122 1122 vtp[i][6]=intmatconcat(vtp[i][6],zwsp); 1123 1123 zwsp=i,ggt; 1124 1124 vtp[j][6]=intmatconcat(vtp[j][6],zwsp); 1125 pairs[i]=delete(pairs[i],k); // finally the pair is deleted 1125 pairs[i]=delete(pairs[i],k); // finally the pair is deleted 1126 1126 // from both sets of pairs 1127 1127 pairs[j]=delete(pairs[j],l); … … 1163 1163 } 1164 1164 } 1165 // 6.3) Order the vertices such that passing from one to the next we 1165 // 6.3) Order the vertices such that passing from one to the next we 1166 1166 // travel along the boundary of the Newton polytope clock wise. 1167 1167 boundaryNSD=findOrientedBoundary(vertices); 1168 1168 list orderedvertices=boundaryNSD[1]; 1169 1169 // 7) Find the unbounded edges emerging from a vertex in the tropical curve. 1170 // For this we check the remaining pairs for the ith NSD. 1170 // For this we check the remaining pairs for the ith NSD. 1171 1171 // Each pair is ordered according 1172 // to the order in which the vertices occur in orderedvertices. 1172 // to the order in which the vertices occur in orderedvertices. 1173 1173 // The direction of the 1174 // unbounded edge is then the outward pointing primitive normal 1174 // unbounded edge is then the outward pointing primitive normal 1175 1175 // vector to the vector 1176 1176 // pointing from the first vertex in a pair to the second one. … … 1182 1182 for (i=1;i<=size(pairs);i++) 1183 1183 { 1184 list ubedges; // stores the unbounded edges 1184 list ubedges; // stores the unbounded edges 1185 1185 k=1; // counts the unbounded edges 1186 1186 for (j=1;j<=size(pairs[i]);j++) 1187 1187 { 1188 1188 // computes the position of the vertices in the 1189 pos1=positionInList(orderedvertices,pairs[i][j][1]); 1189 pos1=positionInList(orderedvertices,pairs[i][j][1]); 1190 1190 // pair in the list orderedvertices 1191 pos2=positionInList(orderedvertices,pairs[i][j][2]); 1191 pos2=positionInList(orderedvertices,pairs[i][j][2]); 1192 1192 if (((pos1>pos2) and !((pos1==size(orderedvertices)) and (pos2==1))) or ((pos2==size(orderedvertices)) and (pos1==1))) // reorders them if necessary 1193 1193 { … … 1197 1197 } 1198 1198 // the vector pointing from vertex 1 in the pair to vertex2 1199 normalvector=pairs[i][j][2]-pairs[i][j][1]; 1199 normalvector=pairs[i][j][2]-pairs[i][j][1]; 1200 1200 ggt=gcd(normalvector[1],normalvector[2]); // the gcd of the entries 1201 1201 zw=normalvector[2]; // create the outward pointing normal vector … … 1229 1229 kill ubedges; 1230 1230 } 1231 // 8) Store the computed information for the ith part 1231 // 8) Store the computed information for the ith part 1232 1232 // of the NSD in the list graph[i]. 1233 1233 list graph,gr; … … 1235 1235 { 1236 1236 // the first coordinate of the ith vertex of the tropical curve 1237 gr[1]=vtp[i][1]; 1237 gr[1]=vtp[i][1]; 1238 1238 // the second coordinate of the ith vertex of the tropical curve 1239 gr[2]=vtp[i][2]; 1239 gr[2]=vtp[i][2]; 1240 1240 // to which vertices is the ith vertex of the tropical curve connected 1241 gr[3]=vtp[i][6]; 1242 // the directions unbounded edges emerging from the ith 1241 gr[3]=vtp[i][6]; 1242 // the directions unbounded edges emerging from the ith 1243 1243 // vertex of the trop. curve 1244 gr[4]=vtp[i][7]; 1245 // the vertices of the boundary of the ith part of the NSD 1246 gr[5]=vtp[i][3]; 1244 gr[4]=vtp[i][7]; 1245 // the vertices of the boundary of the ith part of the NSD 1246 gr[5]=vtp[i][3]; 1247 1247 graph[i]=gr; 1248 1248 } … … 1263 1263 } 1264 1264 } 1265 // 10) Finally store the boundary vertices and 1265 // 10) Finally store the boundary vertices and 1266 1266 // the inner edges as last entry in the list graph. 1267 1267 // This represents the NSD. … … 1280 1280 // the coordinates of the first vertex are graph[1][1],graph[1][2]; 1281 1281 graph[1][1],graph[1][2]; 1282 // the first vertex is connected to the vertices 1282 // the first vertex is connected to the vertices 1283 1283 // graph[1][3][1,1..ncols(graph[1][3])] 1284 1284 intmat M=graph[1][3]; 1285 1285 M[1,1..ncols(graph[1][3])]; 1286 // the weights of the edges to these vertices are 1286 // the weights of the edges to these vertices are 1287 1287 // graph[1][3][2,1..ncols(graph[1][3])] 1288 1288 M[2,1..ncols(graph[1][3])]; 1289 1289 // from the first vertex emerge size(graph[1][4]) unbounded edges 1290 1290 size(graph[1][4]); 1291 // the primitive integral direction vector of the first unbounded edge 1291 // the primitive integral direction vector of the first unbounded edge 1292 1292 // of the first vertex 1293 1293 graph[1][4][1][1]; 1294 1294 // the weight of the first unbounded edge of the first vertex 1295 1295 graph[1][4][1][2]; 1296 // the monomials which are part of the Newton subdivision of the first vertex 1296 // the monomials which are part of the Newton subdivision of the first vertex 1297 1297 graph[1][5]; 1298 // connecting the points in graph[size(graph)][1] we get 1298 // connecting the points in graph[size(graph)][1] we get 1299 1299 // the boundary of the Newton polytope 1300 1300 graph[size(graph)][1]; 1301 // an entry in graph[size(graph)][2] is a pair of points 1301 // an entry in graph[size(graph)][2] is a pair of points 1302 1302 // in the Newton polytope bounding an inner edge 1303 1303 graph[size(graph)][2][1]; … … 1308 1308 proc drawTropicalCurve (def f,list #) 1309 1309 "USAGE: drawTropicalCurve(f[,#]); f poly or list, # optional list 1310 ASSUME: f is list of linear polynomials of the form ax+by+c with 1311 integers a, b and a rational number c representing a tropical 1310 ASSUME: f is list of linear polynomials of the form ax+by+c with 1311 integers a, b and a rational number c representing a tropical 1312 1312 Laurent polynomial defining a tropical plane curve; 1313 alternatively f can be a polynomial in Q(t)[x,y] defining 1314 a tropical plane curve via the valuation map; 1315 the basering must have a global monomial ordering, two 1313 alternatively f can be a polynomial in Q(t)[x,y] defining 1314 a tropical plane curve via the valuation map; 1315 the basering must have a global monomial ordering, two 1316 1316 variables and up to one parameter! 1317 1317 RETURN: NONE 1318 NOTE: - the procedure creates the files /tmp/tropicalcurveNUMBER.tex and 1319 /tmp/tropicalcurveNUMBER.ps, where NUMBER is a random four 1320 digit integer; 1318 NOTE: - the procedure creates the files /tmp/tropicalcurveNUMBER.tex and 1319 /tmp/tropicalcurveNUMBER.ps, where NUMBER is a random four 1320 digit integer; 1321 1321 moreover it displays the tropical curve via kghostview; 1322 if you wish to remove all these files from /tmp, 1322 if you wish to remove all these files from /tmp, 1323 1323 call the procedure cleanTmp 1324 1324 @* - edges with multiplicity greater than one carry this multiplicity 1325 1325 @* - if # is empty, then the tropical curve is computed w.r.t. minimum, 1326 if #[1] is the string 'max', then it is computed w.r.t. maximum 1327 @* - if the last optional argument is 'onlytexfile' then only the 1328 latex file is produced; this option should be used if kghostview 1326 if #[1] is the string 'max', then it is computed w.r.t. maximum 1327 @* - if the last optional argument is 'onlytexfile' then only the 1328 latex file is produced; this option should be used if kghostview 1329 1329 is not installed on your system 1330 @* - note that lattice points in the Newton subdivision which are 1331 black correspond to markings of the marked subdivision, 1330 @* - note that lattice points in the Newton subdivision which are 1331 black correspond to markings of the marked subdivision, 1332 1332 while lattice points in grey are not marked 1333 1333 EXAMPLE: example drawTropicalCurve shows an example" … … 1347 1347 if (typeof(f)=="poly") 1348 1348 { 1349 // exclude the case that the basering has not precisely 1349 // exclude the case that the basering has not precisely 1350 1350 // one parameter and two indeterminates 1351 1351 if ((npars(basering)!=1) or (nvars(basering)!=2)) 1352 1352 { 1353 ERROR("The basering should have precisely one parameter and two indeterminates!"); 1353 ERROR("The basering should have precisely one parameter and two indeterminates!"); 1354 1354 } 1355 1355 texf=texPolynomial(f); // write the polynomial over Q(t) … … 1362 1362 texf="\\mbox{\\tt The defining equation was not handed over!}"; 1363 1363 list graph=f; 1364 } 1364 } 1365 1365 else 1366 1366 { // write the tropical polynomial defined by f 1367 1367 if (size(#)==0) 1368 { 1368 { 1369 1369 texf="\\min\\{"; 1370 1370 } … … 1375 1375 for (j=1;j<=size(f);j++) 1376 1376 { 1377 texf=texf+texPolynomial(f[j]); 1377 texf=texf+texPolynomial(f[j]); 1378 1378 if (j<size(f)) 1379 1379 { … … 1384 1384 texf=texf+"\\}"; 1385 1385 } 1386 } 1386 } 1387 1387 list graph=tropicalCurve(f,#); // the graph of the tropical polynomial f 1388 1388 } … … 1402 1402 \\addtolength{\\topmargin}{-\\headheight} 1403 1403 \\addtolength{\\topmargin}{-\\topskip} 1404 \\setlength{\\textheight}{267mm} 1404 \\setlength{\\textheight}{267mm} 1405 1405 \\addtolength{\\textheight}{\\topskip} 1406 1406 \\addtolength{\\textheight}{-\\footskip} 1407 1407 \\addtolength{\\textheight}{-30pt} 1408 \\setlength{\\oddsidemargin}{-1in} 1408 \\setlength{\\oddsidemargin}{-1in} 1409 1409 \\addtolength{\\oddsidemargin}{20mm} 1410 1410 \\setlength{\\evensidemargin}{\\oddsidemargin} 1411 \\setlength{\\textwidth}{170mm} 1411 \\setlength{\\textwidth}{170mm} 1412 1412 1413 1413 \\begin{document} 1414 1414 \\parindent0cm 1415 1415 \\begin{center} 1416 \\large\\bf The Tropicalisation of 1416 \\large\\bf The Tropicalisation of 1417 1417 1418 1418 \\bigskip … … 1438 1438 1439 1439 \\begin{center} 1440 "+texDrawNewtonSubdivision(graph,#)+" 1440 "+texDrawNewtonSubdivision(graph,#)+" 1441 1441 \\end{center} 1442 1442 \\end{document}"; … … 1445 1445 int rdnum=random(1000,9999); 1446 1446 write(":w /tmp/tropicalcurve"+string(rdnum)+".tex",TEXBILD); 1447 system("sh","cd /tmp; latex /tmp/tropicalcurve"+string(rdnum)+".tex; dvips /tmp/tropicalcurve"+string(rdnum)+".dvi -o; /bin/rm tropicalcurve"+string(rdnum)+".log; /bin/rm tropicalcurve"+string(rdnum)+".aux; /bin/rm tropicalcurve"+string(rdnum)+".ps?; /bin/rm tropicalcurve"+string(rdnum)+".dvi; kghostview tropicalcurve"+string(rdnum)+".ps &"); 1447 system("sh","cd /tmp; latex /tmp/tropicalcurve"+string(rdnum)+".tex; dvips /tmp/tropicalcurve"+string(rdnum)+".dvi -o; /bin/rm tropicalcurve"+string(rdnum)+".log; /bin/rm tropicalcurve"+string(rdnum)+".aux; /bin/rm tropicalcurve"+string(rdnum)+".ps?; /bin/rm tropicalcurve"+string(rdnum)+".dvi; kghostview tropicalcurve"+string(rdnum)+".ps &"); 1448 1448 } 1449 1449 else … … 1472 1472 proc drawNewtonSubdivision (def f,list #) 1473 1473 "USAGE: drawTropicalCurve(f[,#]); f poly, # optional list 1474 ASSUME: f is list of linear polynomials of the form ax+by+c with integers 1475 a, b and a rational number c representing a tropical Laurent 1474 ASSUME: f is list of linear polynomials of the form ax+by+c with integers 1475 a, b and a rational number c representing a tropical Laurent 1476 1476 polynomial defining a tropical plane curve; 1477 alternatively f can be a polynomial in Q(t)[x,y] defining a tropical 1478 plane curve via the valuation map; 1479 the basering must have a global monomial ordering, two variables 1477 alternatively f can be a polynomial in Q(t)[x,y] defining a tropical 1478 plane curve via the valuation map; 1479 the basering must have a global monomial ordering, two variables 1480 1480 and up to one parameter! 1481 1481 RETURN: NONE 1482 NOTE: - the procedure creates the files /tmp/newtonsubdivisionNUMBER.tex, 1483 and /tmp/newtonsubdivisionNUMBER.ps, where NUMBER is a random 1484 four digit integer; 1482 NOTE: - the procedure creates the files /tmp/newtonsubdivisionNUMBER.tex, 1483 and /tmp/newtonsubdivisionNUMBER.ps, where NUMBER is a random 1484 four digit integer; 1485 1485 moreover it desplays the tropical curve defined by f via kghostview; 1486 if you wish to remove all these files from /tmp, call the procedure 1486 if you wish to remove all these files from /tmp, call the procedure 1487 1487 cleanTmp; 1488 @* if # is empty, then the tropical curve is computed w.r.t. minimum, 1489 if #[1] is the string 'max', then it is computed w.r.t. maximum 1490 @* - note that lattice points in the Newton subdivision which are black 1491 correspond to markings of the marked subdivision, while lattice 1488 @* if # is empty, then the tropical curve is computed w.r.t. minimum, 1489 if #[1] is the string 'max', then it is computed w.r.t. maximum 1490 @* - note that lattice points in the Newton subdivision which are black 1491 correspond to markings of the marked subdivision, while lattice 1492 1492 points in grey are not marked 1493 1493 EXAMPLE: example drawNewtonSubdivision; shows an example" … … 1503 1503 { // write the tropical polynomial defined by f 1504 1504 if (size(#)==0) 1505 { 1505 { 1506 1506 texf="\\min\\{"; 1507 1507 } … … 1512 1512 for (j=1;j<=size(f);j++) 1513 1513 { 1514 texf=texf+texPolynomial(f[j]); 1514 texf=texf+texPolynomial(f[j]); 1515 1515 if (j<size(f)) 1516 1516 { … … 1521 1521 texf=texf+"\\}"; 1522 1522 } 1523 } 1523 } 1524 1524 list graph=tropicalCurve(f,#); // the graph of the tropical polynomial f 1525 1525 } … … 1529 1529 \\parindent0cm 1530 1530 \\begin{center} 1531 \\large\\bf The Newtonsubdivison of 1531 \\large\\bf The Newtonsubdivison of 1532 1532 \\begin{displaymath} 1533 1533 f="+texf+" … … 1537 1537 1538 1538 \\begin{center} 1539 "+texDrawNewtonSubdivision(graph)+ 1539 "+texDrawNewtonSubdivision(graph)+ 1540 1540 " \\end{center} 1541 1541 … … 1543 1543 int rdnum=random(1000,9999); 1544 1544 write(":w /tmp/newtonsubdivision"+string(rdnum)+".tex",TEXBILD); 1545 system("sh","cd /tmp; latex /tmp/newtonsubdivision"+string(rdnum)+".tex; dvips /tmp/newtonsubdivision"+string(rdnum)+".dvi -o; /bin/rm newtonsubdivision"+string(rdnum)+".log; /bin/rm newtonsubdivision"+string(rdnum)+".aux; /bin/rm newtonsubdivision"+string(rdnum)+".ps?; /bin/rm newtonsubdivision"+string(rdnum)+".dvi; kghostview newtonsubdivision"+string(rdnum)+".ps &"); 1545 system("sh","cd /tmp; latex /tmp/newtonsubdivision"+string(rdnum)+".tex; dvips /tmp/newtonsubdivision"+string(rdnum)+".dvi -o; /bin/rm newtonsubdivision"+string(rdnum)+".log; /bin/rm newtonsubdivision"+string(rdnum)+".aux; /bin/rm newtonsubdivision"+string(rdnum)+".ps?; /bin/rm newtonsubdivision"+string(rdnum)+".dvi; kghostview newtonsubdivision"+string(rdnum)+".ps &"); 1546 1546 // return(TEXBILD); 1547 1547 } … … 1568 1568 proc tropicalJInvariant (def f,list #) 1569 1569 "USAGE: tropicalJInvariant(f[,#]); f poly or list, # optional list 1570 ASSUME: f is list of linear polynomials of the form ax+by+c with integers 1571 a, b and a rational number c representing a tropical Laurent 1570 ASSUME: f is list of linear polynomials of the form ax+by+c with integers 1571 a, b and a rational number c representing a tropical Laurent 1572 1572 polynomial defining a tropical plane curve; 1573 alternatively f can be a polynomial in Q(t)[x,y] defining a 1574 tropical plane curve via the valuation map; 1575 @* the basering must have a global monomial ordering, two variables 1573 alternatively f can be a polynomial in Q(t)[x,y] defining a 1574 tropical plane curve via the valuation map; 1575 @* the basering must have a global monomial ordering, two variables 1576 1576 and up to one parameter! 1577 RETURN: number, if the graph underlying the tropical curve has precisely 1578 one loop then its weighted lattice length is returned, 1577 RETURN: number, if the graph underlying the tropical curve has precisely 1578 one loop then its weighted lattice length is returned, 1579 1579 otherwise the result will be -1 1580 NOTE: - if the tropical curve is elliptic and its embedded graph has 1581 precisely one loop, then the weigthed lattice length of 1580 NOTE: - if the tropical curve is elliptic and its embedded graph has 1581 precisely one loop, then the weigthed lattice length of 1582 1582 the loop is its tropical j-invariant 1583 @* - the procedure checks if the embedded graph of the tropical 1584 curve has genus one, but it does NOT check if the loop can 1585 be resolved, so that the curve is not a proper tropical 1586 elliptic curve 1587 @* - if the embedded graph of a tropical elliptic curve has more 1588 than one loop, then all but one can be resolved, but this is 1589 not observed by this procedure, so it will not compute 1583 @* - the procedure checks if the embedded graph of the tropical 1584 curve has genus one, but it does NOT check if the loop can 1585 be resolved, so that the curve is not a proper tropical 1586 elliptic curve 1587 @* - if the embedded graph of a tropical elliptic curve has more 1588 than one loop, then all but one can be resolved, but this is 1589 not observed by this procedure, so it will not compute 1590 1590 the j-invariant 1591 1591 @* - if # is empty, then the tropical curve is computed w.r.t. minimum, 1592 if #[1] is the string 'max', then it is computed w.r.t. maximum 1593 @* - the tropicalJInvariant of a plane tropical cubic is the 1594 'cycle length' of the cubic as introduced in the paper: 1595 Eric Katz, Hannah Markwig, Thomas Markwig: The j-invariant 1592 if #[1] is the string 'max', then it is computed w.r.t. maximum 1593 @* - the tropicalJInvariant of a plane tropical cubic is the 1594 'cycle length' of the cubic as introduced in the paper: 1595 Eric Katz, Hannah Markwig, Thomas Markwig: The j-invariant 1596 1596 of a cubic tropical plane curve. 1597 1597 EXAMPLE: example tropicalJInvariant; shows an example" … … 1607 1607 { 1608 1608 if (typeof(f[1])=="list") 1609 { 1609 { 1610 1610 list graph=f; 1611 1611 } … … 1619 1619 { 1620 1620 ERROR("This is no valid input."); 1621 } 1621 } 1622 1622 } 1623 1623 } … … 1636 1636 genus=-genus/2; // we have counted each bounded edge twice 1637 1637 genus=genus+size(graph); // the genus is 1-#bounded_edges+#vertices 1638 // 3) if the embedded graph has not genus one, 1638 // 3) if the embedded graph has not genus one, 1639 1639 // we cannot compute the j-invariant 1640 1640 if(genus!=1) … … 1648 1648 else 1649 1649 { 1650 intmat nullmat[2][1]; // used to set 1651 // 4) find a vertex which has only one bounded edge, 1652 // if none exists zero is returned, 1650 intmat nullmat[2][1]; // used to set 1651 // 4) find a vertex which has only one bounded edge, 1652 // if none exists zero is returned, 1653 1653 // otherwise the number of the vertex in the list graph 1654 int nonloopvertex=findNonLoopVertex(graph); 1654 int nonloopvertex=findNonLoopVertex(graph); 1655 1655 int dv; //checks if vert. has been found to which nonloopvertex is connected 1656 intmat delvert; // takes for a moment graph[i][3] of the vertex 1656 intmat delvert; // takes for a moment graph[i][3] of the vertex 1657 1657 // to which nonloopvertex is connected 1658 // 5) delete successively vertices in the graph which 1658 // 5) delete successively vertices in the graph which 1659 1659 // have only one bounded edge 1660 1660 while (nonloopvertex>0) 1661 1661 { 1662 // find the only vertex to which the nonloopvertex 1662 // find the only vertex to which the nonloopvertex 1663 1663 // is connected, when it is found 1664 1664 // delete the connection in graph[i][3] and set dv=1 … … 1673 1673 { 1674 1674 delvert=graph[i][3]; 1675 delvert=intmatcoldelete(delvert,j); // delete the connection (note 1675 delvert=intmatcoldelete(delvert,j); // delete the connection (note 1676 1676 // there must have been two!) 1677 1677 dv=1; … … 1681 1681 } 1682 1682 } 1683 graph[nonloopvertex][3]=nullmat; // the only connection of nonloopvertex 1683 graph[nonloopvertex][3]=nullmat; // the only connection of nonloopvertex 1684 1684 // is killed 1685 nonloopvertex=findNonLoopVertex(graph); // find the next vertex 1685 nonloopvertex=findNonLoopVertex(graph); // find the next vertex 1686 1686 // which has only one edge 1687 1687 } … … 1689 1689 intvec loop,weights; // encodes the loop and the edges 1690 1690 i=1; 1691 // start by finding some vertex which belongs to the loop 1692 while (loop==0) 1691 // start by finding some vertex which belongs to the loop 1692 while (loop==0) 1693 1693 { 1694 1694 // if graph[i][3] of a vertex in the loop has 2 columns, all others have 1 1695 if (ncols(graph[i][3])==1) 1695 if (ncols(graph[i][3])==1) 1696 1696 { 1697 1697 i++; … … 1699 1699 else 1700 1700 { 1701 loop[1]=i; // a starting vertex is found 1702 loop[2]=graph[i][3][1,1]; // it is connected to vertex with this number 1701 loop[1]=i; // a starting vertex is found 1702 loop[2]=graph[i][3][1,1]; // it is connected to vertex with this number 1703 1703 weights[2]=graph[i][3][2,1]; // and the edge has this weight 1704 1704 } … … 1708 1708 while (j!=i) // the loop ends with the same vertex with which it starts 1709 1709 { 1710 // the first row of graph[j][3] has two entries 1710 // the first row of graph[j][3] has two entries 1711 1711 // corresponding to the two vertices 1712 // to which the active vertex j is connected; 1712 // to which the active vertex j is connected; 1713 1713 // one is loop[k-1], i.e. the one which 1714 1714 // precedes j in the loop; we have to choose the other one 1715 if (graph[j][3][1,1]==loop[k-1]) 1715 if (graph[j][3][1,1]==loop[k-1]) 1716 1716 { 1717 1717 loop[k+1]=graph[j][3][1,2]; … … 1724 1724 } 1725 1725 j=loop[k+1]; // set loop[k+1] the new active vertex 1726 k++; 1727 } 1726 k++; 1727 } 1728 1728 // 7) compute for each edge in the loop the lattice length 1729 poly xcomp,ycomp; // the x- and y-components of the vectors 1729 poly xcomp,ycomp; // the x- and y-components of the vectors 1730 1730 // connecting two vertices of the loop 1731 number nenner; // the product of the denominators of 1731 number nenner; // the product of the denominators of 1732 1732 // the x- and y-components 1733 1733 number jinvariant; // the j-invariant 1734 int eins,zwei,ggt; 1734 int eins,zwei,ggt; 1735 1735 for (i=1;i<=size(loop)-1;i++) // compute the lattice length for each edge 1736 1736 { 1737 xcomp=graph[loop[i]][1]-graph[loop[i+1]][1]; 1738 ycomp=graph[loop[i]][2]-graph[loop[i+1]][2]; 1737 xcomp=graph[loop[i]][1]-graph[loop[i+1]][1]; 1738 ycomp=graph[loop[i]][2]-graph[loop[i+1]][2]; 1739 1739 nenner=denominator(leadcoef(xcomp))*denominator(leadcoef(ycomp)); 1740 execute("eins="+string(numerator(leadcoef(nenner*xcomp)))+";"); 1741 execute("zwei="+string(numerator(leadcoef(nenner*ycomp)))+";"); 1740 execute("eins="+string(numerator(leadcoef(nenner*xcomp)))+";"); 1741 execute("zwei="+string(numerator(leadcoef(nenner*ycomp)))+";"); 1742 1742 ggt=gcd(eins,zwei); // the lattice length is the "gcd" 1743 1743 // of the x-component and the y-component 1744 jinvariant=jinvariant+ggt/(nenner*weights[i+1]); // divided by the 1744 jinvariant=jinvariant+ggt/(nenner*weights[i+1]); // divided by the 1745 1745 // weight of the edge 1746 1746 } 1747 return(jinvariant); 1747 return(jinvariant); 1748 1748 } 1749 1749 } … … 1759 1759 // the curve can have arbitrary degree 1760 1760 tropicalJInvariant(t*(x7+y7+1)+1/t*(x4+y4+x2+y2+x3y+xy3)+1/t7*x2y2); 1761 // the procedure does not realise, if the embedded graph of the tropical 1761 // the procedure does not realise, if the embedded graph of the tropical 1762 1762 // curve has a loop that can be resolved 1763 1763 tropicalJInvariant(1+x+y+xy+tx2y+txy2); 1764 1764 // but it does realise, if the curve has no loop at all ... 1765 1765 tropicalJInvariant(x+y+1); 1766 // or if the embedded graph has more than one loop - even if only one 1766 // or if the embedded graph has more than one loop - even if only one 1767 1767 // cannot be resolved 1768 1768 tropicalJInvariant(1+x+y+xy+tx2y+txy2+t3x5+t3y5+tx2y2+t2xy4+t2yx4); … … 1773 1773 proc weierstrassForm (poly f,list #) 1774 1774 "USAGE: weierstrassForm(wf[,#]); wf poly, # list 1775 ASSUME: wf is a a polynomial whose Newton polygon has precisely one 1776 interior lattice point, so that it defines an elliptic curve 1775 ASSUME: wf is a a polynomial whose Newton polygon has precisely one 1776 interior lattice point, so that it defines an elliptic curve 1777 1777 on the toric surface corresponding to the Newton polygon 1778 1778 RETURN: poly, the Weierstrass normal form of the polynomial … … 1780 1780 to Fernando Rodriguez Villegas, villegas@math.utexas.edu 1781 1781 @* - the characteristic of the base field should not be 2 or 3 1782 @* - if an additional argument # is given, a simplified Weierstrass 1782 @* - if an additional argument # is given, a simplified Weierstrass 1783 1783 form is computed 1784 1784 EXAMPLE: example weierstrassForm; shows an example" … … 1841 1841 1842 1842 proc jInvariant (poly f,list #) 1843 "USAGE: jInvariant(f[,#]); f poly, # list 1844 ASSUME: - f is a a polynomial whose Newton polygon has precisely one 1845 interior lattice point, so that it defines an elliptic curve 1843 "USAGE: jInvariant(f[,#]); f poly, # list 1844 ASSUME: - f is a a polynomial whose Newton polygon has precisely one 1845 interior lattice point, so that it defines an elliptic curve 1846 1846 on the toric surface corresponding to the Newton polygon 1847 @* - it the optional argument # is present the base field should be 1848 Q(t) and the optional argument should be one of the following 1847 @* - it the optional argument # is present the base field should be 1848 Q(t) and the optional argument should be one of the following 1849 1849 strings: 1850 @* 'ord' : then the return value is of type integer, 1850 @* 'ord' : then the return value is of type integer, 1851 1851 namely the order of the j-invariant 1852 @* 'split' : then the return value is a list of two polynomials, 1852 @* 'split' : then the return value is a list of two polynomials, 1853 1853 such that the quotient of these two is the j-invariant 1854 1854 RETURN: poly, the j-invariant of the elliptic curve defined by poly 1855 NOTE: the characteristic of the base field should not be 2 or 3, 1855 NOTE: the characteristic of the base field should not be 2 or 3, 1856 1856 unless the input is a plane cubic 1857 1857 EXAMPLE: example jInvariant; shows an example" … … 1879 1879 echo=2; 1880 1880 ring r=(0,t),(x,y),dp; 1881 // jInvariant computes the j-invariant of a cubic 1881 // jInvariant computes the j-invariant of a cubic 1882 1882 jInvariant(x+y+x2y+y3+1/t*xy); 1883 // if the ground field has one parameter t, then we can instead 1883 // if the ground field has one parameter t, then we can instead 1884 1884 // compute the order of the j-invariant 1885 1885 jInvariant(x+y+x2y+y3+1/t*xy,"ord"); … … 1889 1889 poly h=x22y11+x19y10+x17y9+x16y9+x12y7+x9y6+x7y5+x2y3+x14y8; 1890 1890 // its j-invariant is 1891 jInvariant(h); 1891 jInvariant(h); 1892 1892 } 1893 1893 … … 1908 1908 @* l[6] = a list containing the vertices of the tropical conic f 1909 1909 @* l[7] = a list containing lists with vertices of the tangents 1910 @* l[8] = a string which contains the latex-code to draw the 1910 @* l[8] = a string which contains the latex-code to draw the 1911 1911 tropical conic and its tropicalised tangents 1912 1912 @* l[9] = if # is non-empty, this is the same data for the dual … … 1932 1932 ring LINRING=(0,t),(x,y,a(1..6)),lp; 1933 1933 list points=imap(BASERING,points); 1934 ideal I; // the ideal will contain the linear equations given by the conic 1934 ideal I; // the ideal will contain the linear equations given by the conic 1935 1935 // and the points 1936 1936 for (i=1;i<=5;i++) … … 1953 1953 ring tRING=0,t,ls; 1954 1954 list pointdenom=imap(BASERING,pointdenom); 1955 list pointnum=imap(BASERING,pointnum); 1955 list pointnum=imap(BASERING,pointnum); 1956 1956 intvec pointcoordinates; 1957 1957 for (i=1;i<=size(pointdenom);i++) … … 2029 2029 We consider the concic through the following five points: 2030 2030 \\begin{displaymath} 2031 "; 2031 "; 2032 2032 string texf=texDrawTropical(graphf,list("",scalefactor)); 2033 2033 for (i=1;i<=size(points);i++) … … 2067 2067 \\end{document}"; 2068 2068 setring BASERING; 2069 // If # non-empty, compute the dual conic and the tangents 2070 // through the dual points 2069 // If # non-empty, compute the dual conic and the tangents 2070 // through the dual points 2071 2071 // corresponding to the tangents of the given conic. 2072 2072 if (size(#)>0) 2073 2073 { 2074 2074 list dualpoints; 2075 for (i=1;i<=size(points);i++) 2075 for (i=1;i<=size(points);i++) 2076 2076 { 2077 2077 dualpoints[i]=list(leadcoef(tangents[i])/substitute(tangents[i],x,0,y,0),leadcoef(tangents[i]-lead(tangents[i]))/substitute(tangents[i],x,0,y,0)); … … 2097 2097 // conic[2] is the equation of the conic f passing through the five points 2098 2098 conic[2]; 2099 // conic[3] is a list containing the equations of the tangents 2099 // conic[3] is a list containing the equations of the tangents 2100 2100 // through the five points 2101 2101 conic[3]; 2102 2102 // conic[4] is an ideal representing the tropicalisation of the conic f 2103 2103 conic[4]; 2104 // conic[5] is a list containing the tropicalisation 2104 // conic[5] is a list containing the tropicalisation 2105 2105 // of the five tangents in conic[3] 2106 2106 conic[5]; 2107 // conic[6] is a list containing the vertices of the tropical conic 2107 // conic[6] is a list containing the vertices of the tropical conic 2108 2108 conic[6]; 2109 2109 // conic[7] is a list containing the vertices of the five tangents 2110 2110 conic[7]; 2111 // conic[8] contains the latex code to draw the tropical conic and 2112 // its tropicalised tangents; it can written in a file, processed and 2113 // displayed via kghostview 2114 write(":w /tmp/conic.tex",conic[8]); 2115 system("sh","cd /tmp; latex /tmp/conic.tex; dvips /tmp/conic.dvi -o; 2111 // conic[8] contains the latex code to draw the tropical conic and 2112 // its tropicalised tangents; it can written in a file, processed and 2113 // displayed via kghostview 2114 write(":w /tmp/conic.tex",conic[8]); 2115 system("sh","cd /tmp; latex /tmp/conic.tex; dvips /tmp/conic.dvi -o; 2116 2116 kghostview conic.ps &"); 2117 // with an optional argument the same information for the dual conic is computed 2117 // with an optional argument the same information for the dual conic is computed 2118 2118 // and saved in conic[9] 2119 2119 conic=conicWithTangents(points,1); 2120 2120 conic[9][2]; // the equation of the dual conic 2121 2121 } 2122 2122 2123 2123 /////////////////////////////////////////////////////////////////////////////// 2124 2124 /// Procedures concerned with tropicalisation … … 2129 2129 ASSUME: f is a polynomial in Q(t)[x_1,...,x_n] 2130 2130 RETURN: list, the linear forms of the tropicalisation of f 2131 NOTE: if # is empty, then the valuation of t will be 1, 2132 @* if # is the string 'max' it will be -1; 2133 @* the latter supposes that we consider the maximum of the 2131 NOTE: if # is empty, then the valuation of t will be 1, 2132 @* if # is the string 'max' it will be -1; 2133 @* the latter supposes that we consider the maximum of the 2134 2134 computed linear forms, the former that we consider their minimum 2135 2135 EXAMPLE: example tropicalise; shows an example" … … 2154 2154 { 2155 2155 tropicalf[i]=tropicalf[i]+exp[j]*var(j); 2156 } 2156 } 2157 2157 f=f-lead(f); 2158 2158 } … … 2194 2194 ///////////////////////////////////////////////////////////////////////// 2195 2195 2196 proc tInitialForm (poly f, intvec w) 2196 proc tInitialForm (poly f, intvec w) 2197 2197 "USAGE: tInitialForm(f,w); f a polynomial, w an integer vector 2198 2198 ASSUME: f is a polynomial in Q[t,x_1,...,x_n] and w=(w_0,w_1,...,w_n) 2199 2199 RETURN: poly, the t-initialform of f(t,x) w.r.t. w evaluated at t=1 2200 NOTE: the t-initialform is the sum of the terms with MAXIMAL 2200 NOTE: the t-initialform is the sum of the terms with MAXIMAL 2201 2201 weighted order w.r.t. w 2202 2202 EXAMPLE: example tInitialForm; shows an example" … … 2208 2208 // do the same for the remaining part of f and compare the results 2209 2209 // keep only the smallest ones 2210 int vglgewicht; 2211 f=f-lead(f); 2210 int vglgewicht; 2211 f=f-lead(f); 2212 2212 while (f!=0) 2213 2213 { … … 2224 2224 initialf=initialf+lead(f); 2225 2225 } 2226 } 2226 } 2227 2227 f=f-lead(f); 2228 2228 } … … 2244 2244 proc tInitialIdeal (ideal i,intvec w,list #) 2245 2245 "USAGE: tInitialIdeal(i,w); i ideal, w intvec 2246 ASSUME: i is an ideal in Q[t,x_1,...,x_n] and w=(w_0,...,w_n) 2246 ASSUME: i is an ideal in Q[t,x_1,...,x_n] and w=(w_0,...,w_n) 2247 2247 RETURN: ideal ini, the t-initial ideal of i with respect to w" 2248 2248 { 2249 2249 // THE PROCEDURE WILL BE CALLED FROM OTHER PROCEDURES INSIDE THIS LIBRARY; 2250 // IN THIS CASE THE VARIABLE t WILL INDEED BE THE LAST VARIABLE INSTEAD OF 2250 // IN THIS CASE THE VARIABLE t WILL INDEED BE THE LAST VARIABLE INSTEAD OF 2251 2251 // THE FIRST, 2252 2252 // AND WE THEREFORE HAVE TO MOVE IT BACK TO THE FRONT! … … 2272 2272 // ... and compute a standard basis with 2273 2273 // respect to the homogenised ordering defined by w. Since the generators 2274 // of i will be homogeneous it we can instead take the ordering wp 2275 // with the weightvector (0,w) replaced by (M,...,M)+(0,w) for some 2276 // large M, so that all entries are positive 2277 int M=-minInIntvec(w)[1]+1; // find M such that w[j]+M is 2274 // of i will be homogeneous it we can instead take the ordering wp 2275 // with the weightvector (0,w) replaced by (M,...,M)+(0,w) for some 2276 // large M, so that all entries are positive 2277 int M=-minInIntvec(w)[1]+1; // find M such that w[j]+M is 2278 2278 // strictly positive for all j 2279 2279 intvec whomog=M+1; … … 2283 2283 } 2284 2284 execute("ring WEIGHTRING=("+charstr(basering)+"),("+varstr(basering)+"),(wp("+string(whomog)+"));"); 2285 // map i to the new ring and compute a GB of i, then dehomogenise i, 2286 // so that we can be sure, that the 2285 // map i to the new ring and compute a GB of i, then dehomogenise i, 2286 // so that we can be sure, that the 2287 2287 // initial forms of the generators generate the initial ideal 2288 2288 ideal i=subst(groebner(imap(HOMOGRING,i)),@s,1); … … 2538 2538 proc texDrawBasic (list texdraw) 2539 2539 "USAGE: texDrawBasic(texdraw); list texdraw 2540 ASSUME: texdraw is a list of strings representing texdraw commands 2541 (as produced by texDrawTropical) which should be embedded into 2540 ASSUME: texdraw is a list of strings representing texdraw commands 2541 (as produced by texDrawTropical) which should be embedded into 2542 2542 a texdraw environment 2543 2543 RETURN: string, a texdraw environment enclosing the input … … 2559 2559 \\end{texdraw}"; 2560 2560 return(texdrawtp); 2561 } 2561 } 2562 2562 example 2563 2563 { … … 2576 2576 ASSUME: graph is the output of tropicalCurve 2577 2577 RETURN: string, the texdraw code of the tropical plane curve encoded by graph 2578 NOTE: - if the list # is non-empty, the first entry should be a string; 2579 if this string is 'max', then the tropical curve is considered 2580 with respect to the maximum; otherwise the curve is considered 2581 with respect to the minimum and the string can be used to insert 2578 NOTE: - if the list # is non-empty, the first entry should be a string; 2579 if this string is 'max', then the tropical curve is considered 2580 with respect to the maximum; otherwise the curve is considered 2581 with respect to the minimum and the string can be used to insert 2582 2582 further texdraw commands (e.g. to have a lighter image as when called 2583 from inside conicWithTangents); 2584 @* - the procedure computes a scalefactor for the texdraw command which 2585 should help to display the curve in the right way; this may, 2586 however, be a bad idea if several texDrawTropical outputs are 2587 put together to form one image; the scalefactor can be prescribed 2583 from inside conicWithTangents); 2584 @* - the procedure computes a scalefactor for the texdraw command which 2585 should help to display the curve in the right way; this may, 2586 however, be a bad idea if several texDrawTropical outputs are 2587 put together to form one image; the scalefactor can be prescribed 2588 2588 by a second optional entry of type poly 2589 2589 @* - the list # is optional and may as well be empty … … 2591 2591 { 2592 2592 int i,j; 2593 // deal first with the pathological case that 2593 // deal first with the pathological case that 2594 2594 // the input polynomial was a monomial 2595 // and does therefore not define a tropical curve, 2596 // and check if the Newton polytope is 2595 // and does therefore not define a tropical curve, 2596 // and check if the Newton polytope is 2597 2597 // a line segment so that the curve defines a bunch of lines 2598 2598 int bunchoflines; 2599 2599 // if the boundary of the Newton polytope consists of a single point 2600 if (size(graph[size(graph)][1])==1) 2600 if (size(graph[size(graph)][1])==1) 2601 2601 { 2602 2602 return(string()); … … 2611 2611 } 2612 2612 // then the Newton polytope is a line segment 2613 if ((size(graph[size(graph)][1])-size(syz(M)))==1) 2613 if ((size(graph[size(graph)][1])-size(syz(M)))==1) 2614 2614 { 2615 2615 bunchoflines=1; … … 2618 2618 // go on with the case that a tropical curve is defined 2619 2619 if (size(#)==0) 2620 { 2620 { 2621 2621 string texdrawtp=" 2622 2622 … … 2626 2626 { 2627 2627 if ((#[1]!="max") and (#[1]!="")) 2628 { 2628 { 2629 2629 string texdrawtp=#[1]; 2630 2630 } … … 2675 2675 { 2676 2676 // if the curve is a bunch of lines no vertex has to be drawn 2677 if (bunchoflines==0) 2677 if (bunchoflines==0) 2678 2678 { 2679 2679 texdrawtp=texdrawtp+" … … 2681 2681 } 2682 2682 // draw the bounded edges emerging from the ith vertex 2683 for (j=1;j<=ncols(graph[i][3]);j++) 2684 { 2685 // don't draw it twice - and if there is only one vertex 2683 for (j=1;j<=ncols(graph[i][3]);j++) 2684 { 2685 // don't draw it twice - and if there is only one vertex 2686 2686 // and graph[i][3][1,1] is thus 0, nothing is done 2687 if (i<graph[i][3][1,j]) 2688 { 2687 if (i<graph[i][3][1,j]) 2688 { 2689 2689 texdrawtp=texdrawtp+" 2690 2690 \\move ("+decimal(graph[i][1]-centerx)+" "+decimal(graph[i][2]-centery)+") \\lvec ("+decimal(graph[graph[i][3][1,j]][1]-centerx)+" "+decimal(graph[graph[i][3][1,j]][2]-centery)+")"; 2691 2691 // if the multiplicity is more than one, denote it in the picture 2692 if (graph[i][3][2,j]>1) 2692 if (graph[i][3][2,j]>1) 2693 2693 { 2694 2694 texdrawtp=texdrawtp+" … … 2699 2699 // draw the unbounded edges emerging from the ith vertex 2700 2700 // they should not be too long 2701 for (j=1;j<=size(graph[i][4]);j++) 2702 { 2701 for (j=1;j<=size(graph[i][4]);j++) 2702 { 2703 2703 relxy=shorten(list(decimal(3*graph[i][4][j][1][1]/scalefactor),decimal(3*graph[i][4][j][1][2]/scalefactor),"2.5")); 2704 2704 texdrawtp=texdrawtp+" 2705 2705 \\move ("+decimal(graph[i][1]-centerx)+" "+decimal(graph[i][2]-centery)+") \\rlvec ("+relxy[1]+" "+relxy[2]+")"; 2706 2706 // if the multiplicity is more than one, denote it in the picture 2707 if (graph[i][4][j][2]>1) 2707 if (graph[i][4][j][2]>1) 2708 2708 { 2709 2709 texdrawtp=texdrawtp+" … … 2783 2783 poly scalefactor=minOfPolys(list(12/leadcoef(maxx),12/leadcoef(maxy))); 2784 2784 if (scalefactor<1) 2785 { 2785 { 2786 2786 subdivision=subdivision+" 2787 2787 \\relunitscale"+ decimal(scalefactor); … … 2791 2791 { 2792 2792 subdivision=subdivision+" 2793 \\move ("+string(boundary[i][1])+" "+string(boundary[i][2])+") 2793 \\move ("+string(boundary[i][1])+" "+string(boundary[i][2])+") 2794 2794 \\lvec ("+string(boundary[i+1][1])+" "+string(boundary[i+1][2])+")"; 2795 } 2795 } 2796 2796 subdivision=subdivision+" 2797 \\move ("+string(boundary[size(boundary)][1])+" "+string(boundary[size(boundary)][2])+") 2797 \\move ("+string(boundary[size(boundary)][1])+" "+string(boundary[size(boundary)][2])+") 2798 2798 \\lvec ("+string(boundary[1][1])+" "+string(boundary[1][2])+") 2799 2799 … … 2802 2802 { 2803 2803 subdivision=subdivision+" 2804 \\move ("+string(inneredges[i][1][1])+" "+string(inneredges[i][1][2])+") 2804 \\move ("+string(inneredges[i][1][1])+" "+string(inneredges[i][1][2])+") 2805 2805 \\lvec ("+string(inneredges[i][2][1])+" "+string(inneredges[i][2][2])+")"; 2806 2806 } … … 2822 2822 } 2823 2823 } 2824 // deal with the pathological cases 2824 // deal with the pathological cases 2825 2825 if (size(boundary)==1) // then the Newton polytope is a point 2826 2826 { … … 2877 2877 { 2878 2878 subdivision=subdivision+" 2879 \\move ("+string(markings[i][1])+" "+string(markings[i][2])+") 2879 \\move ("+string(markings[i][1])+" "+string(markings[i][2])+") 2880 2880 \\fcir f:0 r:"+decimal(2/(8*scalefactor),size(string(int(scalefactor)))+1); 2881 2881 } 2882 2882 // enclose subdivision in the texdraw environment 2883 string texsubdivision=" 2883 string texsubdivision=" 2884 2884 \\begin{texdraw} 2885 \\drawdim cm \\relunitscale 1 2885 \\drawdim cm \\relunitscale 1 2886 2886 \\linewd 0.05" 2887 2887 +subdivision+" … … 2897 2897 poly f=x+y+x2y+xy2+1/t*xy; 2898 2898 list graph=tropicalCurve(f); 2899 // compute the texdraw code of the Newton subdivision of the tropical curve 2899 // compute the texdraw code of the Newton subdivision of the tropical curve 2900 2900 texDrawNewtonSubdivision(graph); 2901 2901 } … … 2905 2905 proc texDrawTriangulation (list triang,list polygon) 2906 2906 "USAGE: texDrawTriangulation(triang,polygon); triang,polygon list 2907 ASSUME: polygon is a list of integer vectors describing the 2907 ASSUME: polygon is a list of integer vectors describing the 2908 2908 lattice points of a marked polygon; 2909 triang is a list of integer vectors describing a 2909 triang is a list of integer vectors describing a 2910 2910 triangulation of the marked polygon 2911 in the sense that an integer vector of the form (i,j,k) describes 2911 in the sense that an integer vector of the form (i,j,k) describes 2912 2912 the triangle formed by polygon[i], polygon[j] and polygon[k] 2913 RETURN: string, a texdraw code for the triangulation described 2913 RETURN: string, a texdraw code for the triangulation described 2914 2914 by triang without the texdraw environment 2915 2915 EXAMPLE: example texDrawTriangulation; shows an example" … … 2920 2920 "; 2921 2921 int i,j; // indices 2922 list pairs,markings; // stores edges of the triangulation, respecively 2923 // the marked points for each triangle store the edges and marked 2922 list pairs,markings; // stores edges of the triangulation, respecively 2923 // the marked points for each triangle store the edges and marked 2924 2924 // points of the triangle 2925 2925 for (i=1;i<=size(triang);i++) … … 2930 2930 markings[3*i-2]=triang[i][1]; 2931 2931 markings[3*i-1]=triang[i][2]; 2932 markings[3*i]=triang[i][3]; 2932 markings[3*i]=triang[i][3]; 2933 2933 } 2934 2934 // delete redundant pairs which occur more than once … … 2963 2963 { 2964 2964 latex=latex+" 2965 \\move ("+string(polygon[markings[i]][1])+" "+string(polygon[markings[i]][2])+") 2965 \\move ("+string(polygon[markings[i]][1])+" "+string(polygon[markings[i]][2])+") 2966 2966 \\fcir f:0 r:0.08"; 2967 2967 } … … 2970 2970 { 2971 2971 latex=latex+" 2972 \\move ("+string(polygon[pairs[i][1]][1])+" "+string(polygon[pairs[i][1]][2])+") 2972 \\move ("+string(polygon[pairs[i][1]][1])+" "+string(polygon[pairs[i][1]][2])+") 2973 2973 \\lvec ("+string(polygon[pairs[i][2]][1])+" "+string(polygon[pairs[i][2]][2])+")"; 2974 2974 } … … 2977 2977 { 2978 2978 latex=latex+" 2979 \\move ("+string(polygon[i][1])+" "+string(polygon[i][2])+") 2979 \\move ("+string(polygon[i][1])+" "+string(polygon[i][2])+") 2980 2980 \\fcir f:0.7 r:0.04"; 2981 2981 } … … 2986 2986 "EXAMPLE:"; 2987 2987 echo=2; 2988 // the lattice polygon spanned by the points (0,0), (3,0) and (0,3) 2988 // the lattice polygon spanned by the points (0,0), (3,0) and (0,3) 2989 2989 // with all in