// version="$Id$"; category="Tropical Geometry"; info=" LIBRARY: tropical.lib Computations in Tropical Geometry AUTHORS: Anders Jensen Needergard, email: jensen@math.tu-berlin.de @* Hannah Markwig, email: hannah@uni-math.gwdg.de @* Thomas Markwig, email: keilen@mathematik.uni-kl.de WARNING: - tropicalLifting will only work with LINUX and if in addition gfan is installed. @*- drawTropicalCurve and drawTropicalNewtonSubdivision will only display the @* tropical curve with LINUX and if in addition latex and kghostview @* are installed. @*- For tropicalLifting in the definition of the basering the parameter t @* from the Puiseux series field C{{t}} must be defined as a variable, @* while for all other procedures it must be defined as a parameter. THEORY: Fix some base field K and a bunch of lattice points v0,...,vm in the integer lattice Z^n, then this defines a toric variety as the closure of (K*)^n in the projective space P^m, where the torus is embedded via the map sending a point x in (K*)^n to the point (x^v0,...,x^vm). The generic hyperplane sections are just the images of the hypersurfaces in (K*)^n defined by the polynomials f=a0*x^v0+...+am*x^vm=0. Some properties of these hypersurfaces can be studied via tropicalisation. For this we suppose that K=C{{t}} is the field of Puiseux series over the field of complex numbers (or any other field with a valuation into the real numbers). One associates to the hypersurface given by f=a0*x^v0+...+am*x^vm the tropical hypersurface defined by the tropicalisation trop(f)=min{val(a0)+,...,val(am)+}. Here, denotes the standard scalar product of the integer vector v in Z^n with the vector x=(x1,...,xn) of variables, so that trop(f) is a piecewise linear function on R^n. The corner locus of this function (i.e. the points at which the minimum is attained a least twice) is the tropical hypersurface defined by trop(f). The theorem of Newton-Kapranov states that this tropical hypersurface is the same as if one computes pointwise the valuation of the hypersurface given by f. The analogue holds true if one replaces one equation f by an ideal I. A constructive proof of the theorem is given by an adapted version of the Newton-Puiseux algorithm. The hard part is to find a point in the variety over C{{t}} which corresponds to a given point in the tropical variety. It is the purpose of this library to provide basic means to deal with tropical varieties. Of course we cannot represent the field of Puiseux series over C in its full strength, however, in order to compute interesting examples it will be sufficient to replace the complex numbers C by the rational numbers Q and to replace Puiseux series in t by rational functions in t, i.e. we replace C{{t}} by Q(t), or sometimes even by Q[t]. Note, that this in particular forbids rational exponents for the t's. Moreover, in @sc{Singular} no negative exponents of monomials are allowed, so that the integer vectors vi will have to have non-negative entries. Shifting all exponents by a fixed integer vector does not change the tropicalisation nor does it change the toric variety. Thus this does not cause any restriction. If, however, for some reason you prefer to work with general vi, then you have to pass right away to the tropicalisation of the equations, whereever this is allowed -- these are linear polynomials where the constant coefficient corresponds to the valuation of the original coefficient and where the non-constant coefficient correspond to the exponents of the monomials, thus they may be rational numbers respectively negative numbers: 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}. The main tools provided in this library are as follows: @* - tropicalLifting implements the constructive proof of the Theorem of Newton-Kapranov and constructs a point in the variety over C{{t}} corresponding to a given point in the corresponding tropical variety associated to an ideal I; the generators of I have to be in the polynomial ring Q[t,x1,...,xn] considered as a subring of C{{t}}[x1,...,xn]; a solution will be constructed up to given order; note that several field extensions of Q might be necessary throughout the intermediate computations; the procedures use the external program gfan @* - puiseuxExpansion computes a Newton-Puiseux expansion of a plane curve singularity @* - drawTropicalCurve visualises a tropical plane curve either given by a polynomial in Q(t)[x,y] or by a list of linear polynomials of the form ax+by+c with a,b in Z and c in Q; latex must be installed on your computer @* - tropicalJInvariant computes the tropical j-invaiant of a tropical elliptic curve @* - jInvariant computes the j-invariant of an elliptic curve @* - weierstrassForm computes the Weierstrass form of an elliptic curve PROCEDURES FOR TROPICAL LIFTING: tropicalLifting() computes a point in the tropical variety displayTropicalLifting() displays the output of tropicalLifting puiseuxExpansion() computes a Newton-Puiseux expansion in the plane displayPuiseuxExpansion() displays the output of puiseuxExpansion PROCEDURES FOR DRAWING TROPICAL CURVES: tropicalCurve() computes a tropical curve and its Newton subdivision drawTropicalCurve() produces a post script image of a tropical curve drawNewtonSubdivision() produces a post script image of a Newton subdivision PROCEDURES FOR J-INVARIANTS: tropicalJInvariant() computes the tropical j-invariant of a tropical curve weierstrassForm() computes the Weierstrass form of a cubic polynomial jInvariant() computes the j-invariant of a cubic polynomial GENERAL PROCEDURES: conicWithTangents() computes a conic through five points with tangents tropicalise() computes the tropicalisation of a polynomial tropicaliseSet() computes the tropicalisation several polynomials tInitialForm() computes the tInitial form of a polynomial in Q[t,x_1,...,x_n] tInitialIdeal() computes the tInitial ideal of an ideal in Q[t,x_1,...,x_n] initialForm() computes the initial form of poly in Q[x_1,...,x_n] initialIdeal() computes the initial ideal of an ideal in Q[x_1,...,x_n] PROCEDURES FOR LATEX CONVERSION: texNumber() outputs the texcommand for the leading coefficient of poly texPolynomial() outputs the texcommand for the polynomial poly texMatrix() outputs the texcommand for the matrix texDrawBasic() embeds output of texDrawTropical in a texdraw environment texDrawTropical() computes the texdraw commands for a tropical curve texDrawNewtonSubdivision() computes texdraw commands for a Newton subdivision texDrawTriangulation() computes texdraw commands for a triangulation AUXILARY PROCEDURES: radicalMemberShip() checks radical membership tInitialFormPar() computes the t-initial form of poly in Q(t)[x_1,...,x_n] tInitialFormParMax() same as tInitialFormPar, but uses maximum solveTInitialFormPar() displays approximated solution of a 0-dim ideal detropicalise() computes the detropicalisation of a linear form tDetropicalise() computes the detropicalisation of a linear form dualConic() computes the dual of an affine plane conic parameterSubstitute() substitutes in the polynomial the parameter t by t^N tropicalSubst() makes certain substitutions in a tropical polynomial randomPoly() computes a polynomial with random coefficients cleanTmp() clears /tmp from files created by other procedures KEYWORDS: tropical curves; tropical polynomials "; /////////////////////////////////////////////////////////////////////////////// /// Auxilary Static Procedures in this Library /////////////////////////////////////////////////////////////////////////////// /// - phiOmega /// - cutdown /// - tropicalparametriseNoabs /// - findzeros /// - basictransformideal /// - testw /// - simplifyToOrder /// - scalarproduct /// - intvecdelete /// - invertorder /// - ordermaximalidealsNoabs /// - displaypoly /// - displaycoef /// - choosegfanvector /// - tropicalliftingresubstitute /// - tropicalparametrise /// - eliminatecomponents /// - findzerosAndBasictransform /// - ordermaximalideals /// - verticesTropicalCurve /// - bunchOfLines /// - clearintmat /// - sortintvec /// - sortintmat /// - intmatcoldelete /// - intmatconcat /// - minInIntvec /// - positionInList /// - sortlist /// - minInList /// - vergleiche /// - koeffizienten /// - minOfPolys /// - shorten /// - minOfStringDecimal /// - decimal /// - stringcontainment /// - stringdelete /// - stringinsert /// - texmonomial /// - texcoefficient /// - abs /// - findNonLoopVertex /// - coordinatechange /// - weierstrassFormOfACubic /// - weierstrassFormOfA4x2Curve /// - weierstrassFormOfA2x2Curve /// - jInvariantOfACubic /// - jInvariantOfA4x2Curve /// - jInvariantOfA2x2Curve /// - jInvariantOfAPuiseuxCubic ////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////// LIB "random.lib"; LIB "solve.lib"; LIB "poly.lib"; LIB "elim.lib"; LIB "linalg.lib"; LIB "oldpolymake.lib"; LIB "primdec.lib"; LIB "absfact.lib"; LIB "hnoether.lib"; ////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// /// Procedures concerned with tropical parametrisation /////////////////////////////////////////////////////////////////////////////// proc tropicalLifting (ideal i,intvec w,int ordnung,list #) "USAGE: tropicalLifting(i,w,ord[,opt]); i ideal, w intvec, ord int, opt string ASSUME: - i is an ideal in Q[t,x_1,...,x_n], w=(w_0,w_1,...,w_n) and (w_1/w_0,...,w_n/w_0) is in the tropical variety of i, and ord is the order up to which a point in V(i) over Q{{t}} lying over (w_1/w_0,...,w_n/w_0) shall be computed; w_0 may NOT be ZERO @* - the basering should not have any parameters on its own and it should have a global monomial ordering, e.g. ring r=0,(t,x(1..n)),dp; @* - the first variable of the basering will be treated as the parameter t in the Puiseux series field @* - the optional parameter opt should be one or more strings among the following: @* 'isZeroDimensional' : the dimension i is zero (not to be checked); @* 'isPrime' : the ideal is prime over Q(t)[x_1,...,x_n] (not to be checked); @* 'isInTrop' : (w_1/w_0,...,w_n/w_0) is in the tropical variety (not to be checked); @* 'oldGfan' : uses gfan version 0.2.1 or less @* 'findAll' : find all solutions of a zero-dimensional ideal over (w_1/w_0,...,w_n/w_0) @* 'noAbs' : do NOT use absolute primary decomposition @* 'puiseux' : n=1 and i is generated by one equation @* 'noResubst' : avoids the computation of the resubstitution RETURN: IF THE OPTION 'findAll' WAS NOT SET THEN: @* list, containing one lifting of the given point (w_1/w_0,...,w_n/w_0) in the tropical variety of i to a point in V(i) over Puiseux series field up to the first ord terms; more precisely: @* IF THE OPTION 'noAbs' WAS NOT SET, THEN: @* l[1] = ring Q[a]/m[[t]] @* l[2] = int @* l[3] = intvec @* l[4] = list @* IF THE OPTION 'noAbs' WAS SET, THEN: @* l[1] = ring Q[X(1),...,X(k)]/m[[t]] @* l[2] = int @* l[3] = intvec @* l[4] = list @* l[5] = string @* IF THE OPITON 'findAll' WAS SET, THEN: @* list, containing ALL liftings of the given point ((w_1/w_0,...,w_n/w_0) in the tropical variety of i to a point in V(i) over Puiseux series field up to the first ord terms, if the ideal is zero-dimensional over Q{{t}}; more precisely, each entry of the list is a list l as computed if 'findAll' was NOT set @* WE NOW DESCRIBE THE LIST ENTRIES IF 'findAll' WAS NOT SET: @* - the ring l[1] contains an ideal LIFT, which contains a point in V(i) lying over w up to the first ord terms; @* - and if the integer l[2] is N then t has to be replaced by t^1/N in the lift, or alternatively replace t by t^N in the defining ideal @* - if the k+1st entry of l[3] is non-zero, then the kth component of LIFT has to be multiplied t^(-l[3][k]/l[3][1]) AFTER substituting t by t^1/N @* - unless the option 'noResubst' was set, the kth entry of list l[4] is a string which represents the kth generator of the ideal i where the coordinates have been replaced by the result of the lift; the t-order of the kth entry should in principle be larger than the t-degree of LIFT @* - if the option 'noAbs' was set, then the string in l[5] defines a maximal ideal in the field Q[X(1),...,X(k)], where X(1),...,X(k) are the parameters of the ring in l[1]; the basefield of the ring in l[1] should be considered modulo this ideal REMARK: - it is best to use the procedure displayTropicalLifting to display the result @* - the option 'findAll' cannot be used if 'noAbs' is set @* - if the parameter 'findAll' is set AND the ideal i is zero-dimensional in Q{{t}}[x_1,...,x_n] then ALL points in V(i) lying over w are computed up to order ord; if the ideal is not-zero dimenisonal, then only the points in the ideal after cutting down to dimension zero will be computed @* - the procedure requires that the program GFAN is installed on your computer; if you have GFAN version less than 0.3.0 then you must use the optional parameter 'oldGfan' @* - the procedure requires the @sc{Singular} procedure absPrimdecGTZ to be present in the package primdec.lib, unless the option 'noAbs' is set; but even if absPrimdecGTZ is present it might be necessary to set the option 'noAbs' in order to avoid the costly absolute primary decomposition; the side effect is that the field extension which is computed throughout the recursion might need more than one parameter to be described @* - since Q is infinite, the procedure finishes with probability one @* - you can call the procedure with Z/pZ as base field instead of Q, but there are some problems you should be aware of: @* + the Puiseux series field over the algebraic closure of Z/pZ is NOT algebraicall closed, and thus there may not exist a point in V(i) over the Puiseux series field with the desired valuation; so there is no chance that the procedure produced a sensible output - e.g. if i=tx^p-tx-1 @* + if the dimension of i over Z/pZ(t) is not zero the process of reduction to zero might not work if the characteristic is small and you are unlucky @* + the option 'noAbs' has to be used since absolute primary decomposition in @sc{Singular} only works in characteristic zero @* - the basefield should either be Q or Z/pZ for some prime p; field extensions will be computed if necessary; if you need parameters or field extensions from the beginning they should rather be simulated as variables possibly adding their relations to the ideal; the weights for the additional variables should be zero EXAMPLE: example tropicalLifting; shows an example" { // if the basering has parameters, then exit with an error message if (parstr(basering)!="") { ERROR("The procedure is not implemented for rings with parameters. See: help tropicalLifting; for more information"); } // in order to avoid unpleasent surprises with the names of the variables // we change to a ring where the variables have names t and x(1),...,x(n) def ALTERRING=basering; if (nvars(basering)==2) { execute("ring BASERING=("+charstr(ALTERRING)+"),(t,x(1)),("+ordstr(ALTERRING)+");"); } else { execute("ring BASERING=("+charstr(ALTERRING)+"),(t,x(1.."+string(nvars(ALTERRING)-1)+")),("+ordstr(ALTERRING)+");"); } map altphi=ALTERRING,maxideal(1); ideal i=altphi(i); int j,k,l,jj,kk; // index variables // work through the optionional parameters int isprime,iszerodim,isintrop,gfanold,findall,noabs,nogfan,noresubst,puiseux; for (j=1;j<=size(#);j++) { if (#[j]=="isZeroDimensional") { iszerodim=1; } if (#[j]=="isPrime") { isprime=1; } if (#[j]=="isInTrop") { isintrop=1; } if (#[j]=="oldGfan") { gfanold=1; } if (#[j]=="findAll") { findall=1; } if (#[j]=="noAbs") { noabs=1; } if (#[j]=="puiseux") { int puiseuxgood=0; // test, if the option puiseux makes sense, i.e. we are considering a plane curve if ((size(i)==1) and (nvars(basering)==2)) { puiseuxgood=1; } // the case, where the base field has a parameter and a minimal polynomial is // treated by an additional variable (which should be the last variable) // and an ideal containing the minimal polynomial as first entry if ((size(i)==2) and (nvars(basering)==3)) { // check if the first entry is the minimal polynomial poly mpcheck=i[1]; if (substitute(mpcheck,var(1),0,var(2),0)==mpcheck) { puiseuxgood=1; } kill mpcheck; } if (puiseuxgood==0) { ERROR("The option puiseux is not allowed for this ring. See: help tropicalLifting; for more information"); } puiseux=1; nogfan=1; // if puiseux is set, then gfan should not be used } // this option is not documented -- it prevents the execution of gfan and // just asks for wneu to be inserted -- it can be used to check problems // with the precedure without calling gfan, if wneu is know from previous // computations if (#[j]=="noGfan") { nogfan=1; } if (#[j]=="noResubst") { noresubst=1; } } // if the basering has characteristic not equal to zero, // then absolute factorisation // is not available, and thus we need the option noAbs /* if ((char(basering)!=0) and (noabs!=1)) { ERROR("If the characteristic is not zero the procedure should be called with option 'noAbs'"); } */ int N=1; // we are working with the variable t^1/N - the N may be changed // w_0 must be non-zero! if (w[1]==0) { Error("The first coordinate of your input w must be NON-ZERO, since it is a DENOMINATOR!"); } // if w_0<0, then replace w by -w, so that the "denominator" w_0 is positive if (w[1]<0) { w=-w; } intvec prew=w; // stores w for later reference // for our computations, w[1] represents the weight of t and this // should be -w_0 !!! w[1]=-w[1]; // if w_0!=-1 then replace t by t^-w_0 and w_0 by -1 if (w[1]<-1) { i=subst(i,t,t^(-w[1])); N=-w[1]; w[1]=-1; } // if some entry of w is positive, we have to make a transformation, // which moves it to something non-positive for (j=2;j<=nvars(basering);j++) { if (w[j]>0) { // transform x_j to t^(-w[j])*x_j in the ideal i i=phiOmega(i,j,w[j]); w[j]=0; } } prew[1]=prew[1]+w[1]; prew=prew-w; // this now contains the positive part of the original w, // but the original first comp. of w // pass to a ring which represents Q[t]_[x1,...,xn] // for this, unfortunately, t has to be the last variable !!! ideal variablen; for (j=2;j<=nvars(basering);j++) { variablen=variablen+var(j); } execute("ring GRUNDRING=("+charstr(basering)+"),("+string(variablen)+",t),(dp("+string(nvars(basering)-1)+"),lp(1));"); ideal variablen; for (j=1;j<=nvars(basering)-1;j++) { variablen=variablen+var(j); } map GRUNDPHI=BASERING,t,variablen; ideal i=GRUNDPHI(i); // compute the initial ideal of i and test if w is in the tropical // variety of i // - the last entry 1 only means that t is the last variable in the ring ideal ini=tInitialIdeal(i,w,1); if (isintrop==0) // test if w is in trop(i) only if isInTrop has not been set { poly product=1; for (j=1;j<=nvars(basering)-1;j++) { product=product*var(j); } if (radicalMemberShip(product,ini)==1) // if w is not in Trop(i) - error { ERROR("The integer vector is not in the tropical variety of the ideal."); } } // compute next the dimension of i in K(t)[x] and cut the ideal down to dim 0 if (iszerodim==0) // do so only if is_dim_zero is not set { execute("ring QUOTRING=("+charstr(basering)+",t),("+string(variablen)+"),dp;"); ideal i=groebner(imap(GRUNDRING,i)); int dd=dim(i); setring GRUNDRING; // if the dimension is not zero, we cut the ideal down to dimension zero // and compute the // t-initial ideal of the new ideal at the same time if(dd!=0) { // the procedurce cutdown computes a new ring, in which there lives a // zero-dimensional // ideal which has been computed by cutting down the input with // generic linear forms // of the type x_i1-p_1,...,x_id-p_d for some polynomials // p_1,...,p_d not depending // on the variables x_i1,...,x_id; that way we have reduced // the number of variables by dd !!! // the new zero-dimensional ideal is called i, its t-initial // ideal (with respect to // the new w=CUTDOWN[2]) is ini, and finally there is a list // repl in the ring // which contains at the polynomial p_j at position i_j and //a zero otherwise; if (isprime==0) // the minimal associated primes of i are computed { list CUTDOWN=cutdown(i,w,dd); } else // i is assumed to be prime { list CUTDOWN=cutdown(i,w,dd,"isPrime"); } def CUTDOWNRING=CUTDOWN[1]; intvec precutdownw=w; // save the old w for later reference w=CUTDOWN[2]; // the new w - some components have been eliminated setring CUTDOWNRING; } } else { int dd=0; // the dimension of i } list liftrings; // will contain the final result // if the procedure is called without 'findAll' then it may happen, that no // proper solution is found when dd>0; in that case we have // to start all over again; // this is controlled by the while-loop while (size(liftrings)==0) { // compute lifting for the zero-dimensional ideal via tropicalparametrise if (noabs==1) // do not use absolute primary decomposition { list TP=list(tropicalparametriseNoabs(i,w,ordnung,gfanold,nogfan,puiseux,ini)); } else // use absolute primary decomposition { list TP=tropicalparametrise(i,w,ordnung,gfanold,findall,nogfan,puiseux,ini); } // compute the liftrings by resubstitution kk=1; // counts the liftrings int isgood; // test in the non-zerodimensional case // if the result has the correct valuation for (jj=1;jj<=size(TP);jj++) { // the list TP contains as a first entry the ring over which the // tropical parametrisation // of the (possibly cutdown ideal) i lives def LIFTRING=TP[jj][1]; // if the dimension of i originally was not zero, // then we have to fill in the missing // parts of the parametrisation if (dd!=0) { // we need a ring where the parameters X_1,...,X_k // from LIFTRING are present, // and where also the variables of CUTDOWNRING live execute("ring REPLACEMENTRING=("+charstr(LIFTRING)+"),("+varstr(CUTDOWNRING)+"),dp;"); list repl=imap(CUTDOWNRING,repl); // get the replacement rules // from CUTDOWNRING ideal PARA=imap(LIFTRING,PARA); // get the zero-dim. parametrisatio // from LIFTRING // compute the lift of the solution of the original ideal i ideal LIFT; k=1; // the lift has as many components as GRUNDRING has variables!=t for (j=1;j<=nvars(GRUNDRING)-1;j++) { // if repl[j]=0, then the corresponding variable was not eliminated if (repl[j]==0) { LIFT[j]=PARA[k]; // thus the lift has been // computed by tropicalparametrise k++; // k checks how many entries of PARA have already been used } else // if repl[j]!=0, repl[j] contains replacement rule for the lift { LIFT[j]=repl[j]; // we still have to replace the vars // in repl[j] by the corresp. entries of PARA // replace all variables!=t (from CUTDOWNRING) for (l=1;l<=nvars(CUTDOWNRING)-1;l++) { // substitute the kth variable by PARA[k] LIFT[j]=subst(LIFT[j],var(l),PARA[l]); } } } setring LIFTRING; ideal LIFT=imap(REPLACEMENTRING,LIFT); // test now if the LIFT has the correct valuation !!! // note: it may happen, that when resubstituting PARA into // the replacement rules // there occured some unexpected cancellation; // we only know that for SOME // solution of the zero-dimensional reduction NO // canellation will occur, // but for others this may very well happen; // this in particular means that // we possibly MUST compute all zero-dimensional // solutions when cutting down! intvec testw=precutdownw[1]; for (j=1;j<=size(LIFT);j++) { testw[j+1]=-ord(LIFT[j]); } if (testw==precutdownw) { isgood=1; } else { isgood=0; } } else { setring LIFTRING; ideal LIFT=PARA; isgood=1; } kill PARA; // only if LIFT has the right valuation we have to do something if (isgood==1) { // it remains to reverse the original substitutions, // where appropriate !!! // if some entry of the original w was positive, // we replace the corresponding // variable x_i by t^-w[i]*x_i, so we must now replace // the corresponding entry // LIFT[i] by t^-w[i]*LIFT[i] --- the original w is stored in prew // PROBLEM: THIS CANNOT BE DONE SINCE -w[i] IS NEGATIVE // INSTEAD: RETURN prew IN THE LIST /* for (j=2;j<=size(prew);j++) { if (prew[j]>0) { LIFT[j-1]=t^(-prew[j])*LIFT[j-1]; } } */ // if LIFTRING contains a parameter @a, change it to a if ((noabs==0) and (defined(@a)==-1)) { // pass first to a ring where a and @a // are variables in order to use maps string @mp= string(minpoly); ring INTERRING=char(LIFTRING),(t,@a,a),dp; execute("poly mp=" + @mp + ";"); ideal LIFT=imap(LIFTRING,LIFT); kill LIFTRING; // replace @a by a in minpoly and in LIFT map phi=INTERRING,t,a,a; mp=phi(mp); LIFT=phi(LIFT); // pass now to a ring whithout @a and with a as parameter ring LIFTRING=(char(INTERRING),a),t,ls; minpoly=number(imap(INTERRING,mp)); ideal LIFT=imap(INTERRING,LIFT); kill INTERRING; } // then export LIFT export(LIFT); // test the result by resubstitution setring GRUNDRING; list resubst; if (noresubst==0) { if (noabs==1) { resubst=tropicalliftingresubstitute(substitute(i,t,t^(TP[jj][2])),list(LIFTRING),N*TP[jj][2],TP[jj][3]); } else { resubst=tropicalliftingresubstitute(substitute(i,t,t^(TP[jj][2])),list(LIFTRING),N*TP[jj][2]); } } setring BASERING; // Finally, if t has been replaced by t^N, then we have to change the // third entry of TP by multiplying by N. if (noabs==1) { liftrings[kk]=list(LIFTRING,N*TP[jj][2],prew,resubst,TP[jj][3]); } else { liftrings[kk]=list(LIFTRING,N*TP[jj][2],prew,resubst); } kk++; kill resubst; } setring BASERING; kill LIFTRING; } // if dd!=0 and the procedure was called without the // option findAll, then it might very well // be the case that no solution is found, since // only one solution for the zero-dimensional // reduction was computed and this one might have // had cancellations when resubstituting; // if so we have to restart the process with the option findAll if (size(liftrings)==0) { "The procedure was called without findAll and no proper solution was found."; "The procedure will be restarted with the option 'findAll'."; "Go on by hitting RETURN!"; findall=1; noabs=0; setring CUTDOWNRING; int hadproblems; "i";i; "ini";tInitialIdeal(i,w,1); /* setring GRUNDRING; list repl=imap(CUTDOWNRING,repl); i=imap(CUTDOWNRING,i); for (j=1;j<=nvars(basering)-1;j++) { if (repl[j]!=0) { i=i+ideal(var(j)-repl[j]); } } ini=tInitialIdeal(i,precutdownw,1); w=precutdownw; */ ~ } } // if internally the option findall was set, then return // only the first solution if (defined(hadproblems)!=0) { liftrings=liftrings[1]; } /////////////////////////////////////////////////////////// if (voice+printlevel>=2) { "The procedure has created a list of lists. The jth entry of this list contains a ring, an integer and an intvec. In this ring lives an ideal representing the wanted lifting, if the integer is N then in the parametrisation t has to be replaced by t^1/N, and if the ith component of the intvec is w[i] then the ith component in LIFT should be multiplied by t^-w[i]/N in order to get the parametrisation. Suppose your list has the name L, then you can access the 1st ring via: "; if (findall==1) { "def LIFTRing=L[1][1]; setring LIFTRing; LIFT; "; } else { "def LIFTRing=L[1]; setring LIFTRing; LIFT; "; } } if (findall==1) // if all solutions have been computed, return a list of lists { return(liftrings); } else //if only 1 solution was to be computed, return the 1st list in liftrings { liftrings=liftrings[1]; return(liftrings); } } example { "EXAMPLE:"; echo=2; //////////////////////////////////////////////////////// // 1st EXAMPLE: //////////////////////////////////////////////////////// ring r=0,(t,x),dp; ideal i=(1+t2)*x2+t5x+t2; intvec w=1,-1; list LIST=tropicalLifting(i,w,4); def LIFTRing=LIST[1]; setring LIFTRing; // LIFT contains the first 4 terms of a point in the variety of i // over the Puiseux series field C{{t}} whose order is -w[1]/w[0]=1 LIFT; // Since the computations were done over Q a field extension was necessary, // and the parameter "a" satisfies the equation given by minpoly minpoly; //////////////////////////////////////////////////////// // 2nd EXAMPLE //////////////////////////////////////////////////////// setring r; LIST=tropicalLifting(x12-t11,intvec(12,-11),2,"isPrime","isInTrop"); def LIFTRing2=LIST[1]; setring LIFTRing2; // This time, LIFT contains the lifting of the point -w[1]/w[0]=11/12 // only after we replace in LIFT the variable t by t^1/N with N=LIST[3] LIFT; LIST[3]; /////////////////////////////////////////////////////// // 3rd EXAMPLE //////////////////////////////////////////////////////// ring R=0,(t,x,y,z),dp; ideal i=-y2t4+x2,yt3+xz+y; w=1,-2,0,2; LIST=tropicalLifting(i,w,3); // This time, LIFT contains the lifting of the point v=(-2,0,2) // only after we multiply LIFT[3] by t^k with k=-LIST[4][3]; // NOTE: since the last component of v is positive, the lifting // must start with a negative power of t, which in Singular // is not allowed for a variable. def LIFTRing3=LIST[1]; setring LIFTRing3; LIFT; LIST[4]; // An easier way to display this is via displayTropicalLifting. setring R; displayTropicalLifting(LIST,"subst"); } /////////////////////////////////////////////////////////////////////////////// proc displayTropicalLifting (list troplift,list #) "USAGE: displaytropcallifting(troplift[,#]); troplift list, # list ASSUME: troplift is the output of tropicalLifting; the optional parameter # can be the string 'subst' RETURN: none NOTE: - the procedure displays the output of the procedure tropicalLifting @* - if the optional parameter 'subst' is given, then the lifting is substituted into the ideal and the result is displayed EXAMPLE: example displayTropicalLifting; shows an example" { int j; // if the procedure has to display more than one lifting if (typeof(troplift[1])=="list") { for (j=1;j<=size(troplift);j++) { "============================="; string(j)+". Lifting:"; ""; displayTropicalLifting(troplift[j],#); ""; } } // if the procedure has to display only one lifting else { list variablen; for (j=2;j<=nvars(basering);j++) { variablen[j-1]=string(var(j)); } def LIFTRing=troplift[1]; int N=troplift[2]; intvec wdiff=troplift[3]; string LIFTpar=parstr(LIFTRing); if (char(LIFTRing)==0) { string Kstring="Q"; } else { string Kstring="Z/"+string(char(LIFTRing))+"Z"; } // this means that tropicalLifting was called with // absolute primary decomposition if (size(troplift)==4) { setring LIFTRing; "The lifting of the point in the tropical variety lives in the ring"; if ((size(LIFTpar)==0) and (N==1)) { Kstring+"[[t]]"; } if ((size(LIFTpar)==0) and (N!=1)) { Kstring+"[[t^(1/"+string(N)+")]]"; } if ((size(LIFTpar)!=0) and (N!=1)) { Kstring+"["+LIFTpar+"]/"+string(minpoly)+"[[t^(1/"+string(N)+")]]"; } if ((size(LIFTpar)!=0) and (N==1)) { Kstring+"["+LIFTpar+"]/"+string(minpoly)+"[[t]]"; } } else // tropicalLifting was called with the option noAbs { string m=troplift[5]; setring LIFTRing; "The lifting of the point in the tropical variety lives in the ring"; if ((size(LIFTpar)==0) and (N==1)) { Kstring+"[[t]]"; } if ((size(LIFTpar)==0) and (N!=1)) { Kstring+"[[t^(1/"+string(N)+")]]"; } if ((size(LIFTpar)!=0) and (N!=1)) { Kstring+"["+LIFTpar+"]/M[[t^(1/"+string(N)+")]]"; "where M is the maximal ideal"; "M=<"+m+">"; } if ((size(LIFTpar)!=0) and (N==1)) { Kstring+"["+LIFTpar+"]/M[[t]]"; "where M is the maximal ideal"; "M=<"+m+">"; } } ""; "The lifting has the form:"; for (j=1;j<=size(LIFT);j++) { if (ncols(LIFT)==size(variablen)) { variablen[j]+"="+displaypoly(LIFT[j],N,wdiff[j+1],wdiff[1]); } else { "var("+string(j)+")="+displaypoly(LIFT[j],N,wdiff[j+1],wdiff[1]); } } if (size(#)>0) { if (#[1]=="subst") { ""; "Substituting the solution into the ideal gives:"; for(j=1;j<=size(troplift[4]);j++) { "i["+string(j)+"]="+troplift[4][j]; } } } } } example { "EXAMPLE:"; echo=2; ring r=0,(t,x,y,z),dp; ideal i=-y2t4+x2,yt3+xz+y; intvec w=2,-4,0,4; displayTropicalLifting(tropicalLifting(i,w,3),"subst"); } proc puiseuxExpansion (poly f,int n,list #) "USAGE: puiseuxExpansion(f,n,#); f poly, n int, # list ASSUME: f is a non-constant polynomial in two variables which is not divisible by the first variable and which is squarefree as a power series over the complex numbers; the base field is either the field of rational numbers or a finite extension thereof; monomial ordering is assumed to be local; the optional parameter # can be the string 'subst' RETURN: list, where each entry of the list l describes the Newton-Puiseux @* parametrisations of one branch of the plane curve singularity @* at the origin defined by f; only the first n terms of each @* parametetrisation are computed @* l[i][1] = is a ring @* l[i][2] = int @* l[i][3] = string @* @* WE NOW DESCRIBE THE LIST ENTRIES l[i] IN MORE DETAIL: @* - the ring l[i][1] contains an ideal LIFT and the Newton-Puiseux parametrisation of the branch is given by x=t^N and y=LIFT[1], where N=l[i][2] @* - if the base field had a parameter and a minimal polynomial, then the new base field will have a parameter and a new minimal polynomial, and LIFT[2] describes how the old parameter can be computed from the new one @* - if the option subst was set l[i][3] contains the polynomial where y has been substituted by y(t^{1/N}) as a string REMARK: - it is best to use the procedure displayPuiseuxExpansion to display the result @* - the procedure requires the @sc{Singular} procedure absPrimdecGTZ to be present in the package primdec.lib @* - if f is not squarefree it will be replaced by its squarefree part EXAMPLE: example puiseuxExpansion; shows an example" { if (deg(f)<=0) { ERROR("The input polynomial is not allowed to be constant!"); } if (char(basering)!=0) { ERROR("In positive characteristic a Puiseux expansion will not in general exist!"); } if ((npars(basering)>1) or ((npars(basering)==1) and (minpoly==0))) { ERROR("The procedure is not implemented for this base field!"); } if (nvars(basering)!=2) { ERROR("The base ring should depend on exactly two variables."); } if (f==(f / var(1))*var(1)) { ERROR("The input polynomial is not allowed to be divisible by the first variable."); } // if the ordering is not local, change to a local ordering if ((11) { if (#[1]=="subst") { noResubst=0; } } // replace f by its squarefree part f=squarefree(f); // check if var(1) or var(2) divide f int coordinatebranchtwo; if (f==(f / var(2))*var(2)) { coordinatebranchtwo=1; f=f/var(2); } int jj; // if f is now constant then we should skip the next part if (deg(f)!=0) { // compute the Newton polygon int zw; intvec w; int ggteiler; list NewtP=newtonpoly(f); list tls; list pexp; // if the base field has a minimal polynomial change the base ring and the ideal if (minpoly!=0) { string @mp= string(minpoly); def OLDRING=basering; execute("ring NEWRING=0,("+varstr(basering)+","+parstr(basering)+"),ds;"); execute("ideal I = " + @mp + ";"); I = I,imap(OLDRING,f); } else { ideal I=f; } // for each facet of the Newton polygon compute the corresponding parametrisations // using tropicalLifting with the option "puiseux" which avoids gfan for (jj=1;jj<=size(NewtP)-1;jj++) { w=NewtP[jj]-NewtP[jj+1]; ggteiler=gcd(w[1],w[2]); zw=w[1]/ggteiler; w[1]=w[2]/ggteiler; w[2]=zw; // if we have introduced a third variable for the parameter, then w needs a third component if (nvars(basering)==3) { w[3]=0; } if (noResubst==0) { tls=tropicalLifting(I,w,n,"findAll","puiseux"); } else { tls=tropicalLifting(I,w,n,"findAll","puiseux","noResubst"); } pexp=pexp+tls; } // kill rings that are no longer needed if (defined(NEWRING)) { setring OLDRING; kill NEWRING; } if (defined(GLOBALRING)) { setring GLOBALRING; kill LOCALRING; } // remove the third entry in the list of parametrisations since we know // that none of the exponents in the parametrisation will be negative for (jj=1;jj<=size(pexp);jj++) { pexp[jj]=delete(pexp[jj],3); pexp[jj][3]=pexp[jj][3][size(pexp[jj][3])]; // replace the list in pexp[jj][3] by its first entry } } else // if f was reduced to a constant we still have to introduce the list pexp { list pexp; } // we have to add the parametrisations for the branches var(1) resp. var(2) if // we removed them def BASERING=basering; if (coordinatebranchtwo==1) { ring brring=0,t,ds; ideal LIFT=0; export(LIFT); setring BASERING; list pexpentry; pexpentry[1]=brring; pexpentry[2]=1; pexpentry[3]="0"; pexp[size(pexp)+1]=pexpentry; kill brring; } // check if all branches have been found int numberofbranchesfound; for (jj=1;jj<=size(pexp);jj++) { def countring=pexp[jj][1]; setring countring; if (minpoly==0) { numberofbranchesfound=numberofbranchesfound+1; } else { string @mp= string(minpoly); ring degreering=0,a,dp; execute("poly mp=" + @mp + ";"); numberofbranchesfound=numberofbranchesfound+deg(mp); setring countring; kill degreering; } setring BASERING; kill countring; } // give a warning if not all branches have been found if (numberofbranchesfound!=ord(subst(f,x,0))+coordinatebranchtwo) { "!!!! WARNING: The number of terms computed in the Puiseux expansion were"; "!!!! not enough to find all branches of the curve singularity!"; } return(pexp); } example { "EXAMPLE:"; echo=2; ring r=0,(x,y),ds; poly f=x2-y4+x5y7; puiseuxExpansion(puiseuxExpansion(f,3)); } proc displayPuiseuxExpansion (list puiseux,list #) "USAGE: displayPuiseuxExpansion(puiseux[,#]); puiseux list, # list ASSUME: puiseux is the output of puiseuxExpansion; the optional parameter # can be the string 'subst' RETURN: none NOTE: - the procedure displays the output of the procedure puiseuxExpansion @* - if the optional parameter 'subst' is given, then the expansion is substituted into the polynomial and the result is displayed @* - if the base field had a parameter and a minimal polynomial, then the new base field will have a parameter and a minimal polynomial; var(2) is the old parameter and it is displayed how the old parameter can be computed from the new one EXAMPLE: example displayPuiseuxExpansion; shows an example" { int j; // if the procedure has to display more than one expansion if (typeof(puiseux[1])=="list") { for (j=1;j<=size(puiseux);j++) { "============================="; string(j)+". Expansion:"; ""; displayPuiseuxExpansion(puiseux[j],#); ""; } } // if the procedure has to display only one expansion else { list variablen; for (j=2;j<=nvars(basering);j++) { variablen[j-1]=string(var(j)); } def BASERING=basering; def LIFTRing=puiseux[1]; int N=puiseux[2]; string LIFTpar=parstr(LIFTRing); string Kstring="Q"; setring LIFTRing; "The Puiseux expansion lives in the ring"; if ((size(LIFTpar)==0) and (N==1)) { Kstring+"[[t]]"; } if ((size(LIFTpar)==0) and (N!=1)) { Kstring+"[[t^(1/"+string(N)+")]]"; } if ((size(LIFTpar)!=0) and (N!=1)) { Kstring+"["+LIFTpar+"]/"+string(minpoly)+"[[t^(1/"+string(N)+")]]"; } if ((size(LIFTpar)!=0) and (N==1)) { Kstring+"["+LIFTpar+"]/"+string(minpoly)+"[[t]]"; } ""; "The expansion has the form:"; // treat the case LIFT==0 separately if (size(LIFT)==0) { variablen[1]+"=0"; } else { for (j=1;j<=size(LIFT);j++) { if (ncols(LIFT)==size(variablen)) { variablen[j]+"="+displaypoly(LIFT[j],N,0,1); } else { "var("+string(j)+")="+displaypoly(LIFT[j],N,0,1); } } } if (size(#)>0) { if (#[1]=="subst") { setring BASERING; ""; "Substituting the expansion into the polynomial gives:"; "f="+puiseux[3]; } } } } example { "EXAMPLE:"; echo=2; ring r=0,(x,y),ds; poly f=x2-y4+x5y7; displayPuiseuxExpansion(puiseuxExpansion(f,3)); } /////////////////////////////////////////////////////////////////////////////// /// Procedures concerned with drawing a tropical curve or a Newton subdivision /////////////////////////////////////////////////////////////////////////////// proc tropicalCurve (def tp,list #) "USAGE: tropicalCurve(tp[,#]); tp list, # optional list ASSUME: tp is list of linear polynomials of the form ax+by+c with integers a, b and a rational number c representing a tropical Laurent polynomial defining a tropical plane curve; alternatively tp can be a polynomial in Q(t)[x,y] defining a tropical plane curve via the valuation map; the basering must have a global monomial ordering, two variables and up to one parameter! RETURN: list, each entry i=1,...,size(l)-1 corresponds to a vertex in the tropical plane curve defined by tp l[i][1] = x-coordinate of the ith vertex l[i][2] = y-coordinate of the ith vertex l[i][3] = intmat, if j is an entry in the first row of intmat then the ith vertex of the tropical curve is connected to the jth vertex with multiplicity given by the corresponding entry in the second row l[i][4] = list of lists, the first entry of a list is a primitive integer vector defining the direction of an unbounded edge emerging from the ith vertex of the graph, the corresponding second entry in the list is the multiplicity of the unbounded edge l[i][5] = a polynomial whose monomials mark the vertices in the Newton polygon corresponding to the entries in tp which take the common minimum at the ith vertex -- if some coefficient a or b of the linear polynomials in the input was negative, then each monomial has to be shifted by the values in l[size(l)][3] l[size(l)][1] = list, the entries describe the boundary points of the Newton subdivision l[size(l)][2] = list, the entries are pairs of integer vectors defining an interior edge of the Newton subdivision l[size(l)][3] = intvec, the monmials occuring in l[i][5] have to be shifted by this vector in order to represent marked vertices in the Newton polygon NOTE: here the tropical polynomial is supposed to be the MINIMUM of the linear forms in tp, unless the optional input #[1] is the string 'max' EXAMPLE: example tropicalCurve; shows an example" { // introduce necessary variables int i,j,k,l,d; intvec v,w,u; intmat D[2][2]; // the basering must have a global monomial ordering if (1>var(1)) { ERROR("The basering should have a global monomial ordering, e.g. ring r=(0,t),(x,y),dp;"); } // if you insert a single polynomial instead of an ideal // representing a tropicalised polynomial, // then we compute first the tropicalisation of this polynomial // -- this feature is not documented in the above help string if (typeof(tp)=="poly") { // exclude the case that the basering has not precisely // one parameter and two indeterminates if ((npars(basering)!=1) or (nvars(basering)!=2)) { ERROR("The basering should have precisely one parameter and two indeterminates!"); } poly f=tp; kill tp; list tp=tropicalise(f,#); // the tropicalisation of f } // Exclude the case that the basering has more than 2 indeterminates if (nvars(basering) != 2) { ERROR("The basering should have precisely two indeterminates!"); } // -1) Exclude the pathological case that the defining // tropical polynomial has only one term, // so that the tropical variety is not defined. if (size(tp)==1) { ERROR("A monomial does not define a tropical curve!"); intmat M[2][1]=0,0; return(list(list(0,0,M,list(),detropicalise(tp[1])),list(list(leadexp(detropicalise(tp[1]))),list()))); } // 0) If the input was a list of linear polynomials, // then some coefficient of x or y can be negative, // i.e. the input corresponds to the tropical curve // of a Laurent polynomial. In that case we should // add some ax+by, so that all coefficients are positive. // This does not change the tropical curve. // however, we have to save (a,b), since the Newton // polygone has to be shifted by (-a,-b). poly aa,bb; // koeffizienten for (i=1;i<=size(tp);i++) { if (koeffizienten(tp[i],1)=1;k--) { deleted=0; for (l=size(pairs[j]);l>=1 and deleted==0;l--) { 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]))) { inneredges[d]=pairs[i][k]; // new inner edge is saved in inneredges d++; ggt=abs(gcd(pairs[i][k][1][1]-pairs[i][k][2][1],pairs[i][k][1][2]-pairs[i][k][2][2])); zwsp=j,ggt; // and it is recorded that the ith and jth // vertex should be connected with mult ggt vtp[i][6]=intmatconcat(vtp[i][6],zwsp); zwsp=i,ggt; vtp[j][6]=intmatconcat(vtp[j][6],zwsp); pairs[i]=delete(pairs[i],k); // finally the pair is deleted // from both sets of pairs pairs[j]=delete(pairs[j],l); deleted=1; } } } } } // 5) The entries in vtp[i][6] are ordered, multiple entries are removed, // and the redundant zero is removed as well for (i=1;i<=size(vtp);i++) { vtp[i][6]=clearintmat(vtp[i][6]); } // 6) Declare the orientation of the boundary of the Newton polytope. // 6.1) Collect all potential vertices of the boundary of the Newton polytope. list vertices; // all vertices in the set of pairs k=1; for (i=1;i<=size(pairs);i++) { for (j=1;j<=size(pairs[i]);j++) { vertices[k]=pairs[i][j][1]; k++; vertices[k]=pairs[i][j][2]; k++; } } // 6.2) delete multiple vertices for (i=size(vertices)-1;i>=1;i--) { for (j=size(vertices);j>=i+1;j--) { if (vertices[i]==vertices[j]) { vertices=delete(vertices,j); } } } // 6.3) Order the vertices such that passing from one to the next we // travel along the boundary of the Newton polytope clock wise. boundaryNSD=findOrientedBoundary(vertices); list orderedvertices=boundaryNSD[1]; // 7) Find the unbounded edges emerging from a vertex in the tropical curve. // For this we check the remaining pairs for the ith NSD. // Each pair is ordered according // to the order in which the vertices occur in orderedvertices. // The direction of the // unbounded edge is then the outward pointing primitive normal // vector to the vector // pointing from the first vertex in a pair to the second one. intvec normalvector; // stores the outward pointing normal vector intvec zwspp; // used for intermediate storage int zw,pos1,pos2; // stores the gcd of entries of the // non-normalised normal vector, etc. int gestorben; // tests if unbounded edges are multiple for (i=1;i<=size(pairs);i++) { list ubedges; // stores the unbounded edges k=1; // counts the unbounded edges for (j=1;j<=size(pairs[i]);j++) { // computes the position of the vertices in the pos1=positionInList(orderedvertices,pairs[i][j][1]); // pair in the list orderedvertices pos2=positionInList(orderedvertices,pairs[i][j][2]); if (((pos1>pos2) and !((pos1==size(orderedvertices)) and (pos2==1))) or ((pos2==size(orderedvertices)) and (pos1==1))) // reorders them if necessary { zwspp=pairs[i][j][1]; pairs[i][j][1]=pairs[i][j][2]; pairs[i][j][2]=zwspp; } // the vector pointing from vertex 1 in the pair to vertex2 normalvector=pairs[i][j][2]-pairs[i][j][1]; ggt=gcd(normalvector[1],normalvector[2]); // the gcd of the entries zw=normalvector[2]; // create the outward pointing normal vector normalvector[2]=-normalvector[1]/ggt; normalvector[1]=zw/ggt; if (size(#)==0) // we are computing w.r.t. minimum { ubedges[k]=list(normalvector,ggt); // store outward pointing normal vec. } else // we are computing w.r.t. maximum { ubedges[k]=list(-normalvector,ggt); //store outward pointing normal vec. } k++; } // remove multiple unbounded edges for (j=size(ubedges);j>=1;j--) { gestorben=0; for(k=1;(k<=j-1) and (gestorben==0);k++) { if (ubedges[k][1]==ubedges[j][1]) { ubedges[k][2]=ubedges[k][2]+ubedges[j][2]; // add the multiplicities ubedges=delete(ubedges,j); gestorben=1; } } } vtp[i][7]=ubedges; // store the unbounded edges in vtp[i][7] kill ubedges; } // 8) Store the computed information for the ith part // of the NSD in the list graph[i]. list graph,gr; for (i=1;i<=size(vtp);i++) { // the first coordinate of the ith vertex of the tropical curve gr[1]=vtp[i][1]; // the second coordinate of the ith vertex of the tropical curve gr[2]=vtp[i][2]; // to which vertices is the ith vertex of the tropical curve connected gr[3]=vtp[i][6]; // the directions unbounded edges emerging from the ith // vertex of the trop. curve gr[4]=vtp[i][7]; // the vertices of the boundary of the ith part of the NSD gr[5]=vtp[i][3]; graph[i]=gr; } // 9) Shift the Newton subdivision by (aa,bb) if necessary intvec shiftvector=intvec(int(aa),int(bb)); if ((aa!=0) or (bb!=0)) { for (i=1;i<=size(boundaryNSD[2]);i++) { boundaryNSD[2][i]=boundaryNSD[2][i]+shiftvector; } for (i=1;i<=size(inneredges);i++) { for (j=1;j<=size(inneredges[i]);j++) { inneredges[i][j]=inneredges[i][j]+shiftvector; } } } // 10) Finally store the boundary vertices and // the inner edges as last entry in the list graph. // This represents the NSD. graph[size(vtp)+1]=list(boundaryNSD[2],inneredges,shiftvector); return(graph); } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),dp; poly f=t*(x7+y7+1)+1/t*(x4+y4+x2+y2+x3y+xy3)+1/t7*x2y2; list graph=tropicalCurve(f); // the tropical curve has size(graph)-1 vertices size(graph)-1; // the coordinates of the first vertex are graph[1][1],graph[1][2]; graph[1][1],graph[1][2]; // the first vertex is connected to the vertices // graph[1][3][1,1..ncols(graph[1][3])] intmat M=graph[1][3]; M[1,1..ncols(graph[1][3])]; // the weights of the edges to these vertices are // graph[1][3][2,1..ncols(graph[1][3])] M[2,1..ncols(graph[1][3])]; // from the first vertex emerge size(graph[1][4]) unbounded edges size(graph[1][4]); // the primitive integral direction vector of the first unbounded edge // of the first vertex graph[1][4][1][1]; // the weight of the first unbounded edge of the first vertex graph[1][4][1][2]; // the monomials which are part of the Newton subdivision of the first vertex graph[1][5]; // connecting the points in graph[size(graph)][1] we get // the boundary of the Newton polytope graph[size(graph)][1]; // an entry in graph[size(graph)][2] is a pair of points // in the Newton polytope bounding an inner edge graph[size(graph)][2][1]; } ///////////////////////////////////////////////////////////////////////// proc drawTropicalCurve (def f,list #) "USAGE: drawTropicalCurve(f[,#]); f poly or list, # optional list ASSUME: f is list of linear polynomials of the form ax+by+c with integers a, b and a rational number c representing a tropical Laurent polynomial defining a tropical plane curve; alternatively f can be a polynomial in Q(t)[x,y] defining a tropical plane curve via the valuation map; the basering must have a global monomial ordering, two variables and up to one parameter! RETURN: NONE NOTE: - the procedure creates the files /tmp/tropicalcurveNUMBER.tex and /tmp/tropicalcurveNUMBER.ps, where NUMBER is a random four digit integer; moreover it displays the tropical curve via kghostview; if you wish to remove all these files from /tmp, call the procedure cleanTmp @* - edges with multiplicity greater than one carry this multiplicity @* - if # is empty, then the tropical curve is computed w.r.t. minimum, if #[1] is the string 'max', then it is computed w.r.t. maximum @* - if the last optional argument is 'onlytexfile' then only the latex file is produced; this option should be used if kghostview is not installed on your system @* - note that lattice points in the Newton subdivision which are black correspond to markings of the marked subdivision, while lattice points in grey are not marked EXAMPLE: example drawTropicalCurve shows an example" { // check if the option "onlytexfile" is set, then only a tex file is produced if (size(#)!=0) { if (#[size(#)]=="onlytexfile") { int onlytexfile; #=delete(#,size(#)); } } // start the actual computations string texf; int j; if (typeof(f)=="poly") { // exclude the case that the basering has not precisely // one parameter and two indeterminates if ((npars(basering)!=1) or (nvars(basering)!=2)) { ERROR("The basering should have precisely one parameter and two indeterminates!"); } // if the characteristic of the base field is not 0 then replace the base field // by a field of characteristic zero if (char(basering)!=0) { string polynomstring=string(f); execute("ring drawring=(0,"+parstr(basering)+"),("+varstr(basering)+"),dp;"); execute("poly f="+polynomstring+";"); } texf=texPolynomial(f); // write the polynomial over Q(t) list graph=tropicalCurve(tropicalise(f,#),#); // graph of tropicalis. of f } if (typeof(f)=="list") { if (typeof(f[1])=="list") { texf="\\mbox{\\tt The defining equation was not handed over!}"; list graph=f; } else { // write the tropical polynomial defined by f if (size(#)==0) { texf="\\min\\{"; } else { texf="\\max\\{"; } for (j=1;j<=size(f);j++) { texf=texf+texPolynomial(f[j]); if (j=0) { "The embedded graph of the curve has not genus one."; } return(-1); } else { intmat nullmat[2][1]; // used to set // 4) find a vertex which has only one bounded edge, // if none exists zero is returned, // otherwise the number of the vertex in the list graph int nonloopvertex=findNonLoopVertex(graph); int dv; //checks if vert. has been found to which nonloopvertex is connected intmat delvert; // takes for a moment graph[i][3] of the vertex // to which nonloopvertex is connected // 5) delete successively vertices in the graph which // have only one bounded edge while (nonloopvertex>0) { // find the only vertex to which the nonloopvertex // is connected, when it is found // delete the connection in graph[i][3] and set dv=1 dv=0; for (i=1;i<=size(graph)-1 and dv==0;i++) { if (i!=nonloopvertex) { for (j=1;j<=ncols(graph[i][3]) and dv==0;j++) { if(graph[i][3][1,j]==nonloopvertex) // the vertex is found { delvert=graph[i][3]; delvert=intmatcoldelete(delvert,j); // delete the connection (note // there must have been two!) dv=1; graph[i][3]=delvert; } } } } graph[nonloopvertex][3]=nullmat; // the only connection of nonloopvertex // is killed nonloopvertex=findNonLoopVertex(graph); // find the next vertex // which has only one edge } // 6) find the loop and the weights of the edges intvec loop,weights; // encodes the loop and the edges i=1; // start by finding some vertex which belongs to the loop while (loop==0) { // if graph[i][3] of a vertex in the loop has 2 columns, all others have 1 if (ncols(graph[i][3])==1) { i++; } else { loop[1]=i; // a starting vertex is found loop[2]=graph[i][3][1,1]; // it is connected to vertex with this number weights[2]=graph[i][3][2,1]; // and the edge has this weight } } j=graph[i][3][1,1]; // the active vertex of the loop k=2; // counts the vertices in the loop while (j!=i) // the loop ends with the same vertex with which it starts { // the first row of graph[j][3] has two entries // corresponding to the two vertices // to which the active vertex j is connected; // one is loop[k-1], i.e. the one which // precedes j in the loop; we have to choose the other one if (graph[j][3][1,1]==loop[k-1]) { loop[k+1]=graph[j][3][1,2]; weights[k+1]=graph[j][3][2,2]; } else { loop[k+1]=graph[j][3][1,1]; weights[k+1]=graph[j][3][2,1]; } j=loop[k+1]; // set loop[k+1] the new active vertex k++; } // 7) compute for each edge in the loop the lattice length poly xcomp,ycomp; // the x- and y-components of the vectors // connecting two vertices of the loop number nenner; // the product of the denominators of // the x- and y-components number jinvariant; // the j-invariant int eins,zwei,ggt; for (i=1;i<=size(loop)-1;i++) // compute the lattice length for each edge { xcomp=graph[loop[i]][1]-graph[loop[i+1]][1]; ycomp=graph[loop[i]][2]-graph[loop[i+1]][2]; nenner=denominator(leadcoef(xcomp))*denominator(leadcoef(ycomp)); execute("eins="+string(numerator(leadcoef(nenner*xcomp)))+";"); execute("zwei="+string(numerator(leadcoef(nenner*ycomp)))+";"); ggt=gcd(eins,zwei); // the lattice length is the "gcd" // of the x-component and the y-component jinvariant=jinvariant+ggt/(nenner*weights[i+1]); // divided by the // weight of the edge } return(jinvariant); } } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),dp; // tropcialJInvariant computes the tropical j-invariant of an elliptic curve tropicalJInvariant(t*(x3+y3+1)+1/t*(x2+y2+x+y+x2y+xy2)+1/t2*xy); // the Newton polygone need not be the standard simplex tropicalJInvariant(x+y+x2y+xy2+1/t*xy); // the curve can have arbitrary degree tropicalJInvariant(t*(x7+y7+1)+1/t*(x4+y4+x2+y2+x3y+xy3)+1/t7*x2y2); // the procedure does not realise, if the embedded graph of the tropical // curve has a loop that can be resolved tropicalJInvariant(1+x+y+xy+tx2y+txy2); // but it does realise, if the curve has no loop at all ... tropicalJInvariant(x+y+1); // or if the embedded graph has more than one loop - even if only one // cannot be resolved tropicalJInvariant(1+x+y+xy+tx2y+txy2+t3x5+t3y5+tx2y2+t2xy4+t2yx4); } ///////////////////////////////////////////////////////////////////////// proc weierstrassForm (poly f,list #) "USAGE: weierstrassForm(wf[,#]); wf poly, # list ASSUME: wf is a a polynomial whose Newton polygon has precisely one interior lattice point, so that it defines an elliptic curve on the toric surface corresponding to the Newton polygon RETURN: poly, the Weierstrass normal form of the polynomial NOTE: - the algorithm for the coefficients of the Weierstrass form is due to Fernando Rodriguez Villegas, villegas@math.utexas.edu @* - the characteristic of the base field should not be 2 or 3 @* - if an additional argument # is given, a simplified Weierstrass form is computed EXAMPLE: example weierstrassForm; shows an example" { // Check first if the Newton polygon has precisely one interior point // - compute the Newton polygon list polygon=newtonPolytopeLP(f); // - find the vertices of the Newton cycle and order it clockwise list pg=findOrientedBoundary(polygon)[2]; // - check if there is precisely one interior point in the Newton polygon if (picksFormula(pg)[3]!=1) { ERROR("The Newton polygon of the input has NOT precisely one interior point!"); } // Compute the Normal form of the polygon. list nf=ellipticNF(polygon); // Transform f by the unimodular affine coordinate transformation poly tf=coordinatechange(f,nf[2],nf[3]); // Check if the Newton polygon is of type 4x2, 2x2 or a cubic poly a,b=flatten(coeffs(tf,ideal(var(1)^4,var(1)^2*var(2)^2))); // if it is of type 4x2 if (a!=0) { return(weierstrassFormOfA4x2Curve(tf)); } // if it is of type 2x2 if (b!=0) { return(weierstrassFormOfA2x2Curve(tf)); } // else it is a cubic return(weierstrassFormOfACubic(tf,#)); } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),lp; // f is already in Weierstrass form poly f=y2+yx+3y-x3-2x2-4x-6; weierstrassForm(f); // g is not, but wg is poly g=x+y+x2y+xy2+1/t*xy; poly wg=weierstrassForm(g); wg; // but it is not yet simple, since it still has an xy-term, unlike swg poly swg=weierstrassForm(g,1); swg; // the j-invariants of all three polynomials coincide jInvariant(g); jInvariant(wg); jInvariant(swg); // the following curve is elliptic as well poly h=x22y11+x19y10+x17y9+x16y9+x12y7+x9y6+x7y5+x2y3; // its Weierstrass form is weierstrassForm(h); } ///////////////////////////////////////////////////////////////////////// proc jInvariant (poly f,list #) "USAGE: jInvariant(f[,#]); f poly, # list ASSUME: - f is a a polynomial whose Newton polygon has precisely one interior lattice point, so that it defines an elliptic curve on the toric surface corresponding to the Newton polygon @* - it the optional argument # is present the base field should be Q(t) and the optional argument should be one of the following strings: @* 'ord' : then the return value is of type integer, namely the order of the j-invariant @* 'split' : then the return value is a list of two polynomials, such that the quotient of these two is the j-invariant RETURN: poly, the j-invariant of the elliptic curve defined by poly NOTE: the characteristic of the base field should not be 2 or 3, unless the input is a plane cubic EXAMPLE: example jInvariant; shows an example" { // compute first the Weierstrass form of the cubic and then the j-invariant if (size(#)==0) { return(jInvariantOfACubic(weierstrassForm(f),#)); } else { if (#[1]=="split") { return(jInvariantOfAPuiseuxCubic(weierstrassForm(f))); } else { return(jInvariantOfAPuiseuxCubic(weierstrassForm(f),"ord")); } } } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),dp; // jInvariant computes the j-invariant of a cubic jInvariant(x+y+x2y+y3+1/t*xy); // if the ground field has one parameter t, then we can instead // compute the order of the j-invariant jInvariant(x+y+x2y+y3+1/t*xy,"ord"); // one can compare the order of the j-invariant to the tropical j-invariant tropicalJInvariant(x+y+x2y+y3+1/t*xy); // the following curve is elliptic as well poly h=x22y11+x19y10+x17y9+x16y9+x12y7+x9y6+x7y5+x2y3+x14y8; // its j-invariant is jInvariant(h); } /////////////////////////////////////////////////////////////////////////////// /// Procedures concerned with conics /////////////////////////////////////////////////////////////////////////////// proc conicWithTangents (list points,list #) "USAGE: conicWithTangents(points[,#]); points list, # optional list ASSUME: points is a list of five points in the plane over K(t) RETURN: list, l[1] = the list points of the five given points @* l[2] = the conic f passing through the five points @* l[3] = list of equations of tangents to f in the given points @* l[4] = ideal, tropicalisation of f (i.e. list of linear forms) @* l[5] = a list of the tropicalisation of the tangents @* l[6] = a list containing the vertices of the tropical conic f @* l[7] = a list containing lists with vertices of the tangents @* l[8] = a string which contains the latex-code to draw the tropical conic and its tropicalised tangents @* l[9] = if # is non-empty, this is the same data for the dual conic and the points dual to the computed tangents NOTE: the points must be generic, i.e. no three on a line EXAMPLE: example conicWithTangents; shows an example" { int i; // Compute the tropical coordinates of the given points - part 1 list pointdenom,pointnum; poly pp1,pp2; for (i=1;i<=size(points);i++) { pp1=denominator(leadcoef(points[i][1])); pp2=denominator(leadcoef(points[i][2])); pointdenom[i]=list(pp1,pp2); pp1=numerator(leadcoef(points[i][1])); pp2=numerator(leadcoef(points[i][2])); pointnum[i]=list(pp1,pp2); } // compute the equation of the conic def BASERING=basering; ring LINRING=(0,t),(x,y,a(1..6)),lp; list points=imap(BASERING,points); ideal I; // the ideal will contain the linear equations given by the conic // and the points for (i=1;i<=5;i++) { I=I,substitute(a(1)*x2+a(2)*xy+a(3)*y2+a(4)*x+a(5)*y+a(6),x,points[i][1],y,points[i][2]); } I=std(I); list COEFFICIENTS; ideal DENOMINATORS; for (i=1;i<=6;i++) { COEFFICIENTS[i]=leadcoef(reduce(a(i),I)); DENOMINATORS[i]=denominator(leadcoef(COEFFICIENTS[i])); } // compute the common denominator of the coefficients ring TRING=0,t,dp; ideal DENOMINATORS=imap(LINRING,DENOMINATORS); poly p=lcm(DENOMINATORS); // compute the tropical coordinates of the points - part 2 ring tRING=0,t,ls; list pointdenom=imap(BASERING,pointdenom); list pointnum=imap(BASERING,pointnum); intvec pointcoordinates; for (i=1;i<=size(pointdenom);i++) { pointcoordinates[2*i-1]=ord(pointnum[i][1])-ord(pointdenom[i][1]); pointcoordinates[2*i]=ord(pointnum[i][2])-ord(pointdenom[i][2]); } // multiply the coefficients of the conic by the common denominator setring LINRING; poly p=imap(TRING,p); for (i=1;i<=6;i++) { COEFFICIENTS[i]=COEFFICIENTS[i]*p; } // Compute the tangents to the conic in the given points // and tropicalise the conic and the tangents setring BASERING; list CO=imap(LINRING,COEFFICIENTS); poly f=CO[1]*x2+CO[2]*xy+CO[3]*y2+CO[4]*x+CO[5]*y+CO[6]; poly fx=diff(f,x); poly fy=diff(f,y); list tangents; list tropicaltangents; for (i=1;i<=5;i++) { tangents[i]=substitute(fx,x,points[i][1],y,points[i][2])*x+substitute(fy,x,points[i][1],y,points[i][2])*y; tangents[i]=tangents[i]-substitute(tangents[i],x,points[i][1],y,points[i][2]); tropicaltangents[i]=tropicalise(tangents[i]); } list tropicalf=tropicalise(f); // compute the vertices of the tropcial conic and of the tangents list graphf=tropicalCurve(tropicalf); list eckpunktef; for (i=1;i<=size(graphf)-1;i++) { eckpunktef[i]=list(graphf[i][1],graphf[i][2]); } list eckpunktetangents; for (i=1;i<=size(tropicaltangents);i++) { eckpunktetangents[i]=verticesTropicalCurve(tropicaltangents[i],0); } // find the minimal and maximal coordinates of vertices poly minx,miny,maxx,maxy=graphf[1][1],graphf[1][2],graphf[1][1],graphf[1][2]; for (i=2;i<=size(graphf)-1;i++) { minx=minOfPolys(list(minx,graphf[i][1])); miny=minOfPolys(list(miny,graphf[i][2])); maxx=-minOfPolys(list(-maxx,-graphf[i][1])); maxy=-minOfPolys(list(-maxy,-graphf[i][2])); } // find the scale factor for the texdraw image poly maxdiffx,maxdiffy; maxdiffx=maxx-minx; maxdiffy=maxy-miny; if (maxdiffx==0) { maxdiffx=1; } if (maxdiffy==0) { maxdiffy=1; } poly scalefactor=minOfPolys(list(12/leadcoef(maxdiffx),16/leadcoef(maxdiffy)))/4; // Produce a Latex-Image of the Conic with each of its Tangents string zusatzinfo; string TEXBILD="\\documentclass[12pt]{amsart} \\usepackage{texdraw} \\begin{document} \\parindent0cm \\begin{center} \\large\\bf A Tropical Conic through 5 Points and Tangents \\end{center} We consider the concic through the following five points: \\begin{displaymath} "; string texf=texDrawTropical(graphf,list("",scalefactor)); for (i=1;i<=size(points);i++) { TEXBILD=TEXBILD+"p_"+string(i)+"=("+texNumber(points[i][1])+","+texNumber(points[i][2])+"),\\;\\;\\;"; } TEXBILD=TEXBILD+" \\end{displaymath}"; for (i=1;i<=size(points);i++) { zusatzinfo=" \\move ("+decimal(pointcoordinates[2*i-1])+" "+decimal(pointcoordinates[2*i])+") \\fcir f:0.5 r:0.25 \\rmove (0.4 0.4) \\htext{$p_"+string(i)+"$=("+string(pointcoordinates[2*i-1])+","+string(pointcoordinates[2*i])+")} "; TEXBILD=TEXBILD+" Here we draw the conic and the tangent through point $p_"+string(i)+"$: \\vspace*{1cm} \\begin{center} "+texDrawBasic(list(texf,zusatzinfo,texDrawTropical(tropicalCurve(tropicaltangents[i]),list("",scalefactor," \\setgray 0.8 \\linewd 0.1 \\lpatt (0.1 0.4)"))))+ "\\end{center}"; if (i0) { list dualpoints; for (i=1;i<=size(points);i++) { 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)); } list dualimage=conicWithTangents(dualpoints); return(list(points,f,tangents,tropicalf,tropicaltangents,eckpunktef,eckpunktetangents,TEXBILD,dualimage)); } else { return(list(points,f,tangents,tropicalf,tropicaltangents,eckpunktef,eckpunktetangents,TEXBILD)); } } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),dp; // the input consists of a list of five points in the plane over Q(t) list points=list(1/t2,t),list(1/t,t2),list(1,1),list(t,1/t2),list(t2,1/t); list conic=conicWithTangents(points); // conic[1] is the list of the given five points conic[1]; // conic[2] is the equation of the conic f passing through the five points conic[2]; // conic[3] is a list containing the equations of the tangents // through the five points conic[3]; // conic[4] is an ideal representing the tropicalisation of the conic f conic[4]; // conic[5] is a list containing the tropicalisation // of the five tangents in conic[3] conic[5]; // conic[6] is a list containing the vertices of the tropical conic conic[6]; // conic[7] is a list containing the vertices of the five tangents conic[7]; // conic[8] contains the latex code to draw the tropical conic and // its tropicalised tangents; it can written in a file, processed and // displayed via kghostview write(":w /tmp/conic.tex",conic[8]); system("sh","cd /tmp; latex /tmp/conic.tex; dvips /tmp/conic.dvi -o; kghostview conic.ps &"); // with an optional argument the same information for the dual conic is computed // and saved in conic[9] conic=conicWithTangents(points,1); conic[9][2]; // the equation of the dual conic } /////////////////////////////////////////////////////////////////////////////// /// Procedures concerned with tropicalisation /////////////////////////////////////////////////////////////////////////////// proc tropicalise (poly f,list #) "USAGE: tropicalise(f[,#]); f polynomial, # optional list ASSUME: f is a polynomial in Q(t)[x_1,...,x_n] RETURN: list, the linear forms of the tropicalisation of f NOTE: if # is empty, then the valuation of t will be 1, @* if # is the string 'max' it will be -1; @* the latter supposes that we consider the maximum of the computed linear forms, the former that we consider their minimum EXAMPLE: example tropicalise; shows an example" { int order,j; list tropicalf; intvec exp; int i=0; while (f!=0) { i++; exp=leadexp(f); if (size(#)==0) { tropicalf[i]=simplifyToOrder(f)[1]; } else { tropicalf[i]=-simplifyToOrder(f)[1]; } for (j=1;j<=nvars(basering);j++) { tropicalf[i]=tropicalf[i]+exp[j]*var(j); } f=f-lead(f); } return(tropicalf); } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),dp; tropicalise(2t3x2-1/t*xy+2t3y2+(3t3-t)*x+ty+(t6+1)); } ///////////////////////////////////////////////////////////////////////// proc tropicaliseSet (ideal i) "USAGE: tropicaliseSet(i); i ideal ASSUME: i is an ideal in Q(t)[x_1,...,x_n] RETURN: list, the jth entry is the tropicalisation of the jth generator of i EXAMPLE: example tropicaliseSet; shows an example" { list tropicalid; for (int j=1;j<=size(i);j++) { tropicalid[j]=tropicalise(i[j]); } return(tropicalid); } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),dp; ideal i=txy-y2+1,2t3x2+1/t*y-t6; tropicaliseSet(i); } ///////////////////////////////////////////////////////////////////////// proc tInitialForm (poly f, intvec w) "USAGE: tInitialForm(f,w); f a polynomial, w an integer vector ASSUME: f is a polynomial in Q[t,x_1,...,x_n] and w=(w_0,w_1,...,w_n) RETURN: poly, the t-initialform of f(t,x) w.r.t. w evaluated at t=1 NOTE: the t-initialform is the sum of the terms with MAXIMAL weighted order w.r.t. w EXAMPLE: example tInitialForm; shows an example" { // take in lead(f) only the term of lowest t-order and set t=1 poly initialf=lead(f); // compute the order of lead(f) w.r.t. (1,w) int gewicht=scalarproduct(w,leadexp(f)); // do the same for the remaining part of f and compare the results // keep only the smallest ones int vglgewicht; f=f-lead(f); while (f!=0) { vglgewicht=scalarproduct(w,leadexp(f)); if (vglgewicht>gewicht) { initialf=lead(f); gewicht=vglgewicht; } else { if (vglgewicht==gewicht) { initialf=initialf+lead(f); } } f=f-lead(f); } initialf=substitute(initialf,var(1),1); return(initialf); } example { "EXAMPLE:"; echo=2; ring r=0,(t,x,y),dp; poly f=t4x2+y2-t2xy+t4x-t9; intvec w=-1,-2,-3; tInitialForm(f,w); } ///////////////////////////////////////////////////////////////////////// proc tInitialIdeal (ideal i,intvec w,list #) "USAGE: tInitialIdeal(i,w); i ideal, w intvec ASSUME: i is an ideal in Q[t,x_1,...,x_n] and w=(w_0,...,w_n) RETURN: ideal ini, the t-initial ideal of i with respect to w" { // THE PROCEDURE WILL BE CALLED FROM OTHER PROCEDURES INSIDE THIS LIBRARY; // IN THIS CASE THE VARIABLE t WILL INDEED BE THE LAST VARIABLE INSTEAD OF // THE FIRST, // AND WE THEREFORE HAVE TO MOVE IT BACK TO THE FRONT! // THIS IS NOT DOCUMENTED FOR THE GENERAL USER!!!! def BASERING=basering; int j; if (size(#)>0) { // we first have to move the variable t to the front again ideal variablen=var(nvars(basering)); for (j=1;j=1;j--) { if (w[j]<>0) // find the last non-zero component of w { int r=2; for (int k=1;k<=size(w);k++) { // fill the order matrix O with unit vectors if(k<>j) // except the unit vector of the non-zero component { intvec u; u[size(w)]=0; u[k]=1; O[r,1..size(w)]=u; r=r+1; } } return(O); } } } example { "EXAMPLE:"; echo=2; intvec v=2,3,0,0; weightVectorToOrderMatrix(v); } ///////////////////////////////////////////////////////////////////////// proc initialIdeal (ideal i, intvec w) "USAGE: initialIdeal(i,w); i ideal, w intvec ASSUME: i is an ideal in Q[x_1,...,x_n] and w=(w_1,...,w_n) RETURN: ideal, the initialIdeal of i w.r.t. w NOTE: the initialIdeal consists of the terms with MAXIMAL weighted order w.r.t. w EXAMPLE: example initialIdeal; shows an example" { def BASERING=basering; intmat O=weightVectorToOrderMatrix(w); execute("ring INITIALRING=("+charstr(BASERING)+"),("+varstr(basering)+"),M("+string(O)+");"); ideal i=imap(BASERING,i); i=std(i); ideal ini; for (int j=1;j<=size(i);j++) { ini[j]=initialForm(i[j],w); } setring BASERING; return(imap(INITIALRING,ini)); } example { "EXAMPLE:"; echo=2; ring r=0,(x,y),dp; poly f=x3+y2-xy+x-1; intvec w=2,3; initialIdeal(f,w); } /////////////////////////////////////////////////////////////////////////////// /// PROCEDURES CONCERNED WITH THE LATEX CONVERSION /////////////////////////////////////////////////////////////////////////////// proc texNumber (poly p) "USAGE: texNumber(f); f poly RETURN: string, tex command representing leading coefficient of f using \frac EXAMPLE: example texNumber; shows an example" { number n=leadcoef(p); number den=denominator(n); number num=numerator(n); if (parstr(basering)=="") { if (den==-1) { den=-den; num=-num; } if (den==1) { return(string(num)); } else { return("\\tfrac{"+string(num)+"}{"+string(den)+"}"); } } def BASERING=basering; execute("ring PARAMETERRING=("+string(char(basering))+"),("+parstr(basering)+"),ds;"); poly den=imap(BASERING,den); poly num=imap(BASERING,num); if (den==1) { return(texPolynomial(num)); } else { return("\\tfrac{"+texPolynomial(num)+"}{"+texPolynomial(den)+"}"); } } example { "EXAMPLE:"; echo=2; ring r=(0,t),x,dp; texNumber((3t2-1)/t3); } ///////////////////////////////////////////////////////////////////////// proc texPolynomial (poly f) "USAGE: texPolynomial(f); f poly RETURN: string, the tex command representing f EXAMPLE: example texPolynomial; shows an example" { int altshort=short; short=0; if (f==0) { return("0"); } string texmf; string texco; string texf; string plus; int laenge=size(f); while (f!=0) { if (size(f)=1) { if (typeof(#[size(#)])=="string") { if (#[size(#)]!="max") { texdrawtp=texdrawtp+#[size(#)]; } } } // find the minimal and maximal coordinates of vertices // and find the scale factor for the texdraw image list SCFINPUT; SCFINPUT[1]=graph; list SCF=minScaleFactor(SCFINPUT); poly minx,miny,maxx,maxy=SCF[4],SCF[5],SCF[6],SCF[7]; poly centerx,centery=SCF[8],SCF[9]; int nachkomma=2; // number of decimals for the scalefactor number sf=1; // correction factor for scalefactor // if no scale factor was handed over to the procedure, use the // one computed by minScaleFactor; // check first if a scale factor has been handed over i=1; int scfpresent; while ((i<=size(#)) and (scfpresent==0)) { // if the scalefactor as polynomial was handed over, get it if (typeof(#[i])=="poly") { poly scalefactor=#[2]; scfpresent=1; } // if the procedure is called for drawing more than one tropical curve // then scalefactor,sf,nachkomma,minx,miny,maxx,maxy,centerx,centery // has been handed over to the procedure if (typeof(#[i])=="list") { poly scalefactor=#[i][1]; sf=#[i][2]; nachkomma=#[i][3]; minx=#[i][4]; miny=#[i][5]; maxx=#[i][6]; maxy=#[i][7]; centerx=#[i][8]; centery=#[i][9]; scfpresent=1; } i++; } // if no scalefactor was handed over we take the one computed in SCF if (scfpresent==0) { poly scalefactor=SCF[1]; sf=SCF[2]; nachkomma=SCF[3]; texdrawtp=texdrawtp+" \\relunitscale "+ decimal(scalefactor,nachkomma); } // compute the texdrawoutput for the tropical curve given by f list relxy; for (i=1;i<=size(graph)-1;i++) { // if the curve is a bunch of lines no vertex has to be drawn if (bunchoflines==0) { texdrawtp=texdrawtp+" \\move ("+decimal((graph[i][1]-centerx)/sf)+" "+decimal((graph[i][2]-centery)/sf)+") \\fcir f:0 r:"+decimal(2/(leadcoef(scalefactor)*10),size(string(int(scalefactor)))+1); } // draw the bounded edges emerging from the ith vertex for (j=1;j<=ncols(graph[i][3]);j++) { // don't draw it twice - and if there is only one vertex // and graph[i][3][1,1] is thus 0, nothing is done if (i1) and (noweights==0)) { texdrawtp=texdrawtp+" \\htext ("+decimal((graph[i][1]-centerx+graph[graph[i][3][1,j]][1])/(2*sf))+" "+decimal((graph[i][2]-centery+graph[graph[i][3][1,j]][2])/(2*sf))+"){$"+string(graph[i][3][2,j])+"$}"; } } } // draw the unbounded edges emerging from the ith vertex // they should not be too long for (j=1;j<=size(graph[i][4]);j++) { relxy=shorten(list(decimal((3*graph[i][4][j][1][1]/scalefactor)*sf),decimal((3*graph[i][4][j][1][2]/scalefactor)*sf),string(5*sf/2))); texdrawtp=texdrawtp+" \\move ("+decimal((graph[i][1]-centerx)/sf)+" "+decimal((graph[i][2]-centery)/sf)+") \\rlvec ("+relxy[1]+" "+relxy[2]+")"; // if the multiplicity is more than one, denote it in the picture if ((graph[i][4][j][2]>1) and (noweights==0)) { texdrawtp=texdrawtp+" \\htext ("+decimal(graph[i][1]/sf-centerx/sf+graph[i][4][j][1][1]/scalefactor)+" "+decimal(graph[i][2]/sf-centery/sf+graph[i][4][j][1][2]/scalefactor)+"){$"+string(graph[i][4][j][2])+"$}"; } } } // add lattice points if the scalefactor is >= 1/2 // and was not handed over if ((scalefactor>1/2) and (scfpresent==0)) { int uh=1; if (scalefactor>3) { uh=0; } texdrawtp=texdrawtp+" %% HERE STARTS THE CODE FOR THE LATTICE"; for (i=int(minx)-uh;i<=int(maxx)+uh;i++) { for (j=int(miny)-uh;j<=int(maxy)+uh;j++) { texdrawtp=texdrawtp+" \\move ("+decimal(i-centerx)+" "+decimal(j-centery)+") \\fcir f:0.8 r:"+decimal(1/(10*scalefactor),size(string(int(scalefactor)))+1); } } texdrawtp=texdrawtp+" %% HERE ENDS THE CODE FOR THE LATTICE "; } return(texdrawtp); } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),dp; poly f=x+y+x2y+xy2+1/t*xy; list graph=tropicalCurve(f); // compute the texdraw code of the tropical curve defined by f texDrawTropical(graph); // compute the texdraw code again, but set the scalefactor to 1 texDrawTropical(graph,"",1); } ///////////////////////////////////////////////////////////////////////// proc texDrawNewtonSubdivision (list graph,list #) "USAGE: texDrawNewtonSubdivision(graph[,#]); graph list, # optional list ASSUME: graph is the output of tropicalCurve RETURN: string, the texdraw code of the Newton subdivision of the tropical plane curve encoded by graph NOTE: - the list # may contain optional arguments, of which only one will be considered, namely the first entry of type 'poly'; this entry should be a rational number which specifies the scaling factor to be used; if it is missing, the factor will be computed; the list # may as well be empty @* - note that lattice points in the Newton subdivision which are black correspond to markings of the marked subdivision, while lattice points in grey are not marked EXAMPLE: example texDrawNewtonSubdivision; shows an example" { int i,j,k,l; list boundary=graph[size(graph)][1]; list inneredges=graph[size(graph)][2]; intvec shiftvector=graph[size(graph)][3]; string subdivision; // find maximal and minimal x- and y-coordinates and define the scalefactor poly maxx,maxy=1,1; for (i=1;i<=size(boundary);i++) { maxx=-minOfPolys(list(-maxx,-boundary[i][1])); maxy=-minOfPolys(list(-maxy,-boundary[i][2])); } // compute the scalefactor poly scalefactor=minOfPolys(list(12/leadcoef(maxx),12/leadcoef(maxy))); // Check if the scalefactor was handed over, then take that instead if (size(#)>0) { int scfpresent; i=1; while (i<=size(#) and scfpresent==0) { if (typeof(#[i])=="poly") { scalefactor=#[i]; scfpresent=1; } i++; } } // check if scaling is necessary if (scalefactor<1) { subdivision=subdivision+" \\relunitscale"+ decimal(scalefactor); } // define the texdraw code for the Newton subdivision for (i=1;i<=size(boundary)-1;i++) { subdivision=subdivision+" \\move ("+string(boundary[i][1])+" "+string(boundary[i][2])+") \\lvec ("+string(boundary[i+1][1])+" "+string(boundary[i+1][2])+")"; } subdivision=subdivision+" \\move ("+string(boundary[size(boundary)][1])+" "+string(boundary[size(boundary)][2])+") \\lvec ("+string(boundary[1][1])+" "+string(boundary[1][2])+") "; for (i=1;i<=size(inneredges);i++) { subdivision=subdivision+" \\move ("+string(inneredges[i][1][1])+" "+string(inneredges[i][1][2])+") \\lvec ("+string(inneredges[i][2][1])+" "+string(inneredges[i][2][2])+")"; } // add lattice points if the scalefactor is >= 1/2 if (scalefactor>1/2) { for (i=shiftvector[1];i<=int(maxx);i++) { for (j=shiftvector[2];j<=int(maxy);j++) { if (scalefactor > 2) { subdivision=subdivision+" \\move ("+string(i)+" "+string(j)+") \\fcir f:0.6 r:"+decimal(2/(10*scalefactor),size(string(int(scalefactor)))+1); } else { subdivision=subdivision+" \\move ("+string(i)+" "+string(j)+") \\fcir f:0.6 r:"+decimal(2/(20*scalefactor),size(string(int(scalefactor)))+1); } } } if ((shiftvector[1]!=0) or (shiftvector[2]!=0)) { subdivision=subdivision+" \\htext (-0.3 0.1){{\\tiny $(0,0)$}}"; } } // deal with the pathological cases if (size(boundary)==1) // then the Newton polytope is a point { subdivision=subdivision+" \\move ("+string(boundary[1][1])+" "+string(boundary[1][2])+") \\fcir f:0 r:0.15"; } else { matrix M[2][size(boundary)]; for (i=1;i<=size(boundary);i++) { M[1,i]=boundary[i][1]; M[2,i]=boundary[i][2]; } if ((size(boundary)-size(syz(M)))==1) { for (i=1;i<=size(boundary);i++) { subdivision=subdivision+" \\move ("+string(boundary[i][1])+" "+string(boundary[i][2])+") \\fcir f:0 r:"+decimal(2/(8*scalefactor),size(string(int(scalefactor)))+1); } } } // find the marked points in the subdivision poly pg; list markings; k=1; for (i=1;i=2;i--) { j=i-1; while (j>=1) { if (markings[i]==markings[j]) { markings=delete(markings,i); j=0; } j--; } } for (i=1;i<=size(markings);i++) { if (scalefactor > 2) { subdivision=subdivision+" \\move ("+string(markings[i][1])+" "+string(markings[i][2])+") \\fcir f:0 r:"+decimal(2/(8*scalefactor),size(string(int(scalefactor)))+1); } else { subdivision=subdivision+" \\move ("+string(markings[i][1])+" "+string(markings[i][2])+") \\fcir f:0 r:"+decimal(2/(16*scalefactor),size(string(int(scalefactor)))+1); } } // enclose subdivision in the texdraw environment string texsubdivision=" \\begin{texdraw} \\drawdim cm \\relunitscale 1 \\linewd 0.05" +subdivision+" \\end{texdraw} "; return(texsubdivision); } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),dp; poly f=x+y+x2y+xy2+1/t*xy; list graph=tropicalCurve(f); // compute the texdraw code of the Newton subdivision of the tropical curve texDrawNewtonSubdivision(graph); } ///////////////////////////////////////////////////////////////////////// proc texDrawTriangulation (list triang,list polygon) "USAGE: texDrawTriangulation(triang,polygon); triang,polygon list ASSUME: polygon is a list of integer vectors describing the lattice points of a marked polygon; triang is a list of integer vectors describing a triangulation of the marked polygon in the sense that an integer vector of the form (i,j,k) describes the triangle formed by polygon[i], polygon[j] and polygon[k] RETURN: string, a texdraw code for the triangulation described by triang without the texdraw environment EXAMPLE: example texDrawTriangulation; shows an example" { // the header of the texdraw output string latex=" \\drawdim cm \\relunitscale 1.2 \\arrowheadtype t:V "; int i,j; // indices list pairs,markings; // stores edges of the triangulation, respecively // the marked points for each triangle store the edges and marked // points of the triangle for (i=1;i<=size(triang);i++) { pairs[3*i-2]=intvec(triang[i][1],triang[i][2]); pairs[3*i-1]=intvec(triang[i][2],triang[i][3]); pairs[3*i]=intvec(triang[i][3],triang[i][1]); markings[3*i-2]=triang[i][1]; markings[3*i-1]=triang[i][2]; markings[3*i]=triang[i][3]; } // delete redundant pairs which occur more than once for (i=size(pairs);i>=1;i--) { for (j=1;j=1;i--) { for (j=1;jgewicht) { execute("koef="+leitkoef[2]+";"); initialf=koef*leadmonom(f); gewicht=vglgewicht; } else { if (vglgewicht==gewicht) { execute("koef="+leitkoef[2]+";"); initialf=initialf+koef*leadmonom(f); } } f=f-lead(f); } return(initialf); } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),dp; poly f=t4x2+y2-t2xy+t4x-1/t6; intvec w=2,3; tInitialFormParMax(f,w); } ///////////////////////////////////////////////////////////////////////// proc solveTInitialFormPar (ideal i) "USAGE: solveTInitialFormPar(i); i ideal ASSUME: i is a zero-dimensional ideal in Q(t)[x_1,...,x_n] generated by the (1,w)-homogeneous elements for some integer vector w - i.e. by the (1,w)-initialforms of polynomials RETURN: none NOTE: the procedure just displays complex approximations of the solution set of i EXAMPLE: example solveTInitialFormPar; shows an example" { i=subst(i,t,1); string CHARAKTERISTIK=charstr(basering); CHARAKTERISTIK=CHARAKTERISTIK[1..size(CHARAKTERISTIK)-2]; def BASERING=basering; execute("ring INITIALRING=("+CHARAKTERISTIK+"),("+varstr(basering)+"),("+ordstr(basering)+");"); list l=solve(imap(BASERING,i)); l; setring BASERING; // return(fetch(SUBSTRING,l)); } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),dp; ideal i=t2x2+y2,x-t2; solveTInitialFormPar(i); } ///////////////////////////////////////////////////////////////////////// proc detropicalise (def p) "USAGE: detropicalise(f); f poly or f list ASSUME: if f is of type poly then t is a linear polynomial with an arbitrary constant term and positive integer coefficients as further coefficients; @* if f is of type list then f is a list of polynomials of the type just described in before RETURN: poly, the detropicalisation of f ignoring the constant parts NOTE: the output will be a monomial and the constant coefficient has been ignored EXAMPLE: example detropicalise; shows an example" { poly dtp=1; while (p!=0) { if (leadmonom(p)!=1) { dtp=dtp*leadmonom(p)^int(leadcoef(p)); } p=p-lead(p); } return(dtp); } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),dp; detropicalise(3x+4y-1); } ///////////////////////////////////////////////////////////////////////// proc tDetropicalise (def p) "USAGE: tDetropicalise(f); f poly or f list ASSUME: if f is of type poly then f is a linear polynomial with an integer constant term and positive integer coefficients as further coefficients; @* if f is of type list then it is a list of polynomials of the type just described in before RETURN: poly, the detropicalisation of f over the field Q(t) NOTE: the output will be a term where the coeffiecient is a Laurent monomial in the variable t EXAMPLE: example tDetropicalise; shows an example" { poly dtp; if (typeof(p)=="list") { int i; for (i=1;i<=size(p);i++) { dtp=dtp+tDetropicalise(p[i]); } } else { dtp=1; while (p!=0) { if (leadmonom(p)!=1) { dtp=dtp*leadmonom(p)^int(leadcoef(p)); } else { dtp=t^int(leadcoef(p))*dtp; } p=p-lead(p); } } return(dtp); } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),dp; tDetropicalise(3x+4y-1); } /////////////////////////////////////////////////////////////////////////////// /// Auxilary Procedures concerned with conics /////////////////////////////////////////////////////////////////////////////// proc dualConic (poly f) "USAGE: dualConic(f); f poly ASSUME: f is an affine conic in two variables x and y RETURN: poly, the equation of the dual conic EXAMPLE: example dualConic; shows an example" { if (deg(f)!=2) { "The polynomial is not a conic"; return(0); } matrix A=coeffs(f,ideal(var(1)^2,var(1)*var(2),var(2)^2,var(1),var(2),1)); matrix C[3][3]=A[1,1],A[2,1]/2,A[4,1]/2,A[2,1]/2,A[3,1],A[5,1]/2,A[4,1]/2,A[5,1]/2,A[6,1]; matrix M[3][1]=var(1),var(2),1; matrix Cdual=-adjoint(C); matrix ffd=transpose(M)*Cdual*M; poly fdual=ffd[1,1]; return(fdual); } example { "EXAMPLE:"; echo=2; ring r=0,(x,y),dp; poly conic=2x2+1/2y2-1; dualConic(conic); } /////////////////////////////////////////////////////////////////////////////// /// Auxilary Procedures concerned with substitution /////////////////////////////////////////////////////////////////////////////// proc parameterSubstitute (poly f,int N) "USAGE: parameterSubstitute(f,N); f poly, N int ASSUME: f is a polynomial in Q(t)[x_1,...,x_n] describing a plane curve over Q(t) RETURN: poly f with t replaced by t^N EXAMPLE: example parameterSubstitute; shows an example" { def BASERING=basering; number nenner=1; for (int i=1;i<=size(f);i++) { nenner=nenner*denominator(leadcoef(f[i])); } f=f*nenner; f=subst(f,t,t^N)/number(subst(nenner,t,t^N)); return(f); } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),dp; poly f=t2xy+1/t*y+t3; parameterSubstitute(f,3); parameterSubstitute(f,-1); } ///////////////////////////////////////////////////////////////////////// proc tropicalSubst (poly f,int N,list #) "USAGE: parameterSubstitute(f,N,L); f poly, N int, L list ASSUME: f is a polynomial in Q(t)[x_1,...,x_k] and L is a list of the form var(i_1),poly_1,...,v(i_k),poly_k RETURN: list, the list is the tropical polynomial which we get from f by replacing the i-th variable be the i-th polynomial but in the i-th polynomial the parameter t is replaced by t^1/N EXAMPLE: example tropicalSubst; shows an example" { f=parameterSubstitute(f,N); f=substitute(f,#); list tp=tropicalise(f); poly const; for (int i=1;i<=size(tp);i++) { const=substitute(tp[i],x,0,y,0); tp[i]=tp[i]-const+const/N; } return(tp); } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),dp; poly f=t2x+1/t*y-1; tropicalSubst(f,2,x,x+t,y,tx+y+t2); // The procedure can be used to study the effect of a transformation of // the form x -> x+t^b, with b a rational number, on the tropicalisation and // the j-invariant of a cubic over the Puiseux series. f=t7*y3+t3*y2+t*(x3+xy2+y+1)+xy; // - the j-invariant, and hence its valuation, // does not change under the transformation jInvariant(f,"ord"); // - b=3/2, then the cycle length of the tropical cubic equals -val(j-inv) list g32=tropicalSubst(f,2,x,x+t3,y,y); tropicalJInvariant(g32); // - b=1, then it is still true, but only just ... list g1=tropicalSubst(f,1,x,x+t,y,y); tropicalJInvariant(g1); // - b=2/3, as soon as b<1, the cycle length is strictly less than -val(j-inv) list g23=tropicalSubst(f,3,x,x+t2,y,y); tropicalJInvariant(g23); } ///////////////////////////////////////////////////////////////////////// proc randomPoly (int d,int ug, int og, list #) "USAGE: randomPoly(d,ug,og[,#]); d, ug, og int, # list ASSUME: the basering has a parameter t RETURN: poly, a polynomial of degree d where the coefficients are of the form t^j with j a random integer between ug and og NOTE: if an optional argument # is given, then the coefficients are instead either of the form t^j as above or they are zero, and this is chosen randomly EXAMPLE: example randomPoly; shows an example" { int i,j,k; def BASERING=basering; ring RANDOMRING=(0,t),(x,y,z),dp; ideal m=maxideal(d); int nnn=size(m); for (i=1;i<=nnn;i++) { m[i]=subst(m[i],z,1); } poly randomPolynomial; for (i=1;i<=nnn;i++) { if (size(#)!=0) { k=random(0,1); } if (k==0) { j=random(ug,og); randomPolynomial=randomPolynomial+t^j*m[i]; } } setring BASERING; poly randomPolynomial=imap(RANDOMRING,randomPolynomial); if ((size(randomPolynomial)<3) and (nnn>=3)) { randomPolynomial=randomPoly(d,ug,og,#); } return(randomPolynomial); } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),dp; randomPoly(3,-2,5); randomPoly(3,-2,5,1); } ///////////////////////////////////////////////////////////////////////// proc cleanTmp () "USAGE: cleanTmp() PURPOSE: some procedures create latex and ps-files in the directory /tmp; in order to remove them simply call cleanTmp(); RETURN: none" { system("sh","cd /tmp; /bin/rm -f tropicalcurve*; /bin/rm -f tropicalnewtonsubdivision*"); } ////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////// /// AUXILARY PROCEDURES, WHICH ARE DECLARED STATIC ////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////// /// Procedures used in tropicalparametriseNoabs respectively in tropicalLifting: /// - phiOmega /// - cutdown /// - tropicalparametriseNoabs /// - findzeros /// - basictransformideal /// - testw /// - simplifyToOrder /// - scalarproduct /// - intvecdelete /// - invertorder /// - ordermaximalidealsNoabs /// - displaypoly /// - displaycoef /// - choosegfanvector /// - tropicalliftingresubstitute ////////////////////////////////////////////////////////////////////////////// static proc phiOmega (ideal i,int j,int wj) "USAGE: phiOmega(i,j,wj); i ideal, j int, wj int ASSUME: i an ideal in Q[t,x_1,...,x_n] and 00, in Q(t)[x]; the optional parameter # can contain the string 'isPrime' to indicate that the input ideal is prime and no minimal associated primes have to be computed RETURN: list, the first entry is a ring, namely the basering where some variables have been eliminated, and the ring contains the ideal i (with the same variables eliminated), the t-initial ideal ini of i (w.r.t. the weight vector where the entries corresponding to the deleted variables have been eliminated) and a list repl where for each eliminated variable there is one entry, namely a polynomial in the remaining variables and t that explains how resubstitution of a solution for the new i gives a solution for the old i; the second entry is the weight vector wvec with the components corresponding to the eliminated variables removed NOTE: needs the libraries random.lib and primdec.lib; is called from tropicalLifting" { // IDEA: i is an ideal of dimension d; we want to cut it with d random linear // forms in such a way that the resulting // ideal is 0-dim and still contains w in the tropical variety // NOTE: t is the last variable in the basering ideal pideal; //this is the ideal we want to return ideal cutideal; ideal pini; //this is the initial ideal ideal primini; //the initial of the prime we picl int j1,j2,dimp,winprim,j3; poly product=1; def BASERING=basering; ideal variablen; intvec setvec;intvec wvecp; setvec[nvars(basering)-1]=0; poly randomp=0; poly randomp1=0; poly randomp2=0; list erglini; list ergl; for (j1=1;j1<=nvars(basering)-1;j1++) { variablen=variablen+var(j1); // read the set of variables // (needed to make the quotring later) product=product*var(j1); // make product of all variables // (needed for the initial-monomial-check later } execute("ring QUOTRING=("+charstr(basering)+",t),("+string(variablen)+"),dp;"); setring BASERING; // change to quotring where we want to compute the primary decomposition of i if (size(#)==0) // we only have to do so if isPrime is not set { setring QUOTRING; ideal jideal=imap(BASERING,jideal); list primp=minAssGTZ(jideal); //compute the primary decomposition for (j1=1;j1<=size(primp);j1++) { for(j2=1;j2<=size(primp[j1]);j2++) { // clear all denominators primp[j1][j2]=primp[j1][j2]/content(primp[j1][j2]); } } setring BASERING; list primp=imap(QUOTRING,primp); // if i is not primary itself // go through the list of min. ass. primes and find the first // one which has w in its tropical variety if (size(primp)>1) { j1=1; while(winprim==0) { //compute the t-initial of the associated prime // - the last entry 1 only means that t is the last variable in the ring primini=tInitialIdeal(primp[j1],wvec,1); // check if it contains a monomial (resp if the product of var // is in the radical) if (radicalMemberShip(product,primini)==0) { // if w is in the tropical variety of the prime, we take that jideal=primp[j1]; setring QUOTRING; dimension=dim(groebner(imap(BASERING,jideal)));//compute its dimension setring BASERING; winprim=1; // and stop the checking } j1=j1+1; //else we look at the next associated prime } } else { jideal=primp[1]; //if i is primary itself we take its prime instead } } // now we start as a first try to intersect with a hyperplane parallel to // coordinate axes, because this would make our further computations // a lot easier. // We choose a subset of our n variables of size d=dim(ideal). // For each of these // variables, we want to fix a value: x_i= a_i*t^-w_i. // This will only work if the // projection of the d-dim variety to the other n-d variables // is the whole n-d plane. // Then a general choice for a_i will intersect the variety // in finitely many points. // If the projection is not the whole n-d plane, // then a general choice will not work. // We could determine if we picked a good // d-subset of variables using elimination // (NOTE, there EXIST d variables such that // a random choice of a_i's would work!). // But since this involves many computations, // we prefer to choose randomly and just // try in the end if our intersected ideal // satisfies our requirements. If this does not // work, we give up this try and use our second intersection idea, which // will work for a Zariksi-open subset (i.e. almost always). // // As random subset of d variables we choose // those for which the absolute value of the // wvec-coordinate is smallest, because this will // give us the smallest powers of t and hence // less effort in following computations. // Note that the smallest absolute value have those // which are biggest, because wvec is negative. //print("first try"); intvec wminust=intvecdelete(wvec,1); intmat A[2][size(wminust)]; // make a matrix with first row -wvec (without t) and second row 1..n A[1,1..size(wminust)]=-wminust; A[2,1..size(wminust)]=1..size(wminust); // sort this matrix in order to get // the d biggest entries and their position in wvec A=sortintmat(A); // we construct a vector which has 1 at entry j if j belongs to the list // of the d biggest entries of wvec and a 0 else for (j1=1;j1<=nvars(basering)-1;j1++) //go through the variables { for (j2=1;j2<=dimension;j2++) //go through the list of smallest entries { if (A[2,j2]==j1)//if the variable belongs to this list { setvec[j1]=1;//put a 1 } } } // using this 0/1-vector we produce // a random constant (i.e. coeff in Q times something in t) // for each of the biggest variables, // we add the forms x_i-random constant to the ideal // and we save the constant at the i-th place of // a list we want to return for later computations j3=0; while (j3<=1) { j3++; pideal=jideal; j2=1; for (j1=1;j1<=nvars(basering)-1;j1++) { if(setvec[j1]==1)//if x_i belongs to the biggest variables { if ((j3==1) and ((char(basering)==0) or (char(basering)>3))) { randomp1=random(1,3); randomp=t^(A[1,j2])*randomp1;// make a random constant // --- first we try small numbers } if ((j3==2) and ((char(basering)==0) or (char(basering)>100))) { randomp1=random(1,100); randomp=t^(A[1,j2])*randomp1;// make a random constant // --- next we try bigger numbers } else { randomp1=random(1,char(basering)-1); randomp=t^(A[1,j2])*randomp1;//make a random constant } ergl[A[2,j2]]=randomp;//and save the constant in a list erglini[A[2,j2]]=randomp1; j2=j2+1; } else { ergl[j1]=0; //if the variable is not among the d biggest ones, //save 0 in the list erglini[j1]=0; } } // print(ergl);print(pideal); // now we check if we made a good choice of pideal, i.e. if dim=0 and // wvec is still in the tropical variety // change to quotring where we compute dimension cutideal=pideal; for(j1=1;j1<=nvars(basering)-1;j1++) { if(setvec[j1]==1) { cutideal=cutideal,var(j1)-ergl[j1];//add all forms to the ideal } } setring QUOTRING; ideal cutideal=imap(BASERING,cutideal); dimp=dim(groebner(cutideal)); //compute the dim of p in the quotring //print("dimension"); //print(dimp); kill cutideal; setring BASERING; if (dimp==0) // if it is 0 as we want { // compute the t-initial of the associated prime // - the last 1 just means that the variable t is // the last variable in the ring pini=tInitialIdeal(cutideal,wvec ,1); //print("initial"); //print(pini); // and if the initial w.r.t. t contains no monomial // as we want (checked with // radical-membership of the product of all variables) if (radicalMemberShip(product,pini)==0) { // we made the right choice and now // we substitute the variables in the ideal // to get an ideal in less variables // also we make a projected vector // from wvec only the components of the remaining variables wvecp=wvec; variablen=0; j2=0; for(j1=1;j1<=nvars(basering)-1;j1++) { if(setvec[j1]==1)//if j belongs to the biggest variables { j2=j2+1; pideal=subst(pideal,var(j1),ergl[j1]);//substitute this variable pini=subst(pini,var(j1),erglini[j1]); wvecp=intvecdelete(wvecp,j1+2-j2); } else { variablen=variablen+var(j1); // read the set of remaining variables // (needed to make quotring later) } } // return pideal, the initial and the list ergl which tells us // which variables we replaced by which form execute("ring BASERINGLESS1=("+charstr(BASERING)+"),("+string(variablen)+",t),(dp("+string(ncols(variablen))+"),lp(1));"); ideal i=imap(BASERING,pideal); ideal ini=imap(BASERING,pini); //ideal ini2=tInitialIdeal(i,wvecp,1); list repl=imap(BASERING,ergl); export(i); export(ini); export(repl); return(list(BASERINGLESS1,wvecp)); } } } // this is our second try to cut down, which we only use if the first try // didn't work out. We intersect with d general hyperplanes // (i.e. we don't choose // them to be parallel to coordinate hyperplanes anymore. This works out with // probability 1. // // We choose general hyperplanes, i.e. linear forms which involve all x_i. // Each x_i has to be multiplied bz t^(w_i) in order // to get the same weight (namely 0) // for each term. As we cannot have negative exponents, we multiply // the whole form by t^minimumw. Notice that then in the first form, // there is one term without t- the term of the variable // x_i such that w_i is minimal. That is, we can solve for this variable. // In the second form, we can replace that variable, // and divide by t as much as possible. // Then there is again one term wihtout t - // the term of the variable with second least w. // So we can solve for this one again and also replace it in the first form. // Since all our coefficients are chosen randomly, // we can also from the beginning on // choose the set of variables which belong to the d smallest entries of wvec // (t not counting) and pick random forms g_i(t,x') // (where x' is the set of remaining variables) // and set x_i=g_i(t,x'). // // make a matrix with first row wvec (without t) and second row 1..n //print("second try"); setring BASERING; A[1,1..size(wminust)]=wminust; A[2,1..size(wminust)]=1..size(wminust); // sort this matrix in otder to get the d smallest entries // (without counting the t-entry) A=sortintmat(A); randomp=0; setvec=0; setvec[nvars(basering)-1]=0; // we construct a vector which has 1 at entry j if j belongs to the list of // the d smallest entries of wvec and a 0 else for (j1=1;j1<=nvars(basering)-1;j1++) //go through the variables { for (j2=1;j2<=dimension;j2++) //go through the list of smallest entries { if (A[2,j2]==j1)//if the variable belongs to this list { setvec[j1]=1;//put a 1 } } } variablen=0; wvecp=wvec; j2=0; for(j1=1;j1<=nvars(basering)-1;j1++) { if(setvec[j1]==1)//if j belongs to the biggest variables { j2=j2+1; wvecp=intvecdelete(wvecp,j1+2-j2);// delete the components // we substitute from wvec } else { variablen=variablen+var(j1); // read the set of remaining variables // (needed to make the quotring later) } } setring BASERING; execute("ring BASERINGLESS2=("+charstr(BASERING)+"),("+string(variablen)+",t),(dp("+string(ncols(variablen))+"),lp(1));"); // using the 0/1-vector which tells us which variables belong // to the set of smallest entries of wvec // we construct a set of d random linear // polynomials of the form x_i=g_i(t,x'), // where the set of all x_i is the set of // all variables which are in the list of smallest // entries in wvec, and x' are the other variables. // We add these d random linear polynomials to // the ideal pideal, i.e. we intersect // with these and hope to get something // 0-dim which still contains wvec in its // tropical variety. Also, we produce a list ergl // with g_i at the i-th position. // This is a list we want to return. // we try first with small numbers setring BASERING; pideal=jideal; for(j1=1;j1<=dimension;j1++)//go through the list of variables { // corres to the d smallest in wvec if ((char(basering)==0) or (char(basering)>3)) { randomp1=random(1,3); randomp=randomp1*t^(-A[1,j1]); } else { randomp1=random(1,char(basering)-1); randomp=randomp1*t^(-A[1,j1]); } for(j2=1;j2<=nvars(basering)-1;j2++)//go through all variables { if(setvec[j2]==0)//if x_j belongs to the set x' { // add a random term with the suitable power // of t to the random linear form if ((char(basering)==0) or (char(basering)>3)) { randomp2=random(1,3); randomp1=randomp1+randomp2*var(j2); randomp=randomp+randomp2*t^(wminust[j2]-A[1,j1])*var(j2); } else { randomp2=random(char(basering)-1); randomp1=randomp1+randomp2*var(j2); randomp=randomp+randomp2*t^(wminust[j2]-A[1,j1])*var(j2); } ergl[j2]=0;//if j belongs to x' we want a 0 in our list erglini[j2]=0; } } ergl[A[2,j1]]=randomp;// else we put there the random oly g_i(t,x') erglini[A[2,j1]]=randomp1; randomp=0; randomp1=0; } //print(ergl); // Again, we have to test if we made a good choice // to intersect,i.e. we have to check whether // pideal is 0-dim and contains wvec in the tropical variety. cutideal=pideal; for(j1=1;j1<=nvars(basering)-1;j1++) { if(setvec[j1]==1) { cutideal=cutideal,var(j1)-ergl[j1];//add all forms to the ideal } } setring QUOTRING; ideal cutideal=imap(BASERING,cutideal); dimp=dim(groebner(cutideal)); //compute the dim of p in the quotring //print("dimension"); //print(dimp); kill cutideal; setring BASERING; if (dimp==0) // if it is 0 as we want { // compute the t-initial of the associated prime // - the last 1 just means that the variable t // is the last variable in the ring pini=tInitialIdeal(cutideal,wvec ,1); //print("initial"); //print(pini); // and if the initial w.r.t. t contains no monomial as we want (checked with // radical-membership of the product of all variables) if (radicalMemberShip(product,pini)==0) { // we want to replace the variables x_i by the forms -g_i in // our ideal in order to return an ideal with less variables // first we substitute the chosen variables for(j1=1;j1<=nvars(basering)-1;j1++) { if(setvec[j1]==1)//if j belongs to the biggest variables { pideal=subst(pideal,var(j1),ergl[j1]);//substitute it pini=subst(pini,var(j1),erglini[j1]); } } setring BASERINGLESS2; ideal i=imap(BASERING,pideal); ideal ini=imap(BASERING,pini);//ideal ini2=tInitialIdeal(i,wvecp,1); list repl=imap(BASERING,ergl); export(i); export(ini); export(repl); return(list(BASERINGLESS2,wvecp)); } } // now we try bigger numbers while (1) //a never-ending loop which will stop with prob. 1 { // as we find a suitable ideal with that prob setring BASERING; pideal=jideal; for(j1=1;j1<=dimension;j1++)//go through the list of variables { // corres to the d smallest in wvec randomp1=random(1,100); randomp=randomp1*t^(-A[1,j1]); for(j2=1;j2<=nvars(basering)-1;j2++)//go through all variables { if(setvec[j2]==0)//if x_j belongs to the set x' { // add a random term with the suitable power // of t to the random linear form if ((char(basering)==0) or (char(basering)>100)) { randomp2=random(1,100); randomp1=randomp1+randomp2*var(j2); randomp=randomp+randomp2*t^(wminust[j2]-A[1,j1])*var(j2); } else { randomp2=random(1,char(basering)-1); randomp1=randomp1+randomp2; randomp=randomp+randomp2*t^(wminust[j2]-A[1,j1])*var(j2); } ergl[j2]=0;//if j belongs to x' we want a 0 in our list erglini[j2]=0; } } ergl[A[2,j1]]=randomp;// else we put there the random oly g_i(t,x') erglini[A[2,j1]]=randomp1; randomp=0; randomp1=0; } //print(ergl); // Again, we have to test if we made a good choice to // intersect,i.e. we have to check whether // pideal is 0-dim and contains wvec in the tropical variety. cutideal=pideal; for(j1=1;j1<=nvars(basering)-1;j1++) { if(setvec[j1]==1) { cutideal=cutideal,var(j1)-ergl[j1];//add all forms to the ideal } } setring QUOTRING; ideal cutideal=imap(BASERING,cutideal); dimp=dim(groebner(cutideal)); //compute the dim of p in the quotring //print("dimension"); //print(dimp); kill cutideal; setring BASERING; if (dimp==0) // if it is 0 as we want { // compute the t-initial of the associated prime // - the last 1 just means that the variable t // is the last variable in the ring pini=tInitialIdeal(cutideal,wvec ,1); //print("initial"); //print(pini); // and if the initial w.r.t. t contains no monomial // as we want (checked with // radical-membership of the product of all variables) if (radicalMemberShip(product,pini)==0) { // we want to replace the variables x_i by the forms -g_i in // our ideal in order to return an ideal with less variables //first we substitute the chosen variables for(j1=1;j1<=nvars(basering)-1;j1++) { if(setvec[j1]==1)//if j belongs to the biggest variables { pideal=subst(pideal,var(j1),ergl[j1]);//substitute it pini=subst(pini,var(j1),erglini[j1]); } } // return pideal and the list ergl which tells us // which variables we replaced by which form setring BASERINGLESS2; ideal i=imap(BASERING,pideal); ideal ini=imap(BASERING,pini);//ideal ini2=tInitialIdeal(i,wvecp,1); list repl=imap(BASERING,ergl); export(i); export(ini); export(repl); return(list(BASERINGLESS2,wvecp)); } } } } ///////////////////////////////////////////////////////////////////////// static proc tropicalparametriseNoabs (ideal i,intvec ww,int ordnung,int gfanold,int nogfan,int puiseux,list #) "USAGE: tropicalparametriseNoabs(i,tw,ord,gf,ng,pu[,#]); i ideal, tw intvec, ord int, gf,ng,pu int, # opt. list ASSUME: - i is an ideal in Q[t,x_1,...,x_n,X_1,...,X_k], tw=(w_0,w_1,...,w_n,0,...,0) and (w_0,...,w_n,0,...,0) is in the tropical variety of i, and ord is the order up to which a point in V(i) over C((t)) lying over w shall be computed; - moreover, k should be zero if the procedure is not called recursively; - the point in the tropical variety is supposed to lie in the NEGATIVE orthant; - the ideal is zero-dimensional when considered in (Q(t)[X_1,...,X_k]/m)[x_1,...,x_n], where m=#[2] is a maximal ideal in Q(t)[X_1,...,X_k]; - gf is 0 if version 2.2 or larger is used and it is 1 else - ng is 1 if gfan should not be executed - pu is 1 if n=1 and i is generated by one polynomial and newtonpoly is used to find a point in the tropical variety instead of gfan RETURN: list, l[1] = ring Q(0,X_1,...,X_r)[[t]] l[2] = int l[3] = string NOTE: - the procedure is also called recursively by itself, and if it is called in the first recursion, the list # is empty, otherwise #[1] is an integer, one more than the number of true variables x_1,...,x_n, and #[2] will contain the maximal ideal m in the variables X_1,...X_k by which the ring Q[t,x_1,...,x_n,X_1,...,X_k] should be divided to work correctly in K[t,x_1,...,x_n] where K=Q[X_1,...,X_k]/m is a field extension of Q; - the ring l[1] contains an ideal PARA, which contains the parametrisation of a point in V(i) lying over w up to the first ord terms; - the string m=l[3] contains the code of the maximal ideal m, by which we have to divide Q[X_1,...,X_r] in order to have the appropriate field extension over which the parametrisation lives; - and if the integer l[2] is N then t has to be replaced by t^1/N in the parametrisation, or alternatively replace t by t^N in the defining ideal - the procedure REQUIRES that the program GFAN is installed on your computer" { def ALTRING=basering; int recursively; // is set 1 if the procedure is called recursively by itself int jj,kk; if (size(#)==2) // this means the precedure has been called recursively { // how many variables are true variables, and how many come // from the field extension // only true variables have to be transformed int anzahlvariablen=#[1]; ideal gesamt_m=std(#[2]); // stores all maxideals used for field extensions // find the zeros of the w-initial ideal and transform the ideal i; // findzeros and basictransformideal need to know how // many of the variables are true variables list m_ring=findzeros(i,ww,anzahlvariablen); list btr=basictransformideal(i,ww,m_ring,anzahlvariablen); } else // the procedure has been called by tropicalLifting { // how many variables are true variables, and how many come from // the field extension only true variables have to be transformed int anzahlvariablen=nvars(basering); ideal gesamt_m; // stores all maxideals used for field extensions // the initial ideal of i has been handed over as #[1] ideal ini=#[1]; // find the zeros of the w-initial ideal and transform the ideal i; // we should hand the t-initial ideal ine to findzeros, // since we know it already list m_ring=findzeros(i,ww,ini); list btr=basictransformideal(i,ww,m_ring); } // check if for the transformation a field extension was necessary, i.e. if the // new variables had to be added to the old base ring; if so, change to the new // ring in which the transformed i and m and the zero a of V(m) live; otherwise // retreive i, a and m from btr if (size(btr)==1) { def PREVIOUSRING=basering; def TRRING=btr[1]; setring TRRING; ideal gesamt_m=imap(PREVIOUSRING,gesamt_m)+m; // add the newly found maximal // ideal to the previous ones } else { i=std(btr[1]); list a=btr[2]; ideal m=btr[3]; gesamt_m=gesamt_m+m; // add the newly found maximal // ideal to the previous ones } // check if there is a solution which has the n-th component zero, // if so, then eliminate the n-th variable from sat(i+x_n,t), // otherwise leave i as it is; // then check if the (remaining) ideal has as solution // where the n-1st component is zero, // and procede as before; do the same for the remaining variables; // this way we make sure that the remaining ideal has // a solution which has no component zero; intvec deletedvariables; // the jth entry is set 1, if we eliminate x_j int numberdeletedvariables; // the number of eliminated variables ideal variablen; // will contain the variables which are not eliminated intvec tw=ww; // in case some variables are deleted, // we have to store the old weight vector deletedvariables[anzahlvariablen]=0; ideal I,LI; i=i+m; // if a field extension was necessary, then i has to be extended by m for (jj=anzahlvariablen-1;jj>=1;jj--) // the variable t is the last one !!! { I=sat(ideal(var(jj)+i),t)[1]; LI=subst(I,var(nvars(basering)),0); //size(deletedvariables)=anzahlvariablen(before elim.) for (kk=1;kk<=size(deletedvariables)-1;kk++) { LI=subst(LI,var(kk),0); } if (size(LI)==0) // if no power of t is in lead(I) { // (where the X(i) are considered as field elements) // get rid of var(jj) i=eliminate(I,var(jj)); deletedvariables[jj]=1; anzahlvariablen--; // if a variable is eliminated, // then the number of true variables drops numberdeletedvariables++; } else { variablen=variablen+var(jj); // non-eliminated true variables are stored } } variablen=invertorder(variablen); // store also the additional variables and t, // since they for sure have not been eliminated for (jj=anzahlvariablen+numberdeletedvariables-1;jj<=nvars(basering);jj++) { variablen=variablen+var(jj); } // if some variables have been eliminated, // then pass to a new ring which has less variables, // but if no variables are left, then we are done def BASERING=basering; if ((numberdeletedvariables>0) and (anzahlvariablen>1)) // if only t remains, { // all true variables are gone execute("ring NEURING=("+charstr(basering)+"),("+string(variablen)+"),(dp("+string(size(variablen)-1)+"),lp(1));"); ideal i=imap(BASERING,i); ideal gesamt_m=imap(BASERING,gesamt_m); } // now we have to compute a point ww on the tropical variety // of the transformed ideal i; // of course, we only have to do so, if we have not yet // reached the order up to which we // were supposed to do our computations if ((ordnung>1) and (anzahlvariablen>1)) // if only t remains, { // all true variables are gone // else we have to use gfan def PREGFANRING=basering; if (nogfan!=1) { // pass to a ring which has variables which are suitable for gfan execute("ring GFANRING=("+charstr(basering)+"),(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P,Q,R,S,T,U,V,W,X,Y,Z),dp;"); ideal phiideal=b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P,Q,R,S,T,U,V,W,X,Y,Z; phiideal[nvars(PREGFANRING)]=a; // map t to a map phi=PREGFANRING,phiideal; ideal i=phi(i); // homogenise the ideal i with the first not yet // used variable in our ring, since gfan // only handles homogenous ideals; in principle // for this one has first to compute a // standard basis of i and homogenise that, // but for the tropical variety (says Anders) // it suffices to homogenise an arbitrary system of generators // i=groebner(i); i=homog(i,maxideal(1)[nvars(PREGFANRING)+1]); // if gfan version >= 0.3.0 is used and the characteristic // is not zero, then write the basering to the output if ((gfanold!=1) and (char(GFANRING)!=0)) { string ringvariablen=varstr(GFANRING); ringvariablen=ringvariablen[1..2*nvars(PREGFANRING)+1]; write(":w /tmp/gfaninput","Z/"+string(char(GFANRING))+"Z["+ringvariablen+"]"); // write the ideal to a file which gfan takes as input and call gfan write(":a /tmp/gfaninput","{"+string(i)+"}"); } else { // write the ideal to a file which gfan takes as input and call gfan write(":w /tmp/gfaninput","{"+string(i)+"}"); } if (gfanold==1) { if (charstr(basering)!="0") { system("sh","gfan_tropicalbasis --mod "+charstr(basering)+" < /tmp/gfaninput > /tmp/gfanbasis"); system("sh","gfan_tropicalintersection < /tmp/gfanbasis > /tmp/gfanoutput"); } else { // system("sh","gfan_tropicalstartingcone < /tmp/gfaninput > /tmp/gfantropstcone"); // system("sh","gfan_tropicaltraverse < /tmp/gfantropstcone > /tmp/gfanoutput"); system("sh","gfan_tropicalbasis < /tmp/gfaninput > /tmp/gfanbasis"); system("sh","gfan_tropicalintersection < /tmp/gfanbasis > /tmp/gfanoutput"); } string trop=read("/tmp/gfanoutput"); setring PREGFANRING; intvec wneu=-1; // this integer vector will store // the point on the tropical variety wneu[nvars(basering)]=0; // for the time being simply stop Singular and set wneu by yourself int goon=1; trop; "CHOOSE A RAY IN THE OUTPUT OF GFAN WHICH HAS ONLY"; "NON-POSITIVE ENTRIES AND STARTS WITH A NEGATIVE ONE,"; "E.G. (-3,-4,-1,-5,0,0,0) - the last entry will always be 0 -"; "THEN TYPE THE FOLLOWING COMMAND IN SINGULAR: wneu=-3,-4,-1,-5,0,0;"; "AND HIT THE RETURN BUTTON TWICE (!!!) - NOTE, THE LAST 0 IS OMITTED"; "IF YOU WANT wneu TO BE TESTED, THEN SET goon=0;"; ~ // THIS IS NOT NECESSARY !!!! IF WE COMPUTE NOT ONLY THE // TROPICAL PREVARIETY // test, if wneu really is in the tropical variety while (goon==0) { if (testw(reduce(i,std(gesamt_m)),wneu,anzahlvariablen)==1) { "CHOOSE A DIFFERENT RAY"; ~ } else { goon=1; } } } else { system("sh","gfan_tropicallifting -n "+string(anzahlvariablen)+" --noMult -c < /tmp/gfaninput > /tmp/gfanoutput"); // read the result from gfan and store it to a string, // which in a later version // should be interpreded by Singular intvec wneu=choosegfanvector(read("/tmp/gfanoutput"),0)[1]; setring PREGFANRING; } } else // if gfan is NOT executed { // if puiseux is set, then we are in the case of a plane curve and can use the command newtonpoly if (puiseux==1) { list NewtP=newtonpoly(i[1]); wneu=NewtP[1]-NewtP[2]; int ggteiler=gcd(wneu[1],wneu[2]); wneu[1]=-wneu[1]/ggteiler; wneu[2]=wneu[2]/ggteiler; if (wneu[1]>0) { wneu=-wneu; } kill NewtP,ggteiler; } else // set wneu by hand { "Set intvec wneu!"; ~ } } } // if we have not yet computed our parametrisation // up to the required order and // zero is not yet a solution, then we have // to go on by calling the procedure recursively; // if all variables were deleted, then i=0 and thus anzahlvariablen==0 if ((ordnung>1) and (anzahlvariablen>1)) { // we call the procedure with the transformed ideal i, // the new weight vector, with the // required order lowered by one, and with // additional parameters, namely the number of // true variables and the maximal ideal that // was computed so far to describe the field extension list PARALIST=tropicalparametriseNoabs(i,wneu,ordnung-1,gfanold,nogfan,puiseux,anzahlvariablen,gesamt_m); // the output will be a ring, in which the // parametrisation lives, and a string, which contains // the maximal ideal that describes the necessary field extension def PARARing=PARALIST[1]; int tweight=-tw[1]*PARALIST[2]; string PARAm=PARALIST[3]; setring PARARing; // if some variables have been eliminated in before, // then we have to insert zeros // into the parametrisation for those variables if (numberdeletedvariables>0) { ideal PARAneu=PARA; int k; for (jj=1;jj<=anzahlvariablen+numberdeletedvariables-1;jj++) // t admits { // no parametrisation if (deletedvariables[jj]!=1) { k++; PARA[jj]=PARAneu[k]; } else { PARA[jj]=poly(0); } } } } // otherwise we are done and we can start to compute // the last step of the parametrisation else { // we store the information on m in a string string PARAm=string(gesamt_m); // we define the weight of t, i.e. in the parametrisation t // has to be replaced by t^1/tweight int tweight=-tw[1]; // if additional variables were necessary, // we introduce them now as parameters; // in any case the parametrisation ring will // have only one variable, namely t, // and its order will be local, so that it // displays the lowest term in t first if (anzahlvariablen+numberdeletedvariables (a_j+x_j)*t^w_j l[2] = list, a_1,...,a_n a point in V(m) l[3] = ideal, the maximal ideal m or l[1] = ring which contains the ideals i and m, and the list a NOTE: the procedure is called from inside the recursive procedure tropicalparametriseNoabs; if it is called in the first recursion, the list # is empty, otherwise #[1] is an integer, the number of true variables x_1,...,x_n; during the procedure we check if a field extension is necessary to express a zero (a_1,...,a_n) of m; if so, we have to introduce new variables and a list containing a ring is returned, otherwise the list containing i, a and m is returned; the ideal m will be changed during the procedure since all variables which reduce to a polynomial in X_1,...,X_k modulo m will be eliminated, while the others are replaced by new variables X_k+1,...,X_k'" { def BASERING=basering; // the variable t is acutally the LAST variable !!! int j,k; // counters // store the weighted degrees of the elements of i in an integer vector intvec wdegs; for (j=1;j<=size(i);j++) { wdegs[j]=deg(i[j],intvec(w[2..size(w)],w[1])); } // how many variables are true variables, // and how many come from the field extension // only real variables have to be transformed if (size(#)>0) { int anzahlvariablen=#[1]; } else { int anzahlvariablen=nvars(basering); } // int zaehler; // testvariable // get the information if any new variables are needed from m_ring int neuvar=m_ring[2]; // number of variables which have to be added intvec neuevariablen=m_ring[3]; // [i]=1, if for the ith comp. // of a zero of m a new variable is needed def MRING=m_ring[1]; // MRING contains a and m list a=imap(MRING,a); // (a_1,...,a_n)=(a[2],...,a[n+1]) will be // a common zero of the ideal m ideal m=std(imap(MRING,m)); // m is zero, if neuvar==0, // otherwise we change the ring anyway // if a field extension is needed, then extend the polynomial // ring by new variables X_k+1,...,X_k'; if (neuvar>0) { // change to a ring where for each variable needed // in m a new variable has been introduced ideal variablen; // extract the variables x_1,...,x_n for (j=1;jt^w_0 { a[j]=X(zaehler); zaehler++; } for (k=1;k<=size(i);k++) { if (j!=1) // corresponds to x_(j-1) -- note t is the last variable { i[k]=substitute(i[k],var(j-1),(a[j]+var(j-1))*t^(-w[j])); } else // corresponds to t { i[k]=substitute(i[k],t,t^(-w[j])); } } } // we can divide the jth generator of i by t^-wdegs[j] for (j=1;j<=ncols(i);j++) { if (wdegs[j]<0) // if wdegs[j]==0 there is no need to divide, and { // we made sure that it is no positive i[j]=i[j]/t^(-wdegs[j]); } } // since we want to consider i now in the ring // (Q[X_1,...,X_k']/m)[t,x_1,...,x_n] // we can reduce i modulo m, so that "constant terms" // which are "zero" since they // are in m will disappear; simplify(...,2) then really removes them i=simplify(reduce(i,m),2); // if new variables have been introduced, we have // to return the ring containing i, a and m // otherwise we can simply return a list containing i, a and m if (neuvar>0) { export(m); export(i); export(a); list returnlist=TRANSFORMRING; return(returnlist); } else { return(list(i,a,m)); } } ///////////////////////////////////////////////////////////////////////// static proc testw (ideal i,intvec w,int anzahlvariablen,list #) "USAGE: testw(i,w,n); i ideal, w intvec, n number ASSUME: i is an ideal in Q[t,x_1,...,x_n,X_1,...,X_k] and w=(w_0,...,w_n,0,...,0) RETURN: int b, 0 if the t-initial ideal of i considered in Q(X_1,...,X_k)[t,x_1,...,x_n] is monomial free, 1 else NOTE: the procedure is called by tinitialform" { // compute the t-initial of i // - the last 1 just means that the variable t is the last variable in the ring if (size(#)==0) { ideal tin=tInitialIdeal(i,w); } else { ideal tin=tInitialIdeal(i,w,1); } int j; ideal variablen; poly monom=1; if (size(#)==0) { for (j=2;j<=anzahlvariablen;j++) { variablen=variablen+var(j); monom=monom*var(j); } ideal Parameter; for (j=anzahlvariablen+1;j<=nvars(basering);j++) { Parameter=Parameter+var(j); } } else { for (j=1;j<=anzahlvariablen-1;j++) { variablen=variablen+var(j); monom=monom*var(j); } ideal Parameter; for (j=anzahlvariablen;j<=nvars(basering)-1;j++) { Parameter=Parameter+var(j); } } def BASERING=basering; if (anzahlvariablensize(w)) or (size(w)==1)) { return(w); } if (i==1) { return(w[2..size(w)]); } if (i==size(w)) { return(w[1..size(w)-1]); } else { intvec erg=w[1..i-1],w[i+1..size(w)]; return(erg); } } ///////////////////////////////////////////////////////////////////////// static proc invertorder (ideal i) "USAGE: intvertorder(i); i ideal RETURN: ideal, the order of the entries of i has been inverted NOTE: the procedure is called by tropicalparametriseNoabs" { ideal ausgabe; for (int j=size(i);j>=1;j--) { ausgabe[size(i)-j+1]=i[j]; } return(ausgabe); } ///////////////////////////////////////////////////////////////////////// static proc ordermaximalidealsNoabs (list minassi,int anzahlvariablen) "USAGE: ordermaximalidealsNoabs(minassi); minassi list ASSUME: minassi is a list of maximal ideals RETURN: list, the procedure orders the maximal ideals in minassi according to how many new variables are needed to describe the zeros of the ideal l[j][1] = jth maximal ideal l[j][2] = the number of variables needed l[j][3] = intvec, if for the kth variable a new variable is needed to define the corresponding zero of l[j][1], then the k+1st entry is one l[j][4] = list, if for the kth variable no new variable is needed to define the corresponding zero of l[j][1], then its value is the k+1st entry NOTE: if a maximal ideal contains a variable, it is removed from the list; the procedure is called by findzeros" { int j,k,l; int neuvar; // number of new variables needed (for each maximal ideal) int pruefer; // is set one if a maximal ideal contains a variable list minassisort; // will contain the output for (j=1;j<=size(minassi);j++){minassisort[j]=0;} // initialise minassisort // to fix its initial length list zwischen; // needed for reordering list a;// (a_1,...,a_n)=(a[2],...,a[n+1]) will be a common zero of the ideal m poly nf; // normalform of a variable w.r.t. m intvec neuevariablen; // ith entry is 1, if for the ith component of a zero // of m a new variable is needed intvec testvariablen; // integer vector of length n=number of variables // compute for each maximal ideal the number of new variables, // which are needed to describe // its zeros -- note, a new variable is needed if modulo // the maximal ideal it does not reduce // to something which only depends on the following variables; // if no new variable is needed, then store the value // a variable reduces to in the list a; for (j=size(minassi);j>=1;j--) { a[1]=poly(0); // the first entry in a and in neuevariablen // corresponds to the variable t, neuevariablen[1]=0; // which is not present in the INITIALRING neuvar=0; minassi[j]=std(minassi[j]); for (k=1;(k<=anzahlvariablen-1) and (pruefer==0);k++) { nf=reduce(var(k),minassi[j]); // if a variable reduces to zero, then the maximal ideal // contains a variable and we can delete it if (nf==0) { pruefer=1; } // set the entries of testvariablen corresponding to the first // k variables to 1, and the others to 0 for (l=1;l<=nvars(basering);l++) { if (l<=k) { testvariablen[l]=1; } else { testvariablen[l]=0; } } // if the variable x_j reduces to a polynomial // in x_j+1,...,x_n,X_1,...,X_k modulo m // then we can eliminate x_j from the maximal ideal // (since we do not need any // further field extension to express a_j) and a_j // will just be this normalform, // otherwise we have to introduce a new variable in order to express a_j; if (scalarproduct(leadexp(nf),testvariablen)==0) { a[k+1]=nf; // a_k is the normal form of the kth variable modulo m neuevariablen[k+1]=0; // no new variable is needed } else { a[k+1]=poly(0); // a_k is set to zero until we have // introduced the new variable neuevariablen[k+1]=1; neuvar++; // the number of newly needed variables is raised by one } } // if the maximal ideal contains a variable, we simply delete it if (pruefer==0) { minassisort[j]=list(minassi[j],neuvar,neuevariablen,a); } // otherwise we store the information on a, // neuevariablen and neuvar together with the ideal else { minassi=delete(minassi,j); minassisort=delete(minassisort,j); pruefer=0; } } // sort the maximal ideals ascendingly according to // the number of new variables needed to // express the zero of the maximal ideal for (j=2;j<=size(minassi);j++) { l=j; for (k=j-1;k>=1;k--) { if (minassisort[l][2]1) and (N==1)) { dp=dp+displaycoef(koeffizienten[j-ord(p)+1])+"*t^"+string(j-(N*wj)/w1); } if (((j-(N*wj)/w1)>1) and (N!=1)) { dp=dp+displaycoef(koeffizienten[j-ord(p)+1])+"*t^("+string(j-(N*wj)/w1)+"/"+string(N)+")"; } if (((j-(N*wj)/w1)<-1) and (N==1)) { dp=dp+displaycoef(koeffizienten[j-ord(p)+1])+"*1/t^"+string(wj-j); } if (((j-(N*wj)/w1)<-1) and (N!=1)) { dp=dp+displaycoef(koeffizienten[j-ord(p)+1])+"*1/t^("+string((N*wj)/w1-j)+"/"+string(N)+")"; } if (j0) // noAbs was used { execute("ring NOQRing=("+string(char(LIFTRing))+"),("+varstr(basering)+","+parstr(LIFTRing)+"),dp;"); execute("qring TESTRing=std("+#[1]+");"); ideal i=imap(BASERING,i); } else // absolute primary decomposition was used { setring LIFTRing; string @mp= string(minpoly); execute("ring TESTRing=("+charstr(LIFTRing)+"),("+varstr(BASERING)+"),dp;"); execute("minpoly= " + @mp + ";"); ideal i=imap(BASERING,i); } } else { def TESTRing=BASERING; setring TESTRing; } // map the ideal LIFT to the TESTRing and substitute the solutions in i ideal liftideal=imap(LIFTRing,LIFT); for (int j=1;j<=ncols(liftideal);j++) { i=substitute(i,var(j),liftideal[j]); } // map the resulting i back to LIFTRing; // if noAbs, then reduce i modulo the maximal ideal // before going back to LIFTRing if ((size(parstr(LIFTRing))!=0) and size(#)>0) { i=reduce(i,std(0)); setring LIFTRing; ideal i=imap(TESTRing,i); } else { setring LIFTRing; ideal i=imap(TESTRing,i); } list SUBSTTEST; // replace t by t^1/N for (j=1;j<=ncols(i);j++) { SUBSTTEST[j]=displaypoly(i[j],N,0,1); } return(SUBSTTEST); } /////////////////////////////////////////////////////////////////////////////// /// Procedures used in tropicalLifting when using absolute primary decomposition /// - tropicalparametrise /// - eliminatecomponents /// - findzerosAndBasictransform /// - ordermaximalideals /////////////////////////////////////////////////////////////////////////////// static proc tropicalparametrise (ideal i,intvec ww,int ordnung,int gfanold,int findall,int nogfan,int puiseux,list #) "USAGE: tropicalparametrise(i,tw,ord,gf,fa,ng,pu[,#]); i ideal, tw intvec, ord int, gf,fa,ng,pu int, # opt. list ASSUME: - i is an ideal in Q[t,x_1,...,x_n,@a], tw=(w_0,w_1,...,w_n,0) and (w_0,...,w_n,0) is in the tropical variety of i, and ord is the order up to which a point in V(i) over K{{t}} lying over w shall be computed; - moreover, @a should be not be there if the procedure is not called recursively; - the point in the tropical variety is supposed to lie in the NEGATIVE orthant; - the ideal is zero-dimensional when considered in (K(t)[@a]/m)[x_1,...,x_n], where m=#[2] is a maximal ideal in K[@a]; - gf is 0 if version 2.2 or larger is used and it is 1 else - fa is 1 if all solutions should be found - ng is 1 if gfan should not be executed - pu is 1 if n=1 and i is generated by one polynomial and newtonpoly is used to find a point in the tropical variety instead of gfan RETURN: list, l[1] = ring K[@a]/m[[t]] l[2] = int NOTE: - the procedure is also called recursively by itself, and if it is called in the first recursion, the list # is empty, otherwise #[1] is an integer, one more than the number of true variables x_1,...,x_n, and #[2] will contain the maximal ideal m in the variable @a by which the ring K[t,x_1,...,x_n,@a] should be divided to work correctly in L[t,x_1,...,x_n] where L=Q[@a]/m is a field extension of K - the ring l[1] contains an ideal PARA, which contains the parametrisation of a point in V(i) lying over w up to the first ord terms; - the string m contains the code of the maximal ideal m, by which we have to divide K[@a] in order to have the appropriate field extension over which the parametrisation lives; - and if the integer l[3] is N then t has to be replaced by t^1/N in the parametrisation, or alternatively replace t by t^N in the defining ideal - the procedure REQUIRES that the program GFAN is installed on your computer" { def BASERING=basering; int recursively; // is set 1 if the procedure is called recursively by itself int ii,jj,kk,ll,jjj,kkk,lll; if (typeof(#[1])=="int") // this means the precedure has been { // called recursively // how many variables are true variables, and how many come // from the field extension // only true variables have to be transformed int anzahlvariablen=#[1]; // find the zeros of the w-initial ideal and transform the ideal i; // findzeros and basictransformideal need to know // how many of the variables are true variables // and do the basic transformation as well list zerolist=#[2]; list trring=findzerosAndBasictransform(i,ww,zerolist,findall,anzahlvariablen); } else // the procedure has been called by tropicalLifting { // how many variables are true variables, and // how many come from the field extension // only true variables have to be transformed int anzahlvariablen=nvars(basering); list zerolist; // will carry the zeros which are //computed in the recursion steps // the initial ideal of i has been handed over as #[1] ideal ini=#[1]; // find the zeros of the w-initial ideal and transform the ideal i; // we should hand the t-initial ideal ine to findzeros, // since we know it already; // and do the basic transformation as well list trring=findzerosAndBasictransform(i,ww,zerolist,findall,ini); } list wneulist; // carries all newly computed weight vector intvec wneu; // carries the newly computed weight vector int tweight; // carries the weight of t list PARALIST; // carries the result when tropicalparametrise // is recursively called list eliminaterings; // carries the result of eliminatecomponents intvec deletedvariables; // contains inform. which variables have been deleted deletedvariables[anzahlvariablen-1]=0; // initialise this vector int numberdeletedvariables; // the number of deleted variables int oldanzahlvariablen=anzahlvariablen; // anzahlvariablen for later reference list liftings,partliftings; // the computed liftings (all resp. partly) // consider each ring which has been returned when computing the zeros of the // the t-initial ideal, equivalently, consider // each zero which has been returned for (jj=1;jj<=size(trring);jj++) { def TRRING=trring[jj]; setring TRRING; // check if a certain number of components lead to suitable // solutions with zero components; // compute them all if findall==1 eliminaterings=eliminatecomponents(i,m,oldanzahlvariablen,findall,oldanzahlvariablen-1,deletedvariables); // consider each of the rings returned by eliminaterings ... // there is at least one !!! for (ii=1;ii<=size(eliminaterings);ii++) { // #variables which have been eliminated numberdeletedvariables=oldanzahlvariablen-eliminaterings[ii][2]; // #true variables which remain (including t) anzahlvariablen=eliminaterings[ii][2]; // a 1 in this vector says that variable was eliminated deletedvariables=eliminaterings[ii][3]; setring TRRING; // set TRRING - this is important for the loop // pass the ring computed by eliminatecomponents def PREGFANRING=eliminaterings[ii][1]; setring PREGFANRING; poly m=imap(TRRING,m); // map the maximal ideal to this ring list zero=imap(TRRING,zero); // map the vector of zeros to this ring // now we have to compute a point ww on the tropical // variety of the transformed ideal i; // of course, we only have to do so, if we have // not yet reached the order up to which we // were supposed to do our computations if ((ordnung>1) and (anzahlvariablen>1)) // if only t remains, { // all true variables are gone if (nogfan!=1) { // pass to a ring which has variables which are suitable for gfan execute("ring GFANRING=("+string(char(basering))+"),(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P,Q,R,S,T,U,V,W,X,Y,Z),dp;"); ideal phiideal=b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P,Q,R,S,T,U,V,W,X,Y,Z; phiideal[nvars(PREGFANRING)]=a; // map t to a map phi=PREGFANRING,phiideal; ideal II=phi(i); // homogenise the ideal II with the first not yet // used variable in our ring, since gfan // only handles homogenous ideals; in principle for this // one has first to compute a // standard basis of II and homogenise that, // but for the tropical variety (says Anders) // it suffices to homogenise an arbitrary system of generators // II=groebner(II); II=homog(II,maxideal(1)[nvars(PREGFANRING)+1]); // if gfan version >= 0.3.0 is used and the characteristic // is not zero, then write the basering to the output if ((gfanold!=1) and (char(GFANRING)!=0)) { string ringvariablen=varstr(GFANRING); ringvariablen=ringvariablen[1..2*nvars(PREGFANRING)+1]; write(":w /tmp/gfaninput","Z/"+string(char(GFANRING))+"Z["+ringvariablen+"]"); // write the ideal to a file which gfan takes as input and call gfan write(":a /tmp/gfaninput","{"+string(II)+"}"); } else { // write the ideal to a file which gfan takes as input and call gfan write(":w /tmp/gfaninput","{"+string(II)+"}"); } if (gfanold==1) { if (charstr(basering)!="0") { system("sh","gfan_tropicalbasis --mod "+charstr(basering)+" < /tmp/gfaninput > /tmp/gfanbasis"); system("sh","gfan_tropicalintersection < /tmp/gfanbasis > /tmp/gfanoutput"); } else { // system("sh","gfan_tropicalstartingcone < /tmp/gfaninput > /tmp/gfantropstcone"); // system("sh","gfan_tropicaltraverse < /tmp/gfantropstcone > /tmp/gfanoutput"); system("sh","gfan_tropicalbasis < /tmp/gfaninput > /tmp/gfanbasis"); system("sh","gfan_tropicalintersection < /tmp/gfanbasis > /tmp/gfanoutput"); } string trop=read("/tmp/gfanoutput"); setring PREGFANRING; wneu=-1; // this integer vector will store the point on the tropical variety wneu[nvars(basering)]=0; // for the time being simply stop Singular and set wneu by yourself int goon=1; trop; "CHOOSE A RAY IN THE OUTPUT OF GFAN WHICH HAS ONLY"; "NON-POSITIVE ENTRIES AND STARTS WITH A NEGATIVE ONE,"; "E.G. (-3,-4,-1,-5,0,0,0) - the last entry will always be 0 -"; "THEN TYPE THE FOLLOWING COMMAND IN SINGULAR: wneu=-3,-4,-1,-5,0,0;"; "AND HIT THE RETURN BUTTON - NOTE, THE LAST 0 IS OMITTED"; "IF YOU WANT wneu TO BE TESTED, THEN SET goon=0;"; ~ // THIS IS NOT NECESSARY !!!! IF WE COMPUTE NOT ONLY // THE TROPICAL PREVARIETY // test, if wneu really is in the tropical variety while (goon==0) { if (testw(reduce(i,std(gesamt_m)),wneu,anzahlvariablen)==1) { "CHOOSE A DIFFERENT RAY"; ~ } else { goon=1; } } wneulist[1]=wneu; } else { system("sh","gfan_tropicallifting -n "+string(anzahlvariablen)+" --noMult -c < /tmp/gfaninput > /tmp/gfanoutput"); // read the result from gfan and store it to a string, // which in a later version // should be interpreded by Singular wneulist=choosegfanvector(read("/tmp/gfanoutput"),findall); setring PREGFANRING; } kill GFANRING; } else // if gfan is NOT executed { // if puiseux is set, then we are in the case of a plane curve and can use the command newtonpoly if (puiseux==1) { if (nvars(basering)>2) // if the basering has a parameter { // we have to pass to a ring with two variables to compute the newtonpoly def PRENPRing=basering; poly NPpoly=i[1]; ring NPRING=(0,@a),(x(1),t),ds; poly NPpoly=imap(PRENPRing,NPpoly); list NewtP=newtonpoly(NPpoly); setring PRENPRing; kill NPRING; } else { list NewtP=newtonpoly(i[1]); } for (jjj=1;jjj<=size(NewtP)-1;jjj++) { wneu=NewtP[jjj]-NewtP[jjj+1]; int ggteiler=gcd(wneu[1],wneu[2]); wneu[1]=-wneu[1]/ggteiler; wneu[2]=wneu[2]/ggteiler; if (wneu[1]>0) { wneu=-wneu; } if (nvars(basering)>2) { wneu[3]=0; } wneulist[jjj]=wneu; } kill NewtP,ggteiler; } else // we have to set the points in the tropical variety by hand { if (findall==1) { "Set wneulist!"; ~ } else { "Set intvec wneu!"; ~ wneulist[1]=wneu; } } } } // if we have not yet computed our parametrisation up to // the required order and // zero is not yet a solution, then we have to go on // by calling the procedure recursively; // if all variables were deleted, then i=0 and thus anzahlvariablen==0 lll=0; if ((ordnung>1) and (anzahlvariablen>1)) { partliftings=list(); // initialise partliftings // we call the procedure with the transformed // ideal i, the new weight vector, with the // required order lowered by one, and with // additional parameters, namely the number of // true variables and the maximal ideal that // was computed so far to describe the field extension for (kk=1;kk<=size(wneulist);kk++) { wneu=wneulist[kk]; PARALIST=tropicalparametrise(i,wneu,ordnung-1,gfanold,findall,nogfan,puiseux,anzahlvariablen,zero); // the output will be a ring, in which the // parametrisation lives, and a string, which contains // the maximal ideal that describes the necessary field extension for (ll=1;ll<=size(PARALIST);ll++) { def PARARing=PARALIST[ll][1]; tweight=-ww[1]*PARALIST[ll][2]; setring PARARing; // if some variables have been eliminated // in before, then we have to insert zeros // into the parametrisation for those variables if (numberdeletedvariables>0) { ideal PARAneu=PARA; kkk=0; for (jjj=1;jjj<=anzahlvariablen+numberdeletedvariables-1;jjj++) { // t admits no parametrisation if (deletedvariables[jjj]!=1) { kkk++; PARA[jjj]=PARAneu[kkk]; } else { PARA[jjj]=poly(0); } } } lll++; partliftings[lll]=list(PARARing,tweight,wneu); setring PREGFANRING; kill PARARing; } } } // otherwise we are done and we can start // to compute the last step of the parametrisation else { // we define the weight of t, i.e. in the // parametrisation t has to be replaced by t^1/tweight tweight=-ww[1]; // if additional variables were necessary, // we introduce them now as parameters; // in any case the parametrisation ring will // have only one variable, namely t, // and its order will be local, so that it // displays the lowest term in t first if (anzahlvariablen so that // there is a solution which has strictly positive valuation // ADDENDUM: // if i (without m) has only one polynomial, then we can divide // i by t as long as possible to compute the saturation with respect to t /* // DER NACHFOLGENDE TEIL IST MUELL UND WIRD NICHT MEHR GAMACHT // for the test we simply compute the leading ideal // and set all true variables zero; // since the ordering was an elimination ordering // with respect to (@a if present and) t // there remains something not equal to zero // if and only if there is polynomial which only // depends on t (and @a if present), but that is // a unit in K{{t}}, which would show that there // is no solution with the component lastvar zero; // however, we also have to ensure that if there // exists a solution with the component lastvar // equal to zero then this component has a // valuation with all strictly positive values!!!!; // for this we can either saturate the ideal // after elimination with respect to // and see if the saturated ideal is contained in +m, // or ALTERNATIVELY we can pass to the // ring 0,(t,x_1,...,x_n,@a),(ds(n+1),dp(1)), // compute a standard basis of the elimination // ideal (plus m) there and check if the dimension 1 // (since we have not omitted the variable lastvar, // this means that we have the ideal // generated by t,x_1,...,[no x_lastvar],...,x_n // and this defines NO curve after omitting x_lastvar) I=std(ideal(var(lastvar)+i)); // first test, LI=lead(reduce(I,std(m))); //size(deletedvariables)=anzahlvariablen(before elimination) for (j=1;j<=anzahlvariablen-1;j++) { LI=subst(LI,var(j),0); } if (size(LI)==0) // if no power of t is in lead(I) (where @a is considered as a field element) */ if (size(i)-size(m)!=1) { I=reduce(sat(std(ideal(var(lastvar)+i)),t)[1],std(m)); // get rid of the minimal // polynomial for the test } else { I=subst(i,var(lastvar),0); while (subst(I[1],t,0)==0) { I[1]=I[1]/t; } I=reduce(I,std(m)); } LI=subst(I,var(nvars(basering)),0); //size(deletedvariables)=anzahlvariablen(before elimination) for (j=1;j<=anzahlvariablen-1;j++) { LI=subst(LI,var(j),0); } if (size(LI)==0) // the saturation lives in the maximal { // ideal generated by t and the x_i // get rid of var(lastvar) I=eliminate(I,var(lastvar))+m; // add the minimal polynomial again // store the information which variable has been eliminated intvec newdeletedvariables=deletedvariables; newdeletedvariables[lastvar]=1; // pass to a new ring whith one variable less if (anzahlvariablen>2) { string elring="ring ELIMINATERING=("+charstr(basering)+"),("+string(simplify(reduce(maxideal(1),std(var(lastvar))),2))+"),(dp("+string(anzahlvariablen-2)+"),"; } else { string elring="ring ELIMINATERING=("+charstr(basering)+"),("+string(simplify(reduce(maxideal(1),std(var(lastvar))),2))+"),("; } if (anzahlvariablen1) { eliminatelist=eliminatecomponents(i,imap(BASERING,m),anzahlvariablen-1,findall,lastvar-1,newdeletedvariables); } else { export(i); eliminatelist=list(list(ELIMINATERING,anzahlvariablen-1,newdeletedvariables)); } setring BASERING; } // next we have to test if there is also a solution // which has the variable lastvar non-zero; // to do so, we saturate with respect to this variable; // since when considered over K{{t}} // the ideal is zero dimensional, this really removes // all solutions with that component zero; // the checking is then done as above // ADDENDUM // if the ideal i (without m) is generated by a single polynomial // then we saturate by successively dividing by the variable if (size(i)-size(m)==1) { while (subst(i[1],var(lastvar),0)==0) { i[1]=i[1]/var(lastvar); } } else { i=std(sat(i,var(lastvar))[1]); } I=reduce(i,std(m)); // "remove" m from i // test first, if i still is in the ideal -- if not, then // we know that i has no longer a point in the tropical // variety with all components negative LI=subst(I,var(nvars(basering)),0); for (j=1;j<=anzahlvariablen-1;j++) // set all variables { // t,x_1,...,x_n equal to zero LI=subst(LI,var(j),0); } if (size(LI)==0) // the entries of i have no constant part { // check now, if the leading ideal of i contains an element // which only depends on t LI=lead(I); //size(deletedvariables)=anzahlvariablen(before elimination) for (j=1;j<=anzahlvariablen-1;j++) { LI=subst(LI,var(j),0); } if (size(LI)==0) // if no power of t is in lead(i) { // (where @a is considered as a field element) // if not yet all variables have been tested, go on with the // next smaller variable // else prepare the neccessary information for return if (lastvar>1) { saturatelist=eliminatecomponents(i,m,anzahlvariablen,findall,lastvar-1,deletedvariables); } else { execute("ring SATURATERING=("+charstr(basering)+"),("+varstr(basering)+"),("+ordstr(basering)+");"); ideal i=imap(BASERING,i); export(i); setring BASERING; saturatelist=list(list(SATURATERING,anzahlvariablen,deletedvariables)); } } } return(eliminatelist+saturatelist); } else // only one solution is searched for, we can do a simplified { // version of the above // check if there is a solution which has the n-th component // zero and with positive valuation, // if so, then eliminate the n-th variable from // sat(i+x_n,t), otherwise leave i as it is; // then check if the (remaining) ideal has as // solution where the n-1st component is zero ..., // and procede as before; do the same for the remaining variables; // this way we make sure that the remaining ideal has // a solution which has no component zero; ideal variablen; // will contain the variables which are not eliminated for (j=anzahlvariablen-1;j>=1;j--) // the variable t is the last one !!! { I=sat(ideal(std(var(j)+i)),t)[1]; LI=subst(I,var(nvars(basering)),0); //size(deletedvariables)=anzahlvariablen-1(before elimination) for (k=1;k<=size(deletedvariables);k++) { LI=subst(LI,var(k),0); } if (size(LI)==0) // if no power of t is in lead(I) { // (where the X(i) are considered as field elements) // get rid of var(j) i=eliminate(I,var(j)); deletedvariables[j]=1; anzahlvariablen--; // if a variable is eliminated, // then the number of true variables drops } else { variablen=variablen+var(j); // non-eliminated true variables are stored } } variablen=invertorder(variablen); // store also the additional variable and t, // since they for sure have not been eliminated for (j=size(deletedvariables)+1;j<=nvars(basering);j++) { variablen=variablen+var(j); } // if there are variables left, then pass to a ring which // only realises these variables else we are done if (anzahlvariablen>1) { string elring="ring ELIMINATERING=("+charstr(basering)+"),("+string(variablen)+"),(dp("+string(anzahlvariablen-1)+"),"; } else { string elring="ring ELIMINATERING=("+charstr(basering)+"),("+string(variablen)+"),("; } if (size(deletedvariables)+1 t^-w_l*(zero_l+x_l) for (l=1;l<=anzahlvariablen;l++) { for (k=1;k<=size(i);k++) { if (l!=1) // corresponds to x_(l-1) -- note t is the last variable { i[k]=subst(i[k],var(l-1),(zeroneu[l]+var(l-1))*t^(-w[l])); } else // corresponds to t { i[k]=subst(i[k],var(nvars(basering)),var(nvars(basering))^(-w[l])); } } } // we can divide the lth generator of i by t^-wdegs[l] for (l=1;l<=ncols(i);l++) { if (wdegs[l]<0) // if wdegs[l]==0 there is no need to divide, { // and we made sure that it is no positive i[l]=i[l]/t^(-wdegs[l]); } } // since we want to consider i now in the ring (Q[@a]/m)[t,x_1,...,x_n] // we can reduce i modulo m, so that "constant terms" // which are "zero" since they // are in m will disappear; simplify(...,2) then really removes them; // finally we add the minimal polynomial i=simplify(ideal(reduce(i,std(m))+m),2); export(i); export(m); export(zero); extensionringlist[j]=EXTENSIONRING; kill EXTENSIONRING; setring ABSRING; } return(extensionringlist); } ///////////////////////////////////////////////////////////////////////// static proc ordermaximalideals (list minassi,int anzahlvariablen) "USAGE: ordermaximalideals(minassi); minassi list ASSUME: minassi is a list of maximal ideals (together with the information how many conjugates it has), where the first polynomial is the minimal polynomial of the last variable in the ring which is considered as a parameter RETURN: list, the procedure orders the maximal ideals in minassi according to how many new variables are needed (namely none or one) to describe the zeros of the ideal, and accordingly to the degree of the corresponding minimal polynomial l[j][1] = the minimal polynomial for the jth maximal ideal l[j][2] = list, the k+1st entry is the kth coordinate of the zero described by the maximal ideal in terms of the last variable l[j][3] = poly, the replacement for the old variable @a NOTE: if a maximal ideal contains a variable, it is removed from the list; the procedure is called by findzerosAndBasictransform" { int j,k,l; int pruefer; // is set one if a maximal ideal contains a variable list minassisort; // will contain the output for (j=1;j<=size(minassi);j++){minassisort[j]=0;} // initialise minassisort // to fix its initial length list zwischen; // needed for reordering list zero; // (a_1,...,a_n)=(zero[2],...,zero[n+1]) will be // a common zero of the ideal m poly nf; // normalform of a variable w.r.t. m poly minimalpolynomial; // minimal polynomial for the field extension poly parrep; // the old variable @a possibly has to be replaced by a new one // compute for each maximal ideal the number of new variables, which are // needed to describe its zeros -- note, a new variable is needed // if the first entry of minassi[j][1] is not the last variable // store the value a variable reduces to in the list a; for (j=size(minassi);j>=1;j--) { minimalpolynomial=minassi[j][1][1]; if (minimalpolynomial==var(nvars(basering))) { minimalpolynomial=0; } zero[1]=poly(0); // the first entry in zero and in // neuevariablen corresponds to the variable t, minassi[j][1]=std(minassi[j][1]); for (k=1;(k<=anzahlvariablen-1) and (pruefer==0);k++) { // zero_k+1 is the normal form of the kth variable modulo m zero[k+1]=reduce(var(k),minassi[j][1]); // if a variable reduces to zero, then the maximal // ideal contains a variable and we can delete it if (zero[k+1]==0) { pruefer=1; } } // if anzahlvariablen=1;k--) { if (deg(minassisort[l][1])=2;i--) { for (j=i-1;j>=1;j--) { if ((eckpunkte[i][1]==eckpunkte[j][1]) and (eckpunkte[i][2]==eckpunkte[j][2])) { eckpunkte=delete(eckpunkte,i); j=0; } } } return(eckpunkte); } ///////////////////////////////////////////////////////////////////////// static proc bunchOfLines (def tp,list #) "USAGE: bunchOfLines(tp[,#]); tp list, # optional list ASSUME: tp is represents an ideal in Q[x,y] representing a tropical polynomial (in the form of the output of the procedure tropicalise) defining a bunch of ordinary lines in the plane, i.e. the Newton polygone is a line segment RETURN: list, see the procedure tropicalCurve for an explanation of the syntax of this list NOTE: - the tropical curve defined by tp will consist of a bunch of parallel lines and throughout the procedure a list with the name bunchoflines is computed, which represents these lines and has the following interpretation: list, each entry corresponds to a vertex in the tropical plane curve defined by tp l[i][1] = the equation of the ith line in the tropical curve l[i][2] = a polynomial whose monimials mark the vertices in the Newton polygon corresponding to the entries in tp which take the common minimum at the ith vertex - the information in l[i][2] represents the subdivision of the Newton polygon of tp (respectively a polynomial which defines tp); - if # is non-empty and #[1] is the string 'max', then in the computations minimum and maximum are exchanged; - if # is non-empty and the #[1] is not a string, only the vertices will be computed and the information on the Newton subdivision will be omitted; - here the tropical polynomial is supposed to be the MINIMUM of the linear forms in tp, unless the optional input #[1] is the string 'max' - the procedure is called from tropicalCurve" { // introduce necessary variables list oldtp=tp; // save the old entries of tp int i,j,k,l; ideal punkt; poly px; poly wert,vergleich; poly newton; list bunchoflines; int e=1; option(redSB); // find the direction of the line segment in the Newton polygon intvec direction=leadexp(detropicalise(tp[1]))-leadexp(detropicalise(tp[2])); direction=direction/gcd(direction[1],direction[2]); // change the coordinates in such a way, that the Newton polygon // lies on the x-axis if (direction[1]==0) // there is no x-term - exchange x and y { // and get rid of the new y part for (i=1;i<=size(tp);i++) { tp[i]=substitute(tp[i],var(1),0,var(2),var(1)); } } else { for (i=1;i<=size(tp);i++) { tp[i]=substitute(tp[i],var(2),0); } } // For all tuples (i,j) of entries in tp check if they attain // at their point of intersection // the minimal possible value for all linear forms in tp for (i=1;i<=size(tp)-1;i++) { for(j=i+1;j<=size(tp);j++) { punkt=std(ideal(tp[i]-tp[j])); px=-leadcoef(punkt[1]-lead(punkt[1]))/leadcoef(punkt[1]); l=1; wert=substitute(tp[i],var(1),px); newton=0; vergleich=wert; while((l<=size(tp)) and (vergleiche(wert,vergleich,#))) { vergleich=substitute(tp[l],var(1),px); if (vergleich==wert) { newton=newton+detropicalise(oldtp[l]); } l++; } if ((l==size(tp)+1) and (vergleiche(wert,vergleich,#))) { bunchoflines[e]=list((oldtp[i]-oldtp[j])/leadcoef(oldtp[i]-oldtp[j]),newton); e++; } } } // if a vertex appears several times, only its first occurence will be kept for (i=size(bunchoflines);i>=2;i--) { for (j=i-1;j>=1;j--) { if (bunchoflines[i][1]==bunchoflines[j][1]) { bunchoflines=delete(bunchoflines,i); j=0; } } } // sort the lines in an descending way according to the leading // exponent of the polynomial // defining the Newton polygone list nbol; list maximum; while (size(bunchoflines)!=0) { j=1; maximum=bunchoflines[1]; for (i=2;i<=size(bunchoflines);i++) { if (deg(bunchoflines[i][2])=1;i--) { if (vv[1,i]==vv[1,i+1]) { vv[2,i]=vv[2,i]+vv[2,i+1]; vv=intmatcoldelete(vv,i+1); } if (vv[1,i]==0) { vv=intmatcoldelete(vv,i); } } return(vv); } ///////////////////////////////////////////////////////////////////////// static proc sortintvec (intvec w) "USAGE: sortintvec(v); v intvec RETURN: intvec, the entries of v are ordered in an ascending way NOTE: not called at all" { int j,k,stop; intvec v=w[1]; for (j=2;j<=size(w);j++) { k=1; stop=0; while ((k<=size(v)) and (stop==0)) { if (v[k]ncols(w)) or (ncols(w)==1)) { return(w); } if (i==1) { intmat M[nrows(w)][ncols(w)-1]=w[1..nrows(w),2..ncols(w)]; return(M); } if (i==ncols(w)) { intmat M[nrows(w)][ncols(w)-1]=w[1..nrows(w),1..ncols(w)-1]; return(M); } else { intmat M[nrows(w)][i-1]=w[1..nrows(w),1..i-1]; intmat N[nrows(w)][ncols(w)-i]=w[1..nrows(w),i+1..ncols(w)]; return(intmatconcat(M,N)); } } ///////////////////////////////////////////////////////////////////////// static proc intmatconcat (intmat M,intmat N) "USAGE: intmatconcat(M,N); M,N intmat RETURN: intmat, M and N concatenated NOTE: the procedure is called by intmatcoldelete and sortintmat" { if (nrows(M)>=nrows(N)) { int m=nrows(M); } else { int m=nrows(N); } intmat P[m][ncols(M)+ncols(N)]; P[1..nrows(M),1..ncols(M)]=M[1..nrows(M),1..ncols(M)]; P[1..nrows(N),ncols(M)+1..ncols(M)+ncols(N)]=N[1..nrows(N),1..ncols(N)]; return(P); } ///////////////////////////////////////////////////////////////////////// static proc minInIntvec (intvec v) "USAGE: minInIntvec(v); v intvec RETURN: list, first entry is the minimal value in v, second entry is its position in v NOTE: called by sortintmat" { int min=v[1]; int minpos=1; for (int i=2;i<=size(v);i++) { if (v[i]vglwert); } if (size(#)>0) { if (typeof(#[1])=="string") { ergebnis=1-(wertvglwert); } } setring BASERING; return(ergebnis); } ///////////////////////////////////////////////////////////////////////// static proc koeffizienten (poly f,int k) "USAGE: koeffizienten(f,k) f poly, k int ASSUME: f=a*x+b*y+c is a linear polynomial in two variables, k is either 0, 1 or 2 RETURN: poly, one of the coefficients of f, depending on the value of k: k=0 : c is returned k=1 : a is returned k=2 : b is returned NOTE: the procedure is called from tropicalCurve" { poly c=substitute(f,var(1),0,var(2),0); if (k==0) { return(c); } f=f-c; if (k==1) { return(substitute(f,var(1),1,var(2),0)); } else { return(substitute(f,var(1),0,var(2),1)); } } //////////////////////////////////////////////////////////////////////////// /// Procedures used in tex-procedures and conicWithTangents: /// - minOfPolys /// - shorten /// - minOfStringDecimal /// - decimal /// - stringdelete /// - stringinsert /// - texmonomial /// - texcoefficient /// - abs ///////////////////////////////////////////////////////////////////////////// static proc minOfPolys (list v) "USAGE: minOfPolys(v); v list RETURN: poly, the minimum of the numbers in v NOTE: called by texDrawTropical" { poly min=v[1]; for (int i=2;i<=size(v);i++) { if (v[i]cc) or (aa<-cc)) { ascale=abs(cc/aa); } if ((bb>cc) or (bb<-cc)) { bscale=abs(cc/bb); } if (ascale0) { nachkomma=#[1]; } while ((check==0) and i<=size(s)) { if (s[i]==".") { if (i+nachkomma>size(s)) { nachkomma=size(s)-i; } for (j=i;j<=i+nachkomma;j++) { news=news+s[j]; } check=1; } else { news=news+s[i]; } i++; } setring BASERING; return(news); } ///////////////////////////////////////////////////////////////////////// static proc stringcontainment (string a, string b) "USAGE: stringcontainment(a,b); a,b string ASSUME: a is a string containing a single letter RETURN: int, one if a is one of the letters in b and zero else NOTE: the procedure is called from extremeraysC" { int i; for (i=1;i<=size(b);i++) { if (a==b[i]) { return(1); } } return(0); } ///////////////////////////////////////////////////////////////////////// static proc stringdelete (string w,int i) "USAGE: stringdelete(w,i); w string, i int RETURN: string, the string w with the ith component deleted NOTE: the procedure is called by texNumber and choosegfanvector" { if ((i>size(w)) or (i<=0)) { return(w); } if ((size(w)==1) and (i==1)) { return(""); } if (i==1) { return(w[2..size(w)]); } if (i==size(w)) { return(w[1..size(w)-1]); } else { string erg=w[1..i-1],w[i+1..size(w)]; return(erg); } } ///////////////////////////////////////////////////////////////////////// static proc stringinsert (string w,string v,int i) "USAGE: stringinsert(w,v,i); w string, v string, i int RETURN: string, the string w where at the ith component v has been inserted NOTE: the procedure is called by texmonomial" { if (i==1) { return(v[1]+w); } if (i==size(w)+1) { return(w+v[1]); } else { string anfang=w[1..i-1]; string ende=w[i..size(w)]; string erg=anfang+v+ende; return(erg); } } ///////////////////////////////////////////////////////////////////////// static proc texmonomial (poly f) "USAGE: texmonomial(f); f poly RETURN: string, the tex-command for the leading monomial of f NOTE: the procedure is called by texPolynomial and texNumber" { int altshort=short; short=0; string F=string(leadmonom(f)); int k; for (int j=1;j<=size(F);j++) { if (F[j]=="^") { F=stringinsert(F,"{",j+1); k=j+2; while (k<=size(F) and ((F[k]=="0") or (F[k]=="1") or (F[k]=="2") or (F[k]=="3") or (F[k]=="4") or (F[k]=="5") or (F[k]=="6") or (F[k]=="7") or (F[k]=="8") or (F[k]=="9"))) { k++; } F=stringinsert(F,"}",k); j=k+1; } if (F[j]=="(") { F=stringdelete(F,j); F=stringinsert(F,"_{",j); k=j+2; while (k<=size(F) and ((F[k]=="0") or (F[k]=="1") or (F[k]=="2") or (F[k]=="3") or (F[k]=="4") or (F[k]=="5") or (F[k]=="6") or (F[k]=="7") or (F[k]=="8") or (F[k]=="9"))) { k++; } F=stringdelete(F,k); F=stringinsert(F,"}",k); j=k; } if (F[j]=="*") { F=stringdelete(F,j); j--; } } short=altshort; return(F); } ///////////////////////////////////////////////////////////////////////// static proc texcoefficient (poly f,list #) "USAGE: texcoefficient(f[,#]); f poly, # optional list RETURN: string, the tex command representing the leading coefficient of f using \frac as coefficient NOTE: if # is non-empty, then the cdot at the end is omitted; this is needed for the constant term" { string cdot; if (size(#)==0) { cdot="\\cdot "; } string koeffizient=texNumber(f); if (koeffizient=="1") { if (size(#)==0) { return(""); } else { return(koeffizient); } } if (koeffizient=="-1") { if (size(#)==0) { return("-"); } else { return(koeffizient); } } if (size(koeffizient)>5) { string tfractest=koeffizient[2..6]; if (tfractest=="tfrac") { return(koeffizient+cdot); } } int anzahlplus,anzahlminus; for(int j=1;j<=size(koeffizient);j++) { if (koeffizient[j]=="+") { anzahlplus++; } if (koeffizient[j]=="-") { anzahlminus++; } } if ((anzahlplus==0) and (anzahlminus==1) and (koeffizient[1]=="-")) { return(koeffizient+cdot); } else { if (anzahlplus+anzahlminus>=1) { return("("+koeffizient+")"+cdot); } else { return(koeffizient+cdot); } } } ///////////////////////////////////////////////////////////////////////// static proc abs (def n) "USAGE: abs(n); n poly or int RETURN: poly or int, the absolute value of n" { if (n>=0) { return(n); } else { return(-n); } } /////////////////////////////////////////////////////////////////////////////// /// Procedures only used to compute j-invariants /// - findNonLoopVertex /// - coordinatechange /// - weierstrassFormOfACubic /// - weierstrassFormOfA4x2Curve /// - weierstrassFormOfA2x2Curve /// - jInvariantOfACubic /// - jInvariantOfA4x2Curve /// - jInvariantOfA2x2Curve /// - jInvariantOfAPuiseuxCubic /////////////////////////////////////////////////////////////////////////////// static proc findNonLoopVertex (list graph) "USAGE: findNonLoopVertex(graph); graph list ASSUME: graph is a list as in the output of tropicalCurve RETURN: int, the number of a vertex which has only one edge NOTE: the procedure is called by tropicalJInvariant" { int i,j; for (i=1;i<=size(graph)-1;i++) { if ((ncols(graph[i][3])==1) and (graph[i][3][1,1]!=0)) { return(i); } } return(-1); } ///////////////////////////////////////////////////////////////////////// static proc coordinatechange (poly f,intmat A,intvec v) "USAGE: coordinatechange(f,A,v); f poly, A intmat, v intvec ASSUME: f is a polynomial in two variables, A is a 2x2 integer matrix, and v is an integer vector with 2 entries RETURN: poly, the polynomial transformed by (x,y)->A*(x,y)+v NOTE: the procedure is called by weierstrassForm" { poly g; int i; intmat exp[2][1]; for (i=1;i<=size(f);i++) { exp=leadexp(f[i]); exp=A*exp; g=g+leadcoef(f[i])*var(1)^(exp[1,1]+v[1])*var(2)^(exp[2,1]+v[2]); } return(g); } ///////////////////////////////////////////////////////////////////////// static proc weierstrassFormOfACubic (poly f,list #) "USAGE: weierstrassFormOfACubic(wf[,#]); wf poly, # list ASSUME: poly is a cubic polynomial RETURN: poly, the Weierstrass normal form of the cubic, 0 if poly is not a cubic NOTE: - the algorithm for the coefficients of the Weierstrass form is due to Fernando Rodriguez Villegas, villegas@math.utexas.edu - if an additional argument # is given, the simplified Weierstrass form is computed - the procedure is called by weierstrassForm - the characteristic of the base field should not be 2 or 3 if # is non-empty" { if (deg(f)!=3) { ERROR("The curve is not a cubic!"); } // store the coefficients of f in a specific order poly t0,s1,s0,r2,r1,r0,q3,q2,q1,q0=flatten(coeffs(f,ideal(var(2)^3,var(2)^2,var(2)^2*var(1),var(2),var(2)*var(1),var(2)*var(1)^2,1,var(1),var(1)^2,var(1)^3))); // test, if f is already in Weierstrass form if ((t0==0) and (s1==1) and (s0==0) and (r0==0) and (q0==-1) and (size(#)==0)) { return (f); } // compute the coefficients a1,a2,a3,a4, and a6 of the Weierstrass // form y^2+a1xy+a3y=x^3+a2x^2+a4x+a6 poly a1=r1; poly a2=-(s0*q2+s1*q1+r0*r2); poly a3=(9*t0*q0-s0*r0)*q3+((-t0*q1-s1*r0)*q2+(-s0*r2*q1-s1*r2*q0)); poly a4=((-3*t0*r0+s0^2)*q1+(-3*s1*s0*q0+s1*r0^2))*q3 +(t0*r0*q2^2+(s1*s0*q1+((-3*t0*r2+s1^2)*q0+s0*r0*r2))*q2 +(t0*r2*q1^2+s1*r0*r2*q1+s0*r2^2*q0)); poly a6=(-27*t0^2*q0^2+(9*t0*s0*r0-s0^3)*q0-t0*r0^3)*q3^2+ (((9*t0^2*q0-t0*s0*r0)*q1+((-3*t0*s0*r1+(3*t0*s1*r0+ 2*s1*s0^2))*q0+(t0*r0^2*r1-s1*s0*r0^2)))*q2+(-t0^2*q1^3 +(t0*s0*r1+(2*t0*s1*r0-s1*s0^2))*q1^2+((3*t0*s0*r2+ (-3*t0*s1*r1+2*s1^2*s0))*q0+((2*t0*r0^2-s0^2*r0)*r2+ (-t0*r0*r1^2+s1*s0*r0*r1-s1^2*r0^2)))*q1+((9*t0*s1*r2- s1^3)*q0^2+(((-3*t0*r0+s0^2)*r1-s1*s0*r0)*r2+(t0*r1^3 -s1*s0*r1^2+s1^2*r0*r1))*q0)))*q3+(-t0^2*q0*q2^3+ (-t0*s1*r0*q1+((2*t0*s0*r2+(t0*s1*r1-s1^2*s0))*q0- t0*r0^2*r2))*q2^2+(-t0*s0*r2*q1^2+(-t0*s1*r2*q0+ (t0*r0*r1-s1*s0*r0)*r2)*q1+((2*t0*r0-s0^2)*r2^2+ (-t0*r1^2+s1*s0*r1-s1^2*r0)*r2)*q0)*q2+ (-t0*r0*r2^2*q1^2+(t0*r1-s1*s0)*r2^2*q0*q1- t0*r2^3*q0^2)); poly b2=a1^2+4*a2; poly b4=2*a4+a1*a3; poly b6=a3^2+4*a6; poly b8=a1^2*a6+4*a2*a6-a1*a3*a4+a2*a3^2-a4^2; poly c4=b2^2-24*b4; poly c6=-b2^3+36*b2*b4-216*b6; if (size(#)!=0) { return(substitute(var(2)^2+a1*var(1)*var(2)+a3*var(2)-var(1)^3-a2*var(1)^2-a4*var(1)-a6,var(2),var(2)-(a1*var(1)+a3)/2)); } else { return(var(2)^2+a1*var(1)*var(2)+a3*var(2)-var(1)^3-a2*var(1)^2-a4*var(1)-a6); } } ///////////////////////////////////////////////////////////////////////// static proc weierstrassFormOfA4x2Curve (poly f) "USAGE: weierstrassFormOfA4x2Curve(wf); wf poly ASSUME: poly is a polynomial of type 4x2 RETURN: poly, the Weierstrass normal form of the curve NOTE: - the procedure is called by weierstrassForm - the characteristic of the base field should not be 2 or 3" { poly u11,u40,u30,u20,u10,u00,u21,u01,u02=flatten(coeffs(f,ideal(var(1)*var(2),var(1)^4,var(1)^3,var(1)^2,var(1),1,var(1)^2*var(2),var(2),var(2)^2))); poly c4=(u11^4 - 8*u11^2*u20*u02 - 8*u11^2*u01*u21 + 24*u11*u10*u21*u02 + 24*u11*u30*u01*u02 + 192*u00*u40*u02^2 - 48*u00*u21^2*u02 - 48*u10*u30*u02^2 + 16*u20^2*u02^2 - 16*u20*u01*u21*u02 - 48*u40*u01^2*u02 + 16*u01^2*u21^2)/lead(u02^4); poly c6=(-u11^6 + 12*u11^4*u20*u02 + 12*u11^4*u01*u21 - 36*u11^3*u10*u21*u02 - 36*u11^3*u30*u01*u02 + 576*u11^2*u00*u40*u02^2 + 72*u11^2*u00*u21^2*u02 + 72*u11^2*u10*u30*u02^2 - 48*u11^2*u20^2*u02^2 - 24*u11^2*u20*u01*u21*u02 + 72*u11^2*u40*u01^2*u02 - 48*u11^2*u01^2*u21^2 - 864*u11*u00*u30*u21*u02^2 + 144*u11*u10*u20*u21*u02^2 - 864*u11*u10*u40*u01*u02^2 + 144*u11*u10*u01*u21^2*u02 + 144*u11*u20*u30*u01*u02^2 + 144*u11*u30*u01^2*u21*u02 - 2304*u00*u20*u40*u02^3 + 576*u00*u20*u21^2*u02^2 + 864*u00*u30^2*u02^3 + 1152*u00*u40*u01*u21*u02^2 - 288*u00*u01*u21^3*u02 + 864*u10^2*u40*u02^3 - 216*u10^2*u21^2*u02^2 - 288*u10*u20*u30*u02^3 + 144*u10*u30*u01*u21*u02^2 + 64*u20^3*u02^3 - 96*u20^2*u01*u21*u02^2 + 576*u20*u40*u01^2*u02^2 - 96*u20*u01^2*u21^2*u02 - 216*u30^2*u01^2*u02^2 - 288*u40*u01^3*u21*u02 + 64*u01^3*u21^3)/lead(u02^6); return(var(2)^2-var(1)^3+27*c4*var(1)+54*c6); } ///////////////////////////////////////////////////////////////////////// static proc weierstrassFormOfA2x2Curve (poly f) "USAGE: weierstrassFormOfA2x2Curve(f); f poly ASSUME: poly, is a polynomial defining an elliptic curve of type (2,2) on P^1xP^1 i.e. a polynomial of the form a+bx+cx2+dy+exy+fx2y+gy2+hxy2+ix2y2 RETURN: poly, a Weierstrass form of the elliptic curve defined by poly NOTE: - the algorithm is based on the paper Sang Yook An, Seog Young Kim, David C. Marshall, Susan H. Marshall, William G. McCallum and Alexander R. Perlis: Jacobians of Genus One Curves. Journal of Number Theory 90,2 (2001), 304-315. - the procedure is called by weierstrassForm - the characteristic of the base field should not be 2 or 3" { // get the coefficients of the polynomial poly A00,A10,A20,A01,A11,A21,A02,A12,A22=flatten(coeffs(f,ideal(1,var(1),var(1)^2,var(2),var(1)*var(2),var(1)^2*var(2),var(2)^2,var(1)*var(2)^2,var(1)^2*var(2)^2))); // define P1xP1 as quadric in P3 matrix A[4][4]=0,0,0,1/2,0,0,1/2,0,0,1/2,0,0,1/2,0,0,0; // define the image of the (2,2)-curve under the Segre embedding // as quadric which should be // intersected with the image of P1xP1 matrix B[4][4]=A00,A10/2,A01/2,0,A10/2,A20,A11/2,A21/2,A01/2,A11/2,A02,A12/2,0,A21/2,A12/2,A22; // compute the coefficients of the Weierstrass form of // the input polynomial and its j-invariant poly a=det(x*A+B); poly a0,a1,a2,a3,a4=flatten(coeffs(a,ideal(x4,x3,x2,x,1))); a1,a2,a3=-a1/4,a2/6,-a3/4; poly g2=a0*a4-4*a1*a3+3*a2^2; poly g3=a0*a2*a4+2*a1*a2*a3-a0*a3^2-a4*a1^2-a2^3; return(var(2)^2-var(1)^3+g2/4*var(1)+g3/4); } ///////////////////////////////////////////////////////////////////////// static proc jInvariantOfACubic (poly f,list #) "USAGE: jInvariantOfACubic(f[,#]); f poly, # list ASSUME: poly is a cubic polynomial defining an elliptic curve RETURN: poly, the j-invariant of the elliptic curve defined by poly NOTE: - if the base field is Q(t) an optional argument # may be given; then the algorithm only computes the negative of the order of the j-invariant" { if (deg(f)!=3) { ERROR("The input polynomial is not a cubic!"); } // compute first the Weierstrass form of the cubic // - this does not change the j-invariant f=weierstrassFormOfACubic(f); // compute the coefficients of the Weierstrass form poly a1,a2,a3,a4,a6=flatten(coeffs(f,ideal(var(1)*var(2),var(1)^2,var(2),var(1),1))); a2,a4,a6=-a2,-a4,-a6; // compute the j-invariant poly b2=a1^2+4*a2; poly b4=2*a4+a1*a3; poly b6=a3^2+4*a6; poly b8=a1^2*a6+4*a2*a6-a1*a3*a4+a2*a3^2-a4^2; poly c4=b2^2-24*b4; // poly c6=-b2^3+36*b2*b4-216*b6; poly delta=-b2^2*b8-8*b4^3-27*b6^2+9*b2*b4*b6; if (delta==0) { ERROR("The input is a rational curve and has no j-invariant!"); } if (size(#)>0) // if the optional argument is given, then only the { // negative of the order is computed int zaehler=3*simplifyToOrder(c4)[1]; int nenner=simplifyToOrder(delta)[1]; return(nenner-zaehler); } else { poly invariant=(c4^3/delta); return(invariant); } } ///////////////////////////////////////////////////////////////////////// static proc jInvariantOfA2x2Curve (poly f,list #) "USAGE: jInvariantOfA2x2Curve(f[,#]); f poly, # list ASSUME: poly, is a polynomial defining an elliptic curve of type (2,2) on P^1xP^1 i.e. a polynomial of the form a+bx+cx2+dy+exy+fx2y+gy2+hxy2+ix2y2 RETURN: poly, the j-invariant of the elliptic curve defined by poly NOTE: - if the base field is Q(t) an optional argument # may be given; then the algorithm only computes the negative of the order of the j-invariant - the characteristic should not be 2 or 3 - the procedure is not called at all, it is just stored for the purpose of comparison" { // get the coefficients of the polynomial poly A00,A10,A20,A01,A11,A21,A02,A12,A22=flatten(coeffs(f,ideal(1,var(1),var(1)^2,var(2),var(1)*var(2),var(1)^2*var(2),var(2)^2,var(1)*var(2)^2,var(1)^2*var(2)^2))); // define P1xP1 as quadric in P3 matrix A[4][4]=0,0,0,1/2,0,0,1/2,0,0,1/2,0,0,1/2,0,0,0; // define the image of the (2,2)-curve under the Segre embedding as quadric which should be // intersected with the image of P1xP1 matrix B[4][4]=A00,A10/2,A01/2,0,A10/2,A20,A11/2,A21/2,A01/2,A11/2,A02,A12/2,0,A21/2,A12/2,A22; // compute the coefficients of the Weierstrass form of // the input polynomial and its j-invariant poly a=det(var(1)*A+B); poly a0,a1,a2,a3,a4=flatten(coeffs(a,ideal(var(1)^4,var(1)^3,var(1)^2,var(1),1))); a1,a2,a3=-a1/4,a2/6,-a3/4; poly g2=a0*a4-4*a1*a3+3*a2^2; poly g3=a0*a2*a4+2*a1*a2*a3-a0*a3^2-a4*a1^2-a2^3; poly G2cube=g2^3; poly G3square=g3^2; poly jinvnum=1728*G2cube; poly jinvdenom=-27*G3square+G2cube; // check if the curve is rational if (jinvdenom==0) { ERROR("The input is a rational curve and has no j-invariant!"); } if (size(#)>0) // if the optional argument is given, { // then only the negative of the order is computed int zaehler=simplifyToOrder(jinvnum)[1]; int nenner=simplifyToOrder(jinvdenom)[1]; return(nenner-zaehler); } else { poly invariant=(jinvnum/jinvdenom); return(invariant); } } ///////////////////////////////////////////////////////////////////////// static proc jInvariantOfA4x2Curve (poly f,list #) "USAGE: jInvariantOfA4x2Curve(f[,#]); f poly, # list ASSUME: poly, is a polynomial defining an elliptic curve of type (4,2), i.e. a polynomial of the form a+bx+cx2+dx3+ex4+fy+gxy+hx2y+iy2 RETURN: poly, the j-invariant of the elliptic curve defined by poly NOTE: - if the base field is Q(t) an optional argument # may be given; then the algorithm only computes the negative of the order of the j-invariant - the characteristic should not be 2 or 3 - the procedure is not called at all, it is just stored for the purpose of comparison" { poly u11,u40,u30,u20,u10,u00,u21,u01,u02=flatten(coeffs(f,ideal(var(1)*var(2),var(1)^4,var(1)^3,var(1)^2,var(1),1,var(1)^2*var(2),var(2),var(2)^2))); poly c4=(u11^4 - 8*u11^2*u20*u02 - 8*u11^2*u01*u21 + 24*u11*u10*u21*u02 + 24*u11*u30*u01*u02 + 192*u00*u40*u02^2 - 48*u00*u21^2*u02 - 48*u10*u30*u02^2 + 16*u20^2*u02^2 - 16*u20*u01*u21*u02 - 48*u40*u01^2*u02 + 16*u01^2*u21^2); poly c6=(-u11^6 + 12*u11^4*u20*u02 + 12*u11^4*u01*u21 - 36*u11^3*u10*u21*u02 - 36*u11^3*u30*u01*u02 + 576*u11^2*u00*u40*u02^2 + 72*u11^2*u00*u21^2*u02 + 72*u11^2*u10*u30*u02^2 - 48*u11^2*u20^2*u02^2 - 24*u11^2*u20*u01*u21*u02 + 72*u11^2*u40*u01^2*u02 - 48*u11^2*u01^2*u21^2 - 864*u11*u00*u30*u21*u02^2 + 144*u11*u10*u20*u21*u02^2 - 864*u11*u10*u40*u01*u02^2 + 144*u11*u10*u01*u21^2*u02 + 144*u11*u20*u30*u01*u02^2 + 144*u11*u30*u01^2*u21*u02 - 2304*u00*u20*u40*u02^3 + 576*u00*u20*u21^2*u02^2 + 864*u00*u30^2*u02^3 + 1152*u00*u40*u01*u21*u02^2 - 288*u00*u01*u21^3*u02 + 864*u10^2*u40*u02^3 - 216*u10^2*u21^2*u02^2 - 288*u10*u20*u30*u02^3 + 144*u10*u30*u01*u21*u02^2 + 64*u20^3*u02^3 - 96*u20^2*u01*u21*u02^2 + 576*u20*u40*u01^2*u02^2 - 96*u20*u01^2*u21^2*u02 - 216*u30^2*u01^2*u02^2 - 288*u40*u01^3*u21*u02 + 64*u01^3*u21^3); poly c4cube=c4^3; poly jdenom=c4cube-c6^2; if (size(#)>0) // if the optional argument is given, { // then only the negative of the order is computed int zaehler=3*simplifyToOrder(c4)[1]; int nenner=simplifyToOrder(jinvdenom)[1]; return(nenner-zaehler); } else { return(1728*(c4cube/(jdenom))); } } ///////////////////////////////////////////////////////////////////////// static proc jInvariantOfAPuiseuxCubic (poly f,list #) "USAGE: jInvariantOfAPuiseuxCubic(f[,#]); f poly, # list ASSUME: poly is a cubic polynomial over Q(t) defining an elliptic curve and # is empty RETURN: list, containing two polynomials whose ratio is the j-invariant of the elliptic curve defined by poly ASSUME: poly is a cubic polynomial over Q(t) defining an elliptic curve and # is non-empty RETURN: int, the order of the j-invariant of the elliptic curve defined by poly NOTE: the procedure is called by jInvariant" { if (deg(f)!=3) { ERROR("The input polynomial is not a cubic!"); } // compute first the Weierstrass form of the cubic // - this does not change the j-invariant f=weierstrassFormOfACubic(f); // compute the coefficients of the Weierstrass form poly a1,a2,a3,a4,a6=flatten(coeffs(f,ideal(var(1)*var(2),var(1)^2,var(2),var(1),1))); a2,a4,a6=-a2,-a4,-a6; number dna1,dna2,dna3,dna4,dna6=denominator(leadcoef(a1)),denominator(leadcoef(a2)),denominator(leadcoef(a3)),denominator(leadcoef(a4)),denominator(leadcoef(a6)); number na1,na2,na3,na4,na6=numerator(leadcoef(a1)),numerator(leadcoef(a2)),numerator(leadcoef(a3)),numerator(leadcoef(a4)),numerator(leadcoef(a6)); number hn=dna1*dna2*dna3*dna4*dna6; a1=na1*dna2*dna3*dna4*dna6; a2=dna1*na2*dna3*dna4*dna6; a3=dna1*dna2*na3*dna4*dna6; a4=dna1*dna2*dna3*na4*dna6; a6=dna1*dna2*dna3*dna4*na6; // compute the j-invariant poly b2=a1^2+4*a2*hn; poly b4=2*a4*hn+a1*a3; poly b6=a3^2+4*a6*hn; poly b8=a1^2*a6+4*a2*a6*hn-a1*a3*a4+a2*a3^2-a4^2*hn; poly c4=b2^2-24*b4*hn^2; poly delta=-b2^2*b8-8*b4^3*hn-27*b6^2*hn^3+9*b2*b4*b6*hn; if (delta==0) { ERROR("The input is a rational curve and has no j-invariant!"); } if (size(#)>0) // if the optional argument is given, { // then only the negative of the order is computed int zaehler=3*simplifyToOrder(c4)[1]; int nenner=simplifyToOrder(delta)[1]; int ordhn=simplifyToOrder(hn)[1]; return(zaehler-nenner-5*ordhn); } else { def BASERING=basering; execute("ring TRING="+string(char(BASERING))+",t,ds;"); poly hn=imap(BASERING,hn); poly c4=imap(BASERING,c4); poly delta=imap(BASERING,delta); poly num=c4^3; poly denom=delta*hn^5; poly ggt=gcd(num,denom); num=num/ggt; denom=denom/ggt; setring BASERING; poly num=imap(TRING,num); poly denom=imap(TRING,denom); number teiler=1; if (char(BASERING)==0) { teiler=gcd(leadcoef(num),leadcoef(denom)); } num=num/teiler; denom=denom/teiler; list invariant=num,denom; return(invariant); } } ///////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////// /// Further examples for testing the procedures ///////////////////////////////////////////////////////////////////////////// /* /// Note, only the following procedures need further examples /// (the others are too simple): /// A) tropicalLifting (best tested via displayTropicalLifting) /// B) tropicalCurve (best tested via drawTropicalCurve) /// C) tropicalJInvariant /// D) weierstrassForm /// E) jInvariant //////////////////////////////////////////////////////////////////////////// /// A) Examples tropicalLifting //////////////////////////////////////////////////////////////////////////// // ------------------------------------------------------- // Example 1 // ------------------------------------------------------- ring r0=0,(t,x),dp; poly f=t7-6t4+3t3x+8t3-12t2x+6tx2-x3+t2; f; // The point -2/3 is in the tropical variety. list L=tropicalLifting(f,intvec(3,-2),4); L; displayTropicalLifting(L,"subst"); // -------------------------------------------------------- // Example 2 - a field extension is necessary // -------------------------------------------------------- poly g=(1+t2)*x2+t5x+t2; g; // The poin -1 is in the tropical variety. displayTropicalLifting(tropicalLifting(g,intvec(1,-1),4),"subst"); // -------------------------------------------------------- // Example 3 - the ideal is not zero dimensional // -------------------------------------------------------- ring r1=0,(t,x,y),dp; poly f=(9t27-12t26-5t25+21t24+35t23-51t22-40t21+87t20+56t19-94t18-62t17+92t16+56t15-70t14-42t13+38t12+28t11+18t10-50t9-48t8+40t7+36t6-16t5-12t4+4t2)*x2+(-9t24+12t23-4t22-42t20+28t19+30t18-20t17-54t16+16t15+48t14-16t12+8t11-18t10-26t9+30t8+20t7-24t6+4t5+28t4-8t3-16t2+4)*xy+(6t16-10t15-2t14+16t13+4t12-18t11-10t10+24t9+6t8-20t7+8t6+8t5-20t4-4t3+12t2+4t-4)*y2+(-9t28+3t27+8t26-4t25-45t24-6t23+42t22+30t21-94t20-40t19+68t18+82t17-38t16-60t15+26t14+36t13-22t12-20t11+4t10+4t9+12t8+8t7-8t6-8t5+4t4)*x+(9t27-21t26+16t25+14t24+12t23-61t22+27t21+80t20-19t19-100t18+26t17+96t16-24t15-84t14+44t12-2t11-18t10+2t9+40t8+4t7-32t6-12t5+4t3+12t2-4)*y+(9t27-12t26+4t25+36t23-18t22-28t21-2t20+54t19+14t18-52t17-44t16+24t15+40t14-4t13-12t12+4t11-4t10-4t9+4t8); f; displayTropicalLifting(tropicalLifting(f,intvec(1,-1,-4),3),"subst"); // -------------------------------------------------------- // Example 4 - the ideal has even more equations // -------------------------------------------------------- ring r2=0,(t,x,y,z),dp; ideal i=t-x3+3yz,t2xy-2x2z; i; displayTropicalLifting(tropicalLifting(i,intvec(1,1,3,0),2),"subst"); // -------------------------------------------------------- // Example 5-7 - testing some options // -------------------------------------------------------- setring r0; poly f1=(x2-t3)*(x3-t5)*(x5-t7)*(x7-t11)*(x11-t13); f1; displayTropicalLifting(tropicalLifting(f1,intvec(7,-11),4,"noAbs"),"subst"); poly f2=(1+t2)*x2+t5x+t2; f2; displayTropicalLifting(tropicalLifting(f2,intvec(1,-1),4,"isZeroDimensional","findAll"),"subst"); poly f3=t7-6t4+3t3x+8t3-12t2x+6tx2-x3+t2; f3; displayTropicalLifting(tropicalLifting(f3,intvec(3,-2),4,"isInTrop","findAll"),"subst"); // --------------------------------------------------------- // Even more examples // --------------------------------------------------------- ring r1=0,(t,x),dp; poly f1=(x2-t3)*(x3-t5)*(x5-t7)*(x7-t11)*(x11-t13); displayTropicalLifting(tropicalLifting(f1,intvec(7,-11),4,"noAbs"),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (7,-11) /////////////////////////////////////////////////////////////////// // one field extension needed // x=(X(1))*t+(-1/2*X(1))*t^3+(3/8*X(1)-1/2)*t^5+(-5/16*X(1)+1/2)*t^7+(19/128*X(1)-1/2)*t^9+(-15/256*X(1)+1/2)*t^11+(-9/1024*X(1)-1/2)*t^13 where t->t and X(1)^2+1=0 poly f2=(1+t2)*x2+t5x+t2; displayTropicalLifting(tropicalLifting(f2,intvec(1,-1),4,"isZeroDimensional","findAll"),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (1,-1) /////////////////////////////////////////////////////////////////// // x=t2+2t3+t7 the result is exact poly f3=t7-6t4+3t3x+8t3-12t2x+6tx2-x3+t2; displayTropicalLifting(tropicalLifting(f3,intvec(3,-2),4,"isInTrop","findAll"),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (3,-2) /////////////////////////////////////////////////////////////////// // full parametrisation is given by t->t2+t3 and x=3t7+t9+4t11-t12, this however cannot be computed // so: t->t2 and x=-3t7-21/2t8-239/8t9-78t10-25589/128t11-506t12-1298967/1024t13-3159t14-256667105/32768t15-19371t16-12540259065/262144t17-118066t18-1222177873545/4194304t19-719423t20-59637294735495/33554432t21 poly f4=t12-83t11-88t10-69t9+3t8x+153t8+298t7x-81t7+165t6x-402t5x+3t4x2+189t4x-131t3x2+213t2x2-86tx2+x3+9x2; displayTropicalLifting(tropicalLifting(f4,intvec(2,-7),4,"findAll"),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (2,-7) /////////////////////////////////////////////////////////////////// // t->t6 and x=3t3+t7+4t11-t12 poly f5=t12-4120t11+6t10x-1008t10-96t9x+15t8x2+4453t9-2016t8x-144t7x2+20t6x3-18t8-108t7x-1008t6x2-96t5x3+15t4x4-4681t7-36t6x-162t5x2-24t3x4+6t2x5+243t6-18t4x2-108t3x3+x6+1890t5+486t4x-27tx4+243t2x2-729t3; displayTropicalLifting(tropicalLifting(f5,intvec(2,-1),3,"isInTrop","isZeroDimensional","findAll"),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (2,-1) /////////////////////////////////////////////////////////////////// ring r2=0,(t,x,y),dp; poly f7=(9t27-12t26-5t25+21t24+35t23-51t22-40t21+87t20+56t19-94t18-62t17+92t16+56t15-70t14-42t13+38t12+28t11+18t10-50t9-48t8+40t7+36t6-16t5-12t4+4t2)*x2+(-9t24+12t23-4t22-42t20+28t19+30t18-20t17-54t16+16t15+48t14-16t12+8t11-18t10-26t9+30t8+20t7-24t6+4t5+28t4-8t3-16t2+4)*xy+(6t16-10t15-2t14+16t13+4t12-18t11-10t10+24t9+6t8-20t7+8t6+8t5-20t4-4t3+12t2+4t-4)*y2+(-9t28+3t27+8t26-4t25-45t24-6t23+42t22+30t21-94t20-40t19+68t18+82t17-38t16-60t15+26t14+36t13-22t12-20t11+4t10+4t9+12t8+8t7-8t6-8t5+4t4)*x+(9t27-21t26+16t25+14t24+12t23-61t22+27t21+80t20-19t19-100t18+26t17+96t16-24t15-84t14+44t12-2t11-18t10+2t9+40t8+4t7-32t6-12t5+4t3+12t2-4)*y+(9t27-12t26+4t25+36t23-18t22-28t21-2t20+54t19+14t18-52t17-44t16+24t15+40t14-4t13-12t12+4t11-4t10-4t9+4t8); displayTropicalLifting(tropicalLifting(f7,intvec(1,-1,-4),4,"isPrime"),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (1,-1,-4) /////////////////////////////////////////////////////////////////// ring r3=0,(t,x,y,z),dp; ideal i1=-y2t4+x2,yt3+xz+y; displayTropicalLifting(tropicalLifting(i1,intvec(1,-2,0,2),4,"findAll"),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (1,-2,0,2) /////////////////////////////////////////////////////////////////// ideal i2=-y2t4+x2,yt3+yt2+xz; displayTropicalLifting(tropicalLifting(i2,intvec(1,-3,-1,0),4),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (1,-3,-1,0) //////////////////////////////////////////////////////////////// ring r4=0,(t,x,y,z,u),dp; ideal i3=t3x+ty+u,tx+t2u+y; displayTropicalLifting(tropicalLifting(i3,intvec(1,-1,-2,-1,-3),6),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (1,-1,-2,-1,-3) ///////////////////////////////////////////////////////////////////// setring r2; ideal i5=t3-t2y+t2-tx-ty+xy,t3-t2y+t2-2ty+y2,-t4+2t3x-t2x2-t3+3t2x-3tx2+x3+t2+t-x,-t4+2t3x-t2x2-t3+2t2x-tx2+t2y-2txy+x2y+t2+t-y; displayTropicalLifting(tropicalLifting(i5,intvec(1,-1,-1),3),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (1,-1,-1) //////////////////////////////////////////////////////////////////// ring r5=0,(t,V,W,X,Y,Z),dp; ideal i5=VXt5+3VZt4-2VXt3+WXt3+V2t2+3WZt2-2WXt+VW,-4W2X4+V2t2; displayTropicalLifting(tropicalLifting(i5,intvec(2,-3,-3,-1,-1,-1),3),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (2,-3,-3,-1,-1,-1) /////////////////////////////////////////////////////////////////// ring r6=0,(t,x,y,z),dp; ideal i6=t4xy-y3z3+3t2xy3-txy2+t4y; displayTropicalLifting(tropicalLifting(i6,intvec(3,0,-9,2),4),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (3,0,-9,2) //////////////////////////////////////////////////////////// // very easy example stops after one loop with an exact result ideal i7=t2xy-2x2z; displayTropicalLifting(tropicalLifting(i7,intvec(1,1,3,0),4),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (1,1,3,0) ///////////////////////////////////////////////////////////// displayTropicalLifting(tropicalLifting(i7,intvec(1,-3,-1,0),4),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (1,-3,-1,0) ///////////////////////////////////////////////////////////////// // requires one field extension ideal i8=t-x3+3yz; displayTropicalLifting(tropicalLifting(i8,intvec(1,1,3,0),4),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (1,1,3,0) ///////////////////////////////////////////////////////////////// displayTropicalLifting(tropicalLifting(i8,intvec(1,-3,-1,0),4),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (1,-3,-1,0) ///////////////////////////////////////////////////////////////// // requires one field extension ideal i9=t-x3+3yz,t2xy-2x2z; displayTropicalLifting(tropicalLifting(i9,intvec(1,1,3,0),4),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (1,1,3,0) // takes too much time //displayTropicalLifting(tropicalLifting(i9,intvec(1,-3,-1,0),4),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (1,-3,-1,0) ///////////////////////////////////////////////////////////////////// // too much for Gfan did not work on my computer //ideal i10=t2xy+t4-z2,t-2xy+z2+tz,z2-y; //displayTropicalLifting(tropicalLifting(i10,intvec(2,4,-6,-3),2),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (2,4,-6,-3) //////////////////////////////////////////////////////////////////// ring r7=0,(t,x,y),dp; ideal i11=t4xy-y3+3t2xy3-txy2+t4y; displayTropicalLifting(tropicalLifting(i11,intvec(1,3,1),5),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (1,3,1) //////////////////////////////////////////////////////////////// ring r8=0,(t,x),dp; ideal i12=x32+(16t6)*x30+(8t10)*x29+(120t12)*x28+(112t16)*x27+(32t20+560t18)*x26+(728t22)*x25+(388t26+1820t24)*x24+(80t30+2912t28)*x23+(2160t32+4368t30)*x22+(824t36+8008t34)*x21+(138t40+7304t38+8008t36)*x20+(3840t42+16016t40)*x19+(1180t46+16720t44+11440t42)*x18+(170t50+10680t48+24024t46)*x17+(4480t52+27324t50+12870t48)*x16+(1168t56+19680t54+27456t52)*x15+(152t60+9920t58+32736t56+11440t54)*x14+(3472t62+25200t60+24024t58)*x13+(804t66+14140t64+29040t62+8008t60)*x12+(96t70+5824t68+22848t66+16016t64)*x11+(1776t72+13496t70+19008t68+4368t66)*x10+(364t76+6020t74+14640t72+8008t70)*x9+(42t80+2112t78+8680t76+9020t74+1820t72)*x8+(544t82+3920t80+6480t78+2912t76)*x7+(104t86+1448t84+3680t82+2992t80+560t78)*x6+(12t90+400t88+1568t86+1880t84+728t82)*x5+(94t92+564t90+970t88+648t86+120t84)*x4+(16t96+144t94+352t92+320t90+112t88)*x3+(3t100+36t98+112t96+140t94+80t92+16t90)*x2+(6t102+20t100+34t98+24t96+8t94)*x+(t106+5t104+8t102+8t100+4t98+t96); displayTropicalLifting(tropicalLifting(i12,intvec(1,-3),2),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (1,-3) //////////////////////////////////////////////////////////////// ring r9=3,(t,x),dp; ideal i13=imap(r8,i12); displayTropicalLifting(tropicalLifting(i13,intvec(1,-3),3,"noAbs"),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (1,-3) ///////////////////////////////////////////////////// // without the option findAll we first get only a non-valid solution !!! ring r10=0,(t,x,y,z,u),dp; ideal i14=-y2t5-xyt4-2yzt3+yut2+xut+2zu; displayTropicalLifting(tropicalLifting(i14,intvec(1,-2,-1,-3,-3),3),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (1,-2,-1,-3,-3) /////////////////////////////////////////////////////// // test for positive characteristic /////////////////////////////////////////////////////////////////// ring r11=13,(t,x),dp; poly f1=(x2-t3)*(x3-t5)*(x5-t7)*(x7-t11)*(x11-t13); displayTropicalLifting(tropicalLifting(f1,intvec(7,-11),4,"findAll"),"subst"); /////////////////////////////////////////////////////////////////// // weight was: (7,-11) /////////////////////////////////////////////////////////////////// // one field extension needed // x=(X(1))*t+(-1/2*X(1))*t^3+(3/8*X(1)-1/2)*t^5+(-5/16*X(1)+1/2)*t^7+(19/128*X(1)-1/2)*t^9+(-15/256*X(1)+1/2)*t^11+(-9/1024*X(1)-1/2)*t^13 where t->t and X(1)^2+1=0 poly f2=(1+t2)*x2+t5x+t2; displayTropicalLifting(tropicalLifting(f2,intvec(1,-1),4,"isZeroDimensional","findAll"),"subst"); /////////////////////////////////////////////////////////////////// ring r12=5,(t,x,y),dp; ideal i11=t4xy-y3+3t2xy3-txy2+t4y; displayTropicalLifting(tropicalLifting(i11,intvec(1,3,1),5),"subst"); /////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////// /// B) Examples for tropicalCurve /////////////////////////////////////////////////////////////////// ring r=(0,t),(x,y),dp; poly f=t*(x7+y7+1)+1/t*(x4+y4+x2+y2+x3y+xy3)+1/t7*x2y2; list graph=tropicalCurve(f); graph; size(graph)-1; drawTropicalCurve(graph); poly g=t3*(x7+y7+1)+1/t3*(x4+y4+x2+y2+x3y+xy3)+1/t21*x2y2; list tropical_g=tropicalise(g); tropical_g; drawTropicalCurve(tropical_g); /// further examples poly f1=x+y; poly f2=xy+xy3+xy6+txy9; poly f3=xy+x2y3+x3y6+tx4y9; poly f4=t*(x3+y3+1)+1/t*(x2+y2+x+y+x2y+xy2)+1/t2*xy; poly f5=t8x3+t8y3+t8+t4x2y+t4xy2+t4x2+t4y2+t4x+t4y+t2xy; poly f6=x3+y3+1; poly f7=t*x70+t*y70+1/t*x40+1/t7*x20y20+1/t*y40+1/t*x30y+1/t*xy30+1/t*x20+1/t*y20+t; poly f8=t10*x7+t10*y7+1/t10*x4+1/t70*x2y2+1/t10*y4+1/t10*x3y+1/t10*xy3+1/t10*x2+1/t10*y2+t10; poly f9=x+y+x2y+xy2+1/t*xy; poly f10=1+x+y+xy+tx2y+txy2+t3x5+t3y5+tx2y2+t2xy4+t2yx4; poly f11=1/(t8)*x6y+1/(t16)*x5y2+1/(t9)*x4y3+(t4)*x2y5+1/(t11)*xy6+(t5)*y7+1/(t5)*x6+(t16)*x4y2+(t10)*x3y3+(t19)*xy5+(t14)*x5+1/(t17)*x3y2+(t10)*x2y3+1/(t20)*xy4+(t13)*y5+1/(t17)*x4+(t9)*x3y+(t11)*xy3+1/(t3)*xy2+(t7)*y3+1/(t18)*x2+(t16)*x; poly f12=-x3+(6t13+3t12+t5+3t4+3t3+1)/(t19)*x2+(2t5+1)/(t12)*xy+y2+(9t28-12t23-12t22+3t21+9t19+3t16-8t15-16t14-24t13-9t12-6t11-8t10-3t9-3t8-4t6-8t5-8t4-3t3-3t-2)/(t35)*x+(9t27+6t20-8t14-4t13-6t12-3t11-4t9-t8-4t5-4t4-2t3-1)/(t27)*y+(27t50+27t43-36t37-9t36-27t35-27t34-t32+12t31-18t30-11t29-54t28-36t27-12t25+13t24+28t23+61t22+27t21+18t20+16t19+3t18+14t17+3t16+16t15+26t14+44t13+34t12+21t11+27t10+12t9+7t8+t6+10t5+11t4+7t3+t2+2t+2)/(t50); poly f13=x3+(3t9+1)/(t8)*x2y+(3t10+2t+1)/(t8)*xy2+(t10+t7+t+1)/(t7)*y3+1/(t8)*x2+(2t5+1)/(t12)*xy+(t5+t3+1)/(t11)*y2+1/(t8)*x+(t+1)/(t8)*y+1; poly f14=x3+(3t9+3t6c+1)/(t8)*x2y+(3t12+6t9c+3t6c2+2t3+t2+2c)/(t10)*xy2+(t15+3t12c+t12+3t9c2+t6c3+t6+t5+2t3c+t2c+c2)/(t12)*y3+1/(t8)*x2+(2t5+2t2c+1)/(t12)*xy+(t8+t6+2t5c+t3+t2c2+c)/(t14)*y2+1/(t8)*x+(t3+t2+c)/(t10)*y+1; poly f15=x3+x2y+xy2+(t5)*y3+x2+1/(t4)*xy+y2+x+y+1; poly f16=t*(x7+y7+1)+1/t*(x4+y4+x2+y2+x3y+xy3)+1/t7*x2y2; /////////////////////////////////////////////////////////////////// /// B) Examples for tropicalJInvariant /////////////////////////////////////////////////////////////////// // tropcial_j_invariant computes the tropical j-invariant of the elliptic curve f tropicalJInvariant(t*(x3+y3+1)+1/t*(x2+y2+x+y+x2y+xy2)+1/t2*xy); // the Newton polygone need not be the standard simplex tropicalJInvariant(x+y+x2y+xy2+1/t*xy); // the curve can have arbitrary degree tropicalJInvariant(t*(x7+y7+1)+1/t*(x4+y4+x2+y2+x3y+xy3)+1/t7*x2y2); // the procedure does not realise, if the embedded graph of the tropical curve has // a loop that can be resolved tropicalJInvariant(1+x+y+xy+tx2y+txy2); // but it does realise, if the curve has no loop at all ... tropicalJInvariant(x+y+1); // or if the embedded graph has more than one loop - even if only one cannot be resolved tropicalJInvariant(1+x+y+xy+tx2y+txy2+t3x5+t3y5+tx2y2+t2xy4+t2yx4); // f is already in Weierstrass form weierstrassForm(y2+yx+3y-x3-2x2-4x-6); // g is not, but wg is g=x+y+x2y+xy2+1/t*xy; poly wg=weierstrassForm(g); wg; // ... but it is not yet a simple, since it still has an xy-term, unlike swg poly swg=weierstrassForm(g,1); swg; // the j-invariants of all three polynomials coincide ... jInvariant(g); jInvariant(wg); jInvariant(swg); // the following curve is elliptic as well poly h=x22y11+x19y10+x17y9+x16y9+x12y7+x9y6+x7y5+x2y3; // its Weierstrass form is weierstrassForm(h); jInvariant(x+y+x2y+y3+1/t*xy,"ord"); // see also example f5-f16 in part B) */ ///////////////////////////////////////////////////////////////////////// proc drawTwoTropicalCurves (list ff,list #) "USAGE: drawTropicalCurve(f[,#]); f poly or list, # optional list ASSUME: f is list of linear polynomials of the form ax+by+c with integers a, b and a rational number c representing a tropical Laurent polynomial defining a tropical plane curve; alternatively f can be a polynomial in Q(t)[x,y] defining a tropical plane curve via the valuation map; the basering must have a global monomial ordering, two variables and up to one parameter! RETURN: NONE NOTE: - the procedure creates the files /tmp/tropicalcurveNUMBER.tex and /tmp/tropicalcurveNUMBER.ps, where NUMBER is a random four digit integer; moreover it displays the tropical curve via kghostview; if you wish to remove all these files from /tmp, call the procedure cleanTmp @* - edges with multiplicity greater than one carry this multiplicity @* - if # is empty, then the tropical curve is computed w.r.t. minimum, if #[1] is the string 'max', then it is computed w.r.t. maximum @* - if the last optional argument is 'onlytexfile' then only the latex file is produced; this option should be used if kghostview is not installed on your system @* - note that lattice points in the Newton subdivision which are black correspond to markings of the marked subdivision, while lattice points in grey are not marked EXAMPLE: example drawTropicalCurve shows an example" { // check if the option "onlytexfile" is set, then only a tex file is produced if (size(#)!=0) { if (#[size(#)]=="onlytexfile") { int onlytexfile; #=delete(#,size(#)); } } // store this list # for use in the Newton subdivision list oldsharp=#; // start the actual computations string texf; list graph,graphs,texfs; int i,j; for (i=1;i<=size(ff);i++) { def f=ff[i]; if (typeof(f)=="poly") { // exclude the case that the basering has not precisely // one parameter and two indeterminates if ((npars(basering)!=1) or (nvars(basering)!=2)) { ERROR("The basering should have precisely one parameter and two indeterminates!"); } texf=texPolynomial(f); // write the polynomial over Q(t) graph=tropicalCurve(tropicalise(f,#),#); // graph of tropicalis. of f } if (typeof(f)=="list") { if (size(#)==0) { texf="\\min\\{"; } else { texf="\\max\\{"; } for (j=1;j<=size(f);j++) { texf=texf+texPolynomial(f[j]); if (j1) { TDT=TDT+" "+texDrawTropical(graphs[i],#); } else { #=insert(#,"noweights"); TDT=TDT+" "+texDrawTropical(graphs[i],#); #=delete(#,1); } } // add lattice points if the scalefactor is >= 1/2 // and was not handed over def scalefactor=SCF[1]; if (scalefactor>1/2) { def minx=SCF[4]; def miny=SCF[5]; def maxx=SCF[6]; def maxy=SCF[7]; def centerx=SCF[8]; def centery=SCF[9]; int uh=1; if (scalefactor>3) { uh=0; } TDT=TDT+" %% HERE STARTS THE CODE FOR THE LATTICE"; for (i=int(minx)-uh;i<=int(maxx)+uh;i++) { for (j=int(miny)-uh;j<=int(maxy)+uh;j++) { TDT=TDT+" \\move ("+decimal(i-centerx)+" "+decimal(j-centery)+") \\fcir f:0.8 r:"+decimal(1/(10*scalefactor),size(string(int(scalefactor)))+1); } } TDT=TDT+" %% HERE ENDS THE CODE FOR THE LATTICE "; } TEXBILD=TEXBILD+" \\begin{center} "+texDrawBasic(TDT)+ // write the tropical curve "\\end{center} \\vspace*{0.5cm} "; // compute the scaling factor for the Newton subdivisions oldsharp[size(oldsharp)+1]=nsdScaleFactor(graphs); for (i=1;i<=size(ff);i++) { if (i==1) { fname="f_1\\cdots f_{"+string(size(ff)-1)+"}"; } else { fname="f_{"+string(i-1)+"}"; } TEXBILD=TEXBILD+" The Newton subdivision of the tropical curve $"+fname+"$ is: \\vspace*{0.5cm} \\begin{center} "+texDrawNewtonSubdivision(graphs[i],oldsharp)+" \\end{center} \\vspace*{0.5cm} "; } TEXBILD=TEXBILD+" \\end{document}"; if(defined(onlytexfile)==0) { int rdnum=random(1000,9999); write(":w /tmp/tropicalcurve"+string(rdnum)+".tex",TEXBILD); 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 &"); } else { return(TEXBILD); } } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),dp; // poly f=t*(x3+y3+1)+1/t*(x2+y2+x+y+x2y+xy2)+1/t2*xy; poly f=x+y+1; // the command drawTropicalCurve(f) computes the graph of the tropical curve // given by f and displays a post script image, provided you have kghostview // we can instead apply the procedure to a tropical polynomial and use "maximum" poly g=t3*(x7+y7+1)+1/t3*(x4+y4+x2+y2+x3y+xy3)+1/t21*x2y2; // list tropical_g=tropicalise(g); drawTwoTropicalCurves(list(f,g),"max"); } proc minScaleFactor (list graphs) { list graph; int j,i,k; poly minx,miny,maxx,maxy; poly minX,minY,maxX,maxY; poly maxdiffx,maxdiffy; poly centerx,centery; int nachkomma; number sf,scf; poly scalefactor; list SFCS; for (k=1;k<=size(graphs);k++) { graph=graphs[k]; // find the minimal and maximal coordinates of vertices minx,miny,maxx,maxy=graph[1][1],graph[1][2],graph[1][1],graph[1][2]; for (i=2;i<=size(graph)-1;i++) { minx=minOfPolys(list(minx,graph[i][1])); miny=minOfPolys(list(miny,graph[i][2])); maxx=-minOfPolys(list(-maxx,-graph[i][1])); maxy=-minOfPolys(list(-maxy,-graph[i][2])); } if (k==1) { minX=minx; minY=miny; maxX=maxx; maxY=maxy; } else { if (minxmaxX) { maxX=maxx; } if (maxy>maxY) { maxY=maxy; } } } minx,miny,maxx,maxy=minX,minY,maxX,maxY; // find the scale factor for the texdraw image maxdiffx=maxx-minx; maxdiffy=maxy-miny; centerx,centery=int(minx+maxdiffx/2),int(miny+maxdiffy/2); if (maxdiffx==0) { maxdiffx=1; } if (maxdiffy==0) { maxdiffy=1; } nachkomma=2; // number of decimals for the scalefactor sf=1; // correction factor for scalefactor scalefactor=minOfPolys(list(12/leadcoef(maxdiffx),16/leadcoef(maxdiffy))); // if the scalefactor is less than 1/100, then we need more than 2 decimals if (leadcoef(scalefactor) < 1/100) { scf=leadcoef(scalefactor); while (scf < 1/100) { scf=scf * 10; nachkomma++; } } // if the scalefactor is < 1/100, then we should rather scale the // coordinates directly, since otherwise texdraw gets into trouble if (nachkomma > 2) { for (i=3;i<=nachkomma;i++) { scalefactor=scalefactor * 10; sf=sf*10; } } return(list(scalefactor,sf,nachkomma,minx,miny,maxx,maxy,centerx,centery)); } example { "EXAMPLE:"; echo=2; ring r=(0,t),(x,y),dp; poly f=t*(x3+y3+1)+1/t*(x2+y2+x+y+x2y+xy2)+1/t2*xy; poly g=1/t3*(x7+y7+1)+t3*(x4+y4+x2+y2+x3y+xy3)+t21*x2y2; list graphs; graphs[1]=tropicalCurve(f); graphs[2]=tropicalCurve(g); minScaleFactor(graphs); } proc nsdScaleFactor (list graphs) { int i,j; list graph,boundary,scalefactors; poly maxx,maxy; for (j=1;j<=size(graphs);j++) { graph=graphs[j]; boundary=graph[size(graph)][1]; // find maximal and minimal x- and y-coordinates and define the scalefactor maxx,maxy=1,1; for (i=1;i<=size(boundary);i++) { maxx=-minOfPolys(list(-maxx,-boundary[i][1])); maxy=-minOfPolys(list(-maxy,-boundary[i][2])); } scalefactors[j]=minOfPolys(list(12/leadcoef(maxx),12/leadcoef(maxy))); } return(minOfPolys(scalefactors)); }