////////////////////////////////////////////////////////////////////////////// version="version grobcov.lib 4.1.2. January_2020 "; // $Id$; // version N9; January 2020; info=" LIBRARY: grobcov.lib \"Groebner Cover for parametric ideals.\", Comprehensive Groebner Systems, Groebner Cover, Canonical Forms, Parametric Polynomial Systems, Automatic Deduction of Geometric Theorems, Dynamic Geometry, Loci, Envelope, Constructible sets. See: A. Montes A, M. Wibmer, \"Groebner Bases for Polynomial Systems with parameters\", Journal of Symbolic Computation 45 (2010) 1391-1425. (https://www.mat.upc.edu//en/people/antonio.montes/). IMPORTANT: Recently published book: A. Montes. \" The Groebner Cover\": Springer, Algorithms and Computation in Mathematics 27 (2019) ISSN 1431-1550 ISBN 978-3-030-03903-5 ISBN 978-3-030-03904-2 (e-Book) Springer Nature Switzerland AG 2018 https://www.springer.com/gp/book/9783030039035 The book can also be used as a user manual of all the routines included in this library. It defines and proves all the theoretic results used in the library, and shows examples of all the routines. There are many previous papers realted to the subject, and the book actualices all the contents. AUTHORS: Antonio Montes (Universitat Politecnica de Catalunya), Hans Schoenemann (Technische Universitaet Kaiserslautern). OVERVIEW: In 2010, the library was designed to contain Montes-Wibmer's algorithm for computing the Canonical Groebner Cover of a parametric ideal. The central routine is grobcov. Given a parametric ideal, grobcov outputs its Canonical Groebner Cover, consisting of a set of triplets of (lpp, basis, segment). The basis (after normalization) is the reduced Groebner basis for each point of the segment. The segments are disjoint, locally closed and correspond to constant lpp (leading power product) of the basis, and are represented in canonical representation. The segments cover the whole parameter space. The output is canonical, it only depends on the given parametric ideal and the monomial order, because the segments have different lpph of the homogenized system. This is much more than a simple Comprehensive Groebner System. The algorithm grobcov allows options to solve partially the problem when the whole automatic algorithm does not finish in reasonable time. Its existence was proved for the first time by Michael Wibmer \"Groebner bases for families of affine or projective schemes\", JSC, 42,803-834 (2007). grobcov uses a first algorithm cgsdr that outputs a disjoint reduced Comprehensive Groebner System with constant lpp. For this purpose, in this library, the implemented algorithm is Kapur-Sun-Wang algorithm, because it is actually the most efficient algorithm known for this purpose. D. Kapur, Y. Sun, and D.K. Wang \"A New Algorithm for Computing Comprehensive Groebner Systems\". Proceedings of ISSAC'2010, ACM Press, (2010), 29-36. The library has evolved to include new applications of the Groebner Cover, and new theoretical developments have been done. A routine locus has been included to compute loci of points, and determining the taxonomy of the components. Additional routines to transform the output to string (locusdg, locusto) are also included and used in the Dynamic Geometry software GeoGebra. They were described in: M.A. Abanades, F. Botana, A. Montes, T. Recio: \''An Algebraic Taxonomy for Locus Computation in Dynamic Geometry\''. Computer-Aided Design 56 (2014) 22-33. Routines for determining the generalized envelope of a family of hypersurfaces (envelop, AssocTanToEnv, FamElemsToEnvCompPoints) are also included. It also includes procedures for Automatic Deduction of Geometric Theorems (ADGT). The actual version also includes a routine (ConsLevels) for computing the canonical form of a constructible set, given as a union of locally closed sets. It determines the canonical levels of a constructible set. It is described in: J.M. Brunat, A. Montes, \"Computing the canonical representation of constructible sets\". Math. Comput. Sci. (2016) 19: 165-178. A complementary routine Levels transforms the output of ConsLevels into the proper locally closed sets forming the levels of the constructible. Antoher complementary routine Grob1Levels has been included to select the locally closed sets of the segments of the grobcov that correspond to basis different from 1, add them together and return the canonical form of this constructible set. More recently (2019) given two locally closed sets in canonical form the new routine DifConsLCSets determines a set of locally closed sets equivalent to the difference them. The description of the routine is submitted to the Journal of Symbolic Computation. This routine can be also used internally by ADGT with the option \"neg\",1 . With this option DifConsLCSets is used for the negative hypothesis and thesis in ADGT. This version was finished on 1/2/2020, NOTATIONS: Before calling any routine of the library grobcov, the user must define the ideal Q[a][x], and all the input polynomials and ideals defined on it. Internally the routines define and use also other ideals: Q[a], Q[x,a] and so on. PROCEDURES: grobcov(F); Is the basic routine giving the canonical Groebner Cover of the parametric ideal F. This routine accepts many options, that allow to obtain results even when the canonical computation does not finish in reasonable time. cgsdr(F); Is the procedure for obtaining a first disjoint, reduced Comprehensive Groebner System that is used in grobcov, but can also be used independently if only a CGS is required. It is a more efficient routine than buildtree (the own routine of 2010 that is no more available). Now, Kapur-Sun-Wang (KSW) algorithm is used. pdivi(f,F); Performs a pseudodivision of a parametric polynomial by a parametric ideal. pnormalf(f,E,N); Reduces a parametric polynomial f over V(E) \ V(N). E is the null ideal and N the non-null ideal over the parameters. Crep(N,M); Computes the canonical C-representation of V(N) \ V(M). It can be called in Q[a] or in Q[a][x], but the ideals N,M can only contain parameters of Q[a]. Prep(N,M); Computes the canonical P-representation of V(N) \ V(M). It can be called in Q[a] or in Q[a][x], but the ideals N,M can only contain parameters of Q[a]. PtoCrep(L) Starting from the canonical Prep of a locally closed set computes its Crep. extendpoly(f,p,q); Given the generic representation f of an I-regular function F defined by poly f on V(p) \ V(q) it returns its full representation. extendGC(GC); When the grobcov of an ideal has been computed with the default option (\"ext\",0) and the explicit option (\"rep\",2) (which is not the default), then one can call extendGC(GC) (and options) to obtain the full representation of the bases. With the default option (\"ext\",0) only the generic representation of the bases is computed, and one can obtain the full representation using extendGC. locus(G); Special routine for determining the geometrical locus of points verifying given conditions. Given a parametric ideal J in Q[x1,..,xn][u1,..,um] with parameters (x1,..,xn) and variables (u1,..,um), representing the system determining the n-dimensional locus with tracer point (x1,..,xn) verifying certain properties, one can apply locus to the system F, for obtaining the locus. locus provides all the components of the locus and determines their taxonomy, that can be: \"Normal\", \"Special\", \"Accumulation\", \"Degenerate\". The mover are the u-variables. The user can ventually restrict them to a subset of them for geometrical reasons but this can change the true taxonomy. locusdg(G); Is a special routine that determines the \"Relevant\" components of the locus in dynamic geometry. It is to be called to the output of locus and selects from it the \"Normal\", and\"Accumulation\" components. envelop(F,C); Special routine for determining the envelop of a family of hyper-surfaces F in Q[x1,..,xn][t1,..,tm] depending on an ideal of constraints C in Q[t1,..,tm]. It computes the locus of the envelop, and detemines the different components as well as its taxonomy: \"Normal\", \"Special\", \"Accumulation\", \"Degenerate\". (See help for locus). locusto(L); Transforms the output of locus, locusdg, envelop into a string that can be reed from different computational systems. stdlocus(F); Simple procedure to determine the components of the locus, alternative to locus that uses only standard GB computation. Cannot determine the taxonomy of the irreducible components. AssocTanToEnv(F,C,E); Having computed an envelop component E of a family of hyper-surfaces F, with constraints C, it returns the parameter values of the associated tangent hyper-surface of the family passing at one point of the envelop component E. FamElemsAtEnvCompPoints(F,C,E) Having computed an envelop component E of a family of hyper-surfaces F, with constraints C, it returns the parameter values of all the hyper-surfaces of the family passing at one point of the envelop component E. discrim(f,x); Determines the factorized discriminant of a degree 2 polynomial in the variable x. The polynomial can be defined on any ring where x is a variable. The polynomial f can depend on parameters and variables. WLemma(F,A); Given an ideal F in Q[a][x] and an ideal A in Q[a], it returns the list (lpp,B,S) were B is the reduced Groebner basis of the specialized F over the segment S, subset of V(A) with top A, determined by Wibmer's Lemma. S is determined in P-representation (or optionally in C-representation). The basis is given by I-regular functions. WLCGS(F); Given a parametric ideal F in Q[a][x] determines a CGS in full-representation using WLemma intersectpar(L); Auxiliary routine. Given a list of ideals definend on K[a][x] it determines the intersection of all of them in K[x,a] ADGT(H,T,H1,T1); Given 4 ideals H,T,H1,T1 in Q[a][x], corresponding to a problem of Automatic Deduction of Geometric Theorems, it determines the supplementary conditions over the parameters for the Proposition (H and not H1) => (T and not T1) to be a Theorem. If H1=1 then H1 is not considered, and analogously for T1. ConsLevels(A); Given a list of locally colsed sets, constructs the canonical representation of the levels of A an its complement. Levels(L); Transforms the output of ConsLevels into the proper Levels of the constructible set. Grob1Levels(G); From the output of grobcov, Grob1Levels selects the segments of G with basis different from 1 (having solutions), and determines the levels of the constructible set formed by them. DifConsLCSets(A,B); given the canonical forms of the constructible sets A and B, A=[a1,a2,..,ak], B=[b1,b2,...,bj], DifConsLCSets returns a list of locally closed sets of the set A minus B, that can be transformed into the canonical form of A minus B applying ConsLevels. SEE ALSO: compregb_lib "; LIB "poly.lib"; LIB "primdec.lib"; LIB "qhmoduli.lib"; // ************ Begin of the grobcov library ********************* // Development of the library: // Library grobcov.lib // (Groebner Cover): // Release 0: (public) // Initial data: 21-1-2008 // Uses buildtree for cgsdr // Final data: 3-7-2008 // Release 2: (prived) // Initial data: 6-9-2009 // Last version using buildtree for cgsdr // Final data: 25-10-2011 // Release B: (prived) // Initial data: 1-7-2012 // Uses both buildtree and KSW for cgsdr // Final data: 4-9-2012 // Release D. Includes routine locus // Release G: (public) // Initial data: 4-9-2012 // Uses KSW algorithm for cgsdr // Final data: 21-11-2013 // Release K: Includes routine envelop // Release L: (public) // New routine ConsLevels: 25-1-2016 // New routine Levels: 25-1-2016 // New routine Grob1Levels: 25-1-2016 // Updated locus: 10-7-2015 (uses ConsLevels) // Release M: (public) // New routines for computing the envelope of a family of // hyper-surfaces and associated questions: 22-4-2016: 20-9-2016 // New routine WLemma for computing the result of // Wibmer's Lemma: 19-9-2016 // Final data October 2016 // Updated locus (messages) // Updated locus (messages) // Final data Mars 2017 // Release N4: (public) // New routine ADGT for Automatic Discovery of Geometric Theorems: 21-1-2018 // Final data February 2018 // Release N8: July 2019. Actualized versions of the routines and options // Release N9: December 2019. New routine DifConsLCSets, // Updated also ADGT to use as option DifConsLCSets //*************Auxiliary routines************** // elimintfromideal: elimine the constant numbers from an ideal // (designed for W, nonnull conditions) // Input: ideal J // Output:ideal K with the elements of J that are non constants, in the // ring Q[x1,..,xm] static proc elimintfromideal(ideal J) { int i; int j=0; ideal K; if (size(J)==0){return(ideal(0));} for (i=1;i<=ncols(J);i++){if (size(variables(J[i])) !=0){j++; K[j]=J[i];}} return(K); } // elimfromlistel: elimine the ideal J from the list L // Input: list L; list of ideals // ideal J; a possible element of L // Output:ideal K with the elements of L different from J // ring Q[x1,..,xm] static proc elimidealfromlist(list L,ideal J) { int i; int j=0; list K; if (size(L)==0){return(L);} for (i=1;i<=size(L);i++) { if (not(equalideals(J,L[i]))) { j++; K[j]=L[i]; } } return(K); } // delfromideal: deletes the i-th polynomial from the ideal F // Works in any kind of ideal static proc delfromideal(ideal F, int i) { int j; ideal G; if (size(F)1) { for(j=1;j<=i-1;j++) { L[size(L)+1]=l[j]; } } for(j=i+1;j<=size(l);j++) { L[size(L)+1]=l[j]; } return(L); } // eliminates repeated elements form an ideal or matrix or module or intmat or bigintmat static proc elimrepeated(F) { int i; int nt; if (typeof(F)=="ideal"){nt=ncols(F);} else{nt=size(F);} def FF=F; FF=F[1]; for (i=2;i<=nt;i++) { if (not(memberpos(F[i],FF)[1])) { FF[size(FF)+1]=F[i]; } } return(FF); } // equalideals // Input: ideals F and G; // Output: 1 if they are identical (the same polynomials in the same order) // 0 else static proc equalideals(ideal F, ideal G) { int i=1; int t=1; if (size(F)!=size(G)){return(0);} while ((i<=size(F)) and (t==1)) { if (F[i]!=G[i]){t=0;} i++; } return(t); } // returns 1 if the two lists of ideals are equal and 0 if not static proc equallistideals(list L, list M) { int t; int i; if (size(L)!=size(M)){return(0);} else { t=1; if (size(L)>0) { i=1; while ((t) and (i<=size(L))) { if (equalideals(L[i],M[i])==0){t=0;} i++; } } return(t); } } // idcontains // Input: ideal p, ideal q // Output: 1 if p contains q, 0 otherwise // If the routine is to be called from the top, a previous call to static proc idcontains(ideal p, ideal q) { int t; int i; t=1; i=1; def P=p; def Q=q; attrib(P,"isSB",1); poly r; while ((t) and (i<=size(Q))) { r=reduce(Q[i],P,5); if (r!=0){t=0;} i++; } return(t); } // selectminideals // given a list of ideals returns the list of integers corresponding // to the minimal ideals in the list // Input: L (list of ideals) // Output: the list of integers corresponding to the minimal ideals in L // Works in Q[u_1,..,u_m] static proc selectminideals(list L) { list P; int i; int j; int t; if(size(L)==0){return(L)}; if(size(L)==1){P[1]=1; return(P);} for (i=1;i<=size(L);i++) { t=1; j=1; while ((t) and (j<=size(L))) { if (i!=j) { if(idcontains(L[i],L[j])==1) { t=0; } } j++; } if (t){P[size(P)+1]=i;} } return(P); } static proc memberpos(f,J) //"USAGE: memberpos(f,J); // (f,J) expected (polynomial,ideal) // or (int,list(int)) // or (int,intvec) // or (intvec,list(intvec)) // or (list(int),list(list(int))) // or (ideal,list(ideal)) // or (list(intvec), list(list(intvec))). // Works in any kind of ideals //RETURN: The list (t,pos) t int; pos int; // t is 1 if f belongs to J and 0 if not. // pos gives the position in J (or 0 if f does not belong). //EXAMPLE: memberpos; shows an example" { int pos=0; int i=1; int j; int t=0; int nt; if (typeof(J)=="ideal"){nt=ncols(J);} else{nt=size(J);} if ((typeof(f)=="poly") or (typeof(f)=="int")) { // (poly,ideal) or // (poly,list(poly)) // (int,list(int)) or // (int,intvec) i=1; while(i<=nt) { if (f==J[i]){return(list(1,i));} i++; } return(list(0,0)); } else { if ((typeof(f)=="intvec") or ((typeof(f)=="list") and (typeof(f[1])=="int"))) { // (intvec,list(intvec)) or // (list(int),list(list(int))) i=1; t=0; pos=0; while((i<=nt) and (t==0)) { t=1; j=1; if (size(f)!=size(J[i])){t=0;} else { while ((j<=size(f)) and t) { if (f[j]!=J[i][j]){t=0;} j++; } } if (t){pos=i;} i++; } if (t){return(list(1,pos));} else{return(list(0,0));} } else { if (typeof(f)=="ideal") { // (ideal,list(ideal)) i=1; t=0; pos=0; while((i<=nt) and (t==0)) { t=1; j=1; if (ncols(f)!=ncols(J[i])){t=0;} else { while ((j<=ncols(f)) and t) { if (f[j]!=J[i][j]){t=0;} j++; } } if (t){pos=i;} i++; } if (t){return(list(1,pos));} else{return(list(0,0));} } else { if ((typeof(f)=="list") and (typeof(f[1])=="intvec")) { // (list(intvec),list(list(intvec))) i=1; t=0; pos=0; while((i<=nt) and (t==0)) { t=1; j=1; if (size(f)!=size(J[i])){t=0;} else { while ((j<=size(f)) and t) { if (f[j]!=J[i][j]){t=0;} j++; } } if (t){pos=i;} i++; } if (t){return(list(1,pos));} else{return(list(0,0));} } } } } } //example //{ "EXAMPLE:"; echo = 2; // list L=(7,4,5,1,1,4,9); // memberpos(1,L); //} // Auxiliary routine // pos // Input: intvec p of zeros and ones // Output: intvec W of the positions where p has ones. static proc pos(intvec p) { int i; intvec W; int j=1; for (i=1; i<=size(p); i++) { if (p[i]==1){W[j]=i; j++;} } return(W); } // Input: // A,B: lists of ideals // Output: // 1 if both lists of ideals are equal, or 0 if not static proc equallistsofideals(list A, list B) { int i; int tes=0; if (size(A)!=size(B)){return(tes);} tes=1; i=1; while(tes==1 and i<=size(A)) { if (equalideals(A[i],B[i])==0){tes=0; return(tes);} i++; } return(tes); } // Input: // A,B: lists of P-rep, i.e. of the form [p_i,[p_{i1},..,p_{ij_i}]] // Output: // 1 if both lists of P-reps are equal, or 0 if not static proc equallistsA(list A, list B) { int tes=0; if(equalideals(A[1],B[1])==0){return(tes);} tes=equallistsofideals(A[2],B[2]); return(tes); } // Input: // A,B: lists lists of of P-rep, i.e. of the form [[p_1,[p_{11},..,p_{1j_1}]] .. [p_i,[p_{i1},..,p_{ij_i}]] // Output: // 1 if both lists of lists of P-rep are equal, or 0 if not static proc equallistsAall(list A,list B) { int i; int tes; if(size(A)!=size(B)){return(tes);} tes=1; i=1; while(tes and i<=size(A)) { tes=equallistsA(A[i],B[i]); i++; } return(tes); } // idint: ideal intersection // in the ring @P. // it works in an extended ring // input: two ideals in the ring @P // output the intersection of both (is not a GB) static proc idint(ideal I, ideal J) { def RR=basering; ring T=0,t,lp; def K=T+RR; setring(K); def Ia=imap(RR,I); def Ja=imap(RR,J); ideal IJ; int i; for(i=1;i<=size(Ia);i++){IJ[i]=t*Ia[i];} for(i=1;i<=size(Ja);i++){IJ[size(Ia)+i]=(1-t)*Ja[i];} ideal eIJ=eliminate(IJ,t); setring(RR); return(imap(K,eIJ)); } //purpose ideal intersection called in @R and computed in @P static proc idintR(ideal N, ideal M) { def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); setring(P); def Np=imap(RR,N); def Mp=imap(RR,M); def Jp=idint(Np,Mp); setring(RR); return(imap(P,Jp)); } // Auxiliary routine // comb: the list of combinations of elements (1,..n) of order p static proc comb(int n, int p) { list L; list L0; intvec c; intvec d; int i; int j; int last; if ((n<0) or (n=g static proc lesspol(poly f, poly g) { if (leadmonom(f) static proc incquotient(ideal N, poly f) { poly g; int i; ideal Nb; ideal Na=N; if (size(Na)==1) { g=gcd(Na[1],f); if (g!=1) { Na[1]=Na[1]/g; } attrib(Na,"IsSB",1); return(Na); } def P=basering; poly @t; ring H=0,@t,lp; def HP=H+P; setring(HP); def fh=imap(P,f); def Nh=imap(P,N); ideal Nht; for (i=1;i<=size(Nh);i++) { Nht[i]=Nh[i]*@t; } attrib(Nht,"isSB",1); def fht=(1-@t)*fh; option(redSB); Nht=std(Nht,fht); ideal Nc; ideal v; for (i=1;i<=size(Nht);i++) { v=variables(Nht[i]); if(memberpos(@t,v)[1]==0) { Nc[size(Nc)+1]=Nht[i]/fh; } } setring(P); ideal HH; def Nd=imap(HP,Nc); Nb=Nd; option(redSB); Nb=std(Nd); return(Nb); } // Auxiliary routine to define an order for ideals // Returns 1 if the ideal a is shoud precede ideal b by sorting them in idbefid order // 2 if the the contrary happen // 0 if the order is not relevant static proc idbefid(ideal a, ideal b) { poly fa; poly fb; poly la; poly lb; int te=1; int i; int j; int na=size(a); int nb=size(b); int nm; if (na<=nb){nm=na;} else{nm=nb;} for (i=1;i<=nm; i++) { fa=a[i]; fb=b[i]; while((fa!=0) or (fb!=0)) { la=lead(fa); lb=lead(fb); fa=fa-la; fb=fb-lb; la=leadmonom(la); lb=leadmonom(lb); if(leadmonom(la+lb)!=la){return(1);} else{if(leadmonom(la+lb)!=lb){return(2);}} } } if(nanb){return(2);} else{return(0);} } } // sort a list of ideals using idbefid static proc sortlistideals(list L) { int i; int j; int n; ideal a; ideal b; list LL=L; list NL; int k; int te; i=1; while(size(LL)>0) { k=1; for(j=2;j<=size(LL);j++) { te=idbefid(LL[k],LL[j]); if (te==2){k=j;} } NL[size(NL)+1]=LL[k]; n=size(LL); if (n>1){LL=elimfromlist(LL,k);} else{LL=list();} } return(NL); } // Crep // Computes the C-representation of V(N) \ V(M). // input: // ideal N (null ideal) (not necessarily radical nor maximal) // ideal M (hole ideal) (not necessarily containing N) // output: // the list (a,b) of the canonical ideals // the Crep of V(N) \ V(M) // Assumed to be called in the ring Q[a] or Q[x] proc Crep(ideal N, ideal M) "USAGE: Crep(ideal N,ideal M); ideal N (null ideal) (not necessarily radical nor maximal) in Q[a]. (a=parameters, x=variables). ideal M (hole ideal) (not necessarily containing N) in Q[a]. To be called in a ring Q[a][x] or a ring Q[a]. But the ideals can contain only the parameters in Q[a]. RETURN: The canonical C-representation [P,Q] of the locally closed set, formed by a pair of radical ideals with P included in Q, representing the set V(P) \ V(Q) = V(N) \ V(M) KEYWORDS: locally closed set; canoncial form EXAMPLE: Crep; shows an example" { int te; def RR=basering; def Rx=ringlist(RR); if(size(Rx[1])==4) { te=1; def P=ring(Rx[1]); } if(te==1) { setring(P); ideal Np=imap(RR,N); ideal Mp=imap(RR,M); } else {te=0; def Np=N; def Mp=M;} def La=Crep0(Np,Mp); if(size(La)==0) { if(te==1) {setring(RR); list LL;} if(te==0){list LL;} return(LL); } else { if(te==1) {setring(RR); def LL=imap(P,La);} if(te==0){def LL=La;} return(LL); } } example { "EXAMPLE:"; echo = 2; short=0; if(defined(R)){kill R;} ring R=0,(a,b,c),lp; ideal p=a*b; ideal q=a,b-2; // C-representation of V(p) \ V(q) Crep(p,q); } // Crep0 // Computes the C-representation of V(N) \ V(M). // input: // ideal N (null ideal) (not necessarily radical nor maximal) // ideal M (hole ideal) (not necessarily containing N) // output: // the list (a,b) of the canonical ideals // the Crep0 of V(N) \ V(M) // Assumed to be called in a ring Q[x] (example @P) static proc Crep0(ideal N, ideal M) { list l; ideal Np=std(N); if (equalideals(Np,ideal(1))) { l=ideal(1),ideal(1); return(l); } int i; list L; ideal Q=Np+M; ideal P=ideal(1); L=minGTZ(Np); for(i=1;i<=size(L);i++) { L[i]=std(L[i]); if(idcontains(L[i],Q)==0) { P=intersect(P,L[i]); } } P=std(P); Q=std(radical(Q+P)); if(equalideals(P,Q)){return(l);} list T=P,Q; return(T); } // Prep // Computes the P-representation of V(N) \ V(M). // input: // ideal N (null ideal) (not necessarily radical nor maximal) // ideal M (hole ideal) (not necessarily containing N) // output: // the ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r))); // the Prep of V(N) \ V(M) // Assumed to be called in the ring ring Q[a][x]. But the data must only contain parameters. proc Prep(ideal N, ideal M) "USAGE: Prep(ideal N,ideal M); ideal N (null ideal) (not necessarily radical nor maximal) in Q[a]. (a=parameters, x=variables). ideal M (hole ideal) (not necessarily containing N) in Q[a]. To be called in a ring Q[a][x] or a ring Q[a]. But the ideals can contain only the parameters in Q[a]. RETURN: The canonical P-representation of the locally closed set V(N) \ V(M) Output: [Comp_1, .. , Comp_s ] where Comp_i=[p_i,[p_i1,..,p_is_i]] KEYWORDS: locally closed set; canoncial form EXAMPLE: Prep; shows an example" { int te; def RR=basering; def Rx=ringlist(RR); if(size(Rx[1])==4) { def P=ring(Rx[1]); te=1; setring(P); ideal Np=imap(RR,N); ideal Mp=imap(RR,M); } else {te=0; def Np=N; def Mp=M;} def La=Prep0(Np,Mp); if(te==1) {setring(RR); def LL=imap(P,La); } if(te==0){def LL=La;} return(LL); } example { "EXAMPLE:"; echo = 2; short=0; if(defined(R)){kill R;} ring R=0,(a,b,c),lp; ideal p=a*b;; ideal q=a,b-1; // P-representation of V(p) \ V(q) Prep(p,q); } // Prep0 // Computes the P-representation of V(N) \ V(M). // input: // ideal N (null ideal) (not necessarily radical nor maximal) // ideal M (hole ideal) (not necessarily containing N) // output: // the ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r))); // the Prep of V(N) \ V(M) // Assumed to be called in a ring Q[x] (example @P) static proc Prep0(ideal N, ideal M) { int te; if (N[1]==1) { return(list(list(ideal(1),list(ideal(1))))); } int i; int j; list L0; list Ni=minGTZ(N); list prep; for(j=1;j<=size(Ni);j++) { option(redSB); Ni[j]=std(Ni[j]); } list Mij; for (i=1;i<=size(Ni);i++) { Mij=minGTZ(Ni[i]+M); for(j=1;j<=size(Mij);j++) { option(redSB); Mij[j]=std(Mij[j]); } if ((size(Mij)==1) and (equalideals(Ni[i],Mij[1])==1)){;} else { prep[size(prep)+1]=list(Ni[i],Mij); } } //"T_before="; prep; if (size(prep)==0){prep=list(list(ideal(1),list(ideal(1))));} //"T_Prep="; prep; //def Lout=CompleteA(prep,prep); //"T_Lout="; Lout; return(prep); } // PtoCrep // Computes the C-representation from the P-representation. // input: // list ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r))); // the P-representation of V(N) \ V(M) // output: // list (ideal ida, ideal idb) // the C-representaion of V(N) \ V(M) = V(ida) \ V(idb) // Assumed to be called in the ring Q[a] or Q[x] proc PtoCrep(list L) "USAGE: PtoCrep(list L) list L= [ Comp_1, .. , Comp_s ] where list Comp_i=[p_i,[p_i1,..,p_is_i] ], is the P-representation of a locally closed set V(N) \ V(M). To be called in a ring Q[a][x] or a ring Q[a]. But the ideals can contain only the parameters in Q[a]. RETURN:The canonical C-representation [P,Q] of the locally closed set. A pair of radical ideals with P included in Q, representing the set V(P) \ V(Q) KEYWORDS: locally closed set; canoncial form EXAMPLE: PtoCrep; shows an example" { int te; def RR=basering; def Rx=ringlist(RR); if(size(Rx[1])==0){return(PtoCrep0(L));} else { def P=ring(Rx[1]); setring(P); list Lp=imap(RR,L); def LLp=PtoCrep0(Lp); setring(RR); def LL=imap(P,LLp); return(LL); } } example { "EXAMPLE:"; echo = 2; if(defined(R)){kill R;} ring R=0,(a,b,c),lp; short=0; ideal p=a*(a^2+b^2+c^2-25); ideal q=a*(a-3),b-4; // C-representaion of V(p) \ V(q) def Cr=Crep(p,q); Cr; // C-representation of V(p) \ V(q) def L=Prep(p,q); L; PtoCrep(L); } // PtoCrep0 // Computes the C-representation from the P-representation. // input: // list ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r))); // the P-representation of V(N) \ V(M) // output: // list (ideal ida, ideal idb) // the C-representation of V(N) \ V(M) = V(ida) \ V(idb) // Assumed to be called in a ring Q[x] (example @P) static proc PtoCrep0(list L) { int te=0; def Lp=L; int i; int j; ideal ida=ideal(1); ideal idb=ideal(1); list Lb; ideal N; for (i=1;i<=size(Lp);i++) { option(returnSB); //"T_Lp[i]="; Lp[i]; N=Lp[i][1]; Lb=Lp[i][2]; //"T_Lb="; Lb; ida=intersect(ida,N); for(j=1;j<=size(Lb);j++) { idb=intersect(idb,Lb[j]); } } //idb=radical(idb); def La=list(ida,idb); return(La); } // input: F a parametric ideal in Q[a][x] // output: a disjoint and reduced Groebner System. // It uses Kapur-Sun-Wang algorithm, and with the options // can compute the homogenization before (('can',0) or ( 'can',1)) // and dehomogenize the result. proc cgsdr(ideal F, list #) "USAGE: cgsdr(ideal F); F: ideal in Q[a][x] (a=parameters, x=variables) to be discussed. Computes a disjoint, reduced Comprehensive Groebner System (CGS). cgsdr is the starting point of the fundamental routine grobcov. The basering R, must be of the form Q[a][x], (a=parameters, x=variables), and should be defined previously. RETURN: Returns a list T describing a reduced and disjoint Comprehensive Groebner System (CGS). The output is a list of (full,hole,basis), where the ideals full and hole represent the segment V(full) \ V(hole). With option (\"out\",0) the segments are grouped by leading power products (lpp) of the reduced Groebner basis and given in P-representation. The returned list is of the form: [ [lpp, [num,basis,segment],..., [num,basis,segment],lpph], ... , [lpp, [num,basis,segment],..., [num,basis,segment],lpph] ]. The bases are the reduced Groebner bases (after normalization) for each point of the corresponding segment. The third element lpph of each lpp segment is the lpp of the homogenized ideal used ideal in the CGS as a string, that is shown only when option (\"can\",1) is used. With option (\"can\",0) the homogenized basis is used. With option (\"can\",1) the homogenized ideal is used. With option (\"can\",2) the given basis is used. With option (\"out\",1) (default) only KSW is applied and segments are given as difference of varieties and are not grouped The returned list is of the form: [[E,N,B],..[E,N,B]] E is the top variety N is the hole variety. Segment = V(E) \ V(N) B is the reduced Groebner basis OPTIONS: An option is a pair of arguments: string, integer. To modify the default options, pairs of arguments -option name, value- of valid options must be added to the call. Inside grobcov the default option is \"can\",1. It can be used also with option \"can\",0 but then the output is not the canonical Groebner Cover. grobcov cannot be used with option \"can\",2. When cgsdr is called directly, the options are \"can\",0-1-2: The default value is \"can\",2. In this case no homogenization is done. With option (\"can\",0) the given basis is homogenized, and with option (\"can\",1) the whole given ideal is homogenized before computing the cgs and dehomogenized after. With option (\"can\",0) the homogenized basis is used. With option (\"can\",1) the homogenized ideal is used. With option (\"can\",2) the given basis is used. \"null\",ideal E: The default is (\"null\",ideal(0)). \"nonnull\",ideal N: The default (\"nonnull\",ideal(1)). When options \"null\" and/or \"nonnull\" are given, then the parameter space is restricted to V(E) \ V(N). \"comment\",0-1: The default is (\"comment\",0). Setting (\"comment\",1) will provide information about the development of the computation. \"out\",0-1: (default is 1) the output segments are given as as difference of varieties. With option \"out\",0 the output segments are given in P-representation and the segments grouped by lpp. With options (\"can\",0) and (\"can\",1) the option (\"out\",1) is set to (\"out\",0) because it is not compatible. One can give none or whatever of these options. With the default options (\"can\",2,\"out\",1), only the Kapur-Sun-Wang algorithm is computed. The algorithm used is: D. Kapur, Y. Sun, and D.K. Wang \"A New Algorithm for Computing Comprehensive Groebner Systems\". Proceedings of ISSAC'2010, ACM Press, (2010), 29-36. It is very efficient but is only the starting point for the computation of grobcov. When grobcov is computed, the call to cgsdr inside uses specific options that are more expensive (\"can\",0-1,\"out\",0). KEYWORDS: CGS; disjoint; reduced; Comprehensive Groebner System EXAMPLE: cgsdr; shows an example" { int te; def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); // if(size(Rx[1])==4){te=1; Rx[1]=0; def D=ring(Rx); def RP=D+P;} // else{te=0;} //setglobalrings();} // INITIALIZING OPTIONS int i; int j; def E=ideal(0); def N=ideal(1); int comment=0; int can=2; int out=1; poly f; ideal B; int start=timer; list L=#; for(i=1;i<=size(L) div 2;i++) { if(L[2*i-1]=="null"){E=L[2*i];} else { if(L[2*i-1]=="nonnull"){N=L[2*i];} else { if(L[2*i-1]=="comment"){comment=L[2*i];} else { if(L[2*i-1]=="can"){can=L[2*i];} else { if(L[2*i-1]=="out"){out=L[2*i];} } } } } } //if(can==2){out=1;} B=F; if ((printlevel) and (comment==0)){comment=printlevel;} if((can<2) and (out>0)){"Option out,1 is not compatible with can,0,1"; out=0;} // DEFINING OPTIONS list LL; LL[1]="can"; LL[2]=can; LL[3]="comment"; LL[4]=comment; LL[5]="out"; LL[6]=out; LL[7]="null"; LL[8]=E; LL[9]="nonnull"; LL[10]=N; if(comment>=1) { " "; string("Begin cgsdr with options: ",LL); } int ish; for (i=1;i<=size(B);i++){ish=ishomog(B[i]); if(ish==0){break;};} if (ish) { if(comment>0){" "; string("The given system is homogneous");} def GS=KSW(B,LL); //can=0; } else { // ACTING DEPENDING ON OPTIONS if(can==2) { // WITHOUT HOMOHGENIZING if(comment>0){" "; string("Option of cgsdr: do not homogenize");} def GS=KSW(B,LL); // setglobalrings(); } else { if(can==1) { // COMPUTING THE HOMOGOENIZED IDEAL if(comment>0){" "; string("Homogenizing the whole ideal: option can=1");} list RRL=ringlist(RR); RRL[3][1][1]="dp"; def Pa=ring(RRL[1]); list Lx; Lx[1]=0; Lx[2]=RRL[2]+RRL[1][2]; Lx[3]=RRL[1][3]; Lx[4]=RRL[1][4]; RRL[1]=0; def D=ring(RRL); def RP=D+Pa; setring(RP); def B1=imap(RR,B); option(redSB); if(comment>0){" ";string("Basis before computing its std basis="); B1;} B1=std(B1); if(comment>0){" ";string("Basis after computing its std basis="); B1;} setring(RR); def B2=imap(RP,B1); } else { // (can=0) if(comment>0){" "; string( "Homogenizing the basis: option can=0");} def B2=B; } // COMPUTING HOMOGENIZED CGS poly @t; ring H=0,@t,dp; def RH=RR+H; setring(RH); // setglobalrings(); def BH=imap(RR,B2); def LH=imap(RR,LL); //"T_BH="; BH; //"T_LH="; LH; for (i=1;i<=size(BH);i++) { BH[i]=homog(BH[i],@t); } if (comment>0){" "; string("Homogenized system = "); BH;} def RHx=ringlist(RH); def PH=ring(RHx[1]); RHx[1]=0; def DH=ring(RHx); def RPH=DH+PH; def GSH=KSW(BH,LH); //"T_GSH="; GSH; //setglobalrings(); // DEHOMOGENIZING THE RESULT if(out==0) { for (i=1;i<=size(GSH);i++) { GSH[i][1]=subst(GSH[i][1],@t,1); for(j=1;j<=size(GSH[i][2]);j++) { GSH[i][2][j][2]=subst(GSH[i][2][j][2],@t,1); } } } else { for (i=1;i<=size(GSH);i++) { GSH[i][3]=subst(GSH[i][3],@t,1); GSH[i][7]=subst(GSH[i][7],@t,1); } } setring(RR); def GS=imap(RH,GSH); } // setglobalrings(); if(out==0) { for (i=1;i<=size(GS);i++) { GS[i][1]=postredgb(mingb(GS[i][1])); for(j=1;j<=size(GS[i][2]);j++) { GS[i][2][j][2]=postredgb(mingb(GS[i][2][j][2])); } } } else { for (i=1;i<=size(GS);i++) { if(GS[i][2]==1) { GS[i][3]=postredgb(mingb(GS[i][3])); if (typeof(GS[i][7])=="ideal") { GS[i][7]=postredgb(mingb(GS[i][7]));} } } } } return(GS); } example { "EXAMPLE:"; echo = 2; // Casas conjecture for degree 4: if(defined(R)){kill R;} ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp; short=0; ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0), x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1), x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0), x2^2+(2*a3)*x2+(a2), x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0), x3+(a3); cgsdr(F); } // input: internal routine called by cgsdr at the end to group the // lpp segments and improve the output // output: grouped segments by lpp obtained in cgsdr static proc grsegments(list T) { int i; list L; list lpp; list lp; list ls; int n=size(T); lpp[1]=T[n][1]; L[1]=list(lpp[1],list(list(T[n][2],T[n][3],T[n][4]))); if (n>1) { for (i=1;i<=size(T)-1;i++) { lp=memberpos(T[n-i][1],lpp); if(lp[1]==1) { ls=L[lp[2]][2]; ls[size(ls)+1]=list(T[n-i][2],T[n-i][3],T[n-i][4]); L[lp[2]][2]=ls; } else { lpp[size(lpp)+1]=T[n-i][1]; L[size(L)+1]=list(T[n-i][1],list(list(T[n-i][2],T[n-i][3],T[n-i][4]))); } } } return(L); } // LCUnion // Given a list of the P-representations of locally closed segments // for which we know that the union is also locally closed // it returns the P-representation of its union // input: L list of segments in P-representation // ((p_j^i,(p_j1^i,...,p_jk_j^i | j=1..t_i)) | i=1..s ) // where i represents a segment // output: P-representation of the union // ((P_j,(P_j1,...,P_jk_j | j=1..t))) static proc LCUnion(list LL) { def RR=basering; def Rx=ringlist(RR); def PP=ring(Rx[1]); setring(PP); def L=imap(RR,LL); int i; int j; int k; list H; list C; list T; list L0; list P0; list P; list Q0; list Q; for (i=1;i<=size(L);i++) { for (j=1;j<=size(L[i]);j++) { P0[size(P0)+1]=L[i][j][1]; L0[size(L0)+1]=intvec(i,j); } } Q0=selectminideals(P0); for (i=1;i<=size(Q0);i++) { Q[i]=L0[Q0[i]]; P[i]=L[Q[i][1]][Q[i][2]]; } //"T_P="; P; // P is the list of the maximal components of the union // with the corresponding initial holes. // Q is the list of intvec positions in L of the first element of the P's // Its elements give (num of segment, num of max component (=min ideal)) for (k=1;k<=size(Q);k++) { H=P[k][2]; // holes of P[k][1] for (i=1;i<=size(L);i++) { if (i!=Q[k][1]) { for (j=1;j<=size(L[i]);j++) { C[size(C)+1]=L[i][j]; } } } T[size(T)+1]=list(Q[k],P[k][1],addpart(H,C)); } setring(RR); def TT=imap(PP,T); return(TT); } // LCUnionN // Given a list of the P-representations of locally closed segments // for which we know that the union is also locally closed // it returns the P-representation of its union // input: L list of segments in P-representation // ((p_j^i,(p_j1^i,...,p_jk_j^i | j=1..t_i)) | i=1..s ) // where i represents a segment // output: P-representation of the union // ((P_j,(P_j1,...,P_jk_j | j=1..t))) static proc LCUnionN(list L) { int i; int j; int k; list H; list C; list T; list L0; list P0; list P; list Q0; list Q; //"T_L="; L; for (i=1;i<=size(L);i++) { P0[size(P0)+1]=L[i][1]; for (j=1;j<=size(L[i]);j++) { L0[size(L0)+1]=intvec(i,j); } } //"T_P0="; P0; Q0=selectminideals(P0); //"T_Q0="; Q0; for (i=1;i<=size(Q0);i++) { //Q[i]=L0[Q0[i]]; P[i]=L[Q0[i][1]];// [Q[i][2]]; } //"T_P="; P; // P is the list of the maximal components of the union // with the corresponding initial holes. // Q is the list of intvec positions in L of the first element of the P's // Its elements give (num of segment, num of max component (=min ideal)) // list C; for (k=1;k<=size(Q0);k++) { kill C; list C; H=P[k][2]; // holes of P[k][1] for (i=1;i<=size(L);i++) { if (i!=Q0[k]) // (i!=Q0[k]) { //for (j=1;j<=size(L[i]);j++) //{ C[size(C)+1]=L[i]; //} } } T[size(T)+1]=list(P[k][1],addpart(H,C)); // Q0[k], } return(T); } // Auxiliary routine // called by LCUnion to modify the holes of a primepart of the union // by the addition of the segments that do not correspond to that part // Works on Q[a] ring. // Input: // H=(p_i1,..,p_is) the holes of a component to be transformed by the addition of // the segments C that do not correspond to that component // C=((q_1,(q_11,..,q_1l_1),pos1),..,(q_k,(q_k1,..,q_kl_k),posk)) // posi=(i,j) position of the component // the list of segments to be added to the holes static proc addpart(list H, list C) { list Q; int i; int j; int k; int l; int t; int t1; Q=H; intvec notQ; list QQ; list addq; // plus those of the components added to the holes. ideal q; i=1; while (i<=size(Q)) { if (memberpos(i,notQ)[1]==0) { q=Q[i]; t=1; j=1; while ((t) and (j<=size(C))) { if (equalideals(q,C[j][1])) { t=0; for (k=1;k<=size(C[j][2]);k++) { t1=1; l=1; while((t1) and (l<=size(Q))) { if ((l!=i) and (memberpos(l,notQ)[1]==0)) { if (idcontains(C[j][2][k],Q[l])) { t1=0; } } l++; } if (t1) { addq[size(addq)+1]=C[j][2][k]; } } if((size(notQ)==1) and (notQ[1]==0)){notQ[1]=i;} else {notQ[size(notQ)+1]=i;} } j++; } if (size(addq)>0) { for (k=1;k<=size(addq);k++) { Q[size(Q)+1]=addq[k]; } kill addq; list addq; } } i++; } for (i=1;i<=size(Q);i++) { if(memberpos(i,notQ)[1]==0) { QQ[size(QQ)+1]=Q[i]; } } if (size(QQ)==0){QQ[1]=ideal(1);} return(addpartfine(QQ,C)); } // Auxiliary routine called by addpart to finish the modification of the holes of a primepart // of the union by the addition of the segments that do not correspond to // that part. // Works on Q[a] ring. static proc addpartfine(list H, list C0) { //"T_H="; H; int i; int j; int k; int te; intvec notQ; int l; list sel; intvec jtesC; if ((size(H)==1) and (equalideals(H[1],ideal(1)))){return(H);} if (size(C0)==0){return(H);} list newQ; list nQ; list Q; list nQ1; list Q0; def Q1=H; //Q1=sortlistideals(Q1,idbefid); def C=C0; while(equallistideals(Q0,Q1)==0) { Q0=Q1; i=0; Q=Q1; kill notQ; intvec notQ; while(i1){tes=1;} else { tes=0; for(k=1;k<=size(B);k++){B[k]=pnormalf(B[k],crep[1],crep[2]);} } if(tes) { // combine is needed kill B; ideal B; for (j=1;j<=size(LCU);j++) { LL[j]=LCU[j][2]; } FF=precombine(LL); for (k=1;k<=size(lpp);k++) { kill L; list L; for (j=1;j<=size(LCU);j++) { L[j]=list(LCU[j][2],LCU[j][1][k]); } B[k]=combine(L,FF); } } for(j=1;j<=size(B);j++) { B[j]=pnormalf(B[j],crep[1],crep[2]); } S[i]=list(lpp,B,prep,crep,lpph); if(comment>=1) { lpi[size(lpi)+1]=string("[",i,"]"); lpi[size(lpi)+1]=S[i][1]; } } if(comment>=1) { string("Time in LCUnion + combine = ",timer-start); if(comment>=2){string("lpp=",lpi)}; } return(S); } // grobcov // input: // ideal F: a parametric ideal in Q[a][x], (a=parameters, x=variables). // list #: (options) list("null",N,"nonnull",W,"can",0-1,ext",0-1, "rep",0-1-2) // where // N is the null conditions ideal (if desired) // W is the ideal of non-null conditions (if desired) // The value of \"can\" is 1 by default and can be set to 0 if we do not // need to obtain the canonical GC, but only a GC. // The value of \"ext\" is 0 by default and so the generic representation // of the bases is given. It can be set to 1, and then the full // representation of the bases is given. // The value of \"rep\" is 0 by default, and then the segments // are given in canonical P-representation. It can be set to 1 // and then they are given in canonical C-representation. // If it is set to 2, then both representations are given. // output: // list S: ((lpp,basis,(idp_1,(idp_11,..,idp_1s_1))), .. // (lpp,basis,(idp_r,(idp_r1,..,idp_rs_r))) ) where // each element of S corresponds to a lpp-segment // given by the lpp, the basis, and the P-representation of the segment proc grobcov(ideal F,list #) "USAGE: grobcov(ideal F[,options]); F: ideal in Q[a][x] (a=parameters, x=variables) to be discussed.This is the fundamental routine of the library. It computes the Groebner Cover of a parametric ideal F in Q[a][x]. See A. Montes , M. Wibmer, \"Groebner Bases for Polynomial Systems with parameters\". JSC 45 (2010) 1391-1425.) or the not yet published book A. Montes. \" The Groebner Cover\" (Discussing Parametric Polynomial Systems). The Groebner Cover of a parametric ideal F consist of a set of pairs(S_i,B_i), where the S_i are disjoint locally closed segments of the parameter space, and the B_i are the reducedGroebner bases of the ideal on every point of S_i. The ideal F must be defined on a parametric ring Q[a][x] (a=parameters, x=variables). RETURN: The list [[lpp_1,basis_1,segment_1], ..., [lpp_s,basis_s,segment_s]] optionally [[ lpp_1,basis_1,segment_1,lpph_1], ..., [lpp_s,basis_s,segment_s,lpph_s]] The lpp are constant over a segment and correspond to the set of lpp of the reduced Groebner basis for each point of the segment. With option (\"showhom\",1) the lpph will be shown: The lpph corresponds to the lpp of the homogenized ideal and is different for each segment. It is given as a string, and shown only for information. With the default option \"can\",1, the segments have different lpph. Basis: to each element of lpp corresponds an I-regular function given in full representation (by option (\"ext\",1)) or in generic representation (default option (\"ext\",0)). The I-regular function is the corresponding element of the reduced Groebner basis for each point of the segment with the given lpp. For each point in the segment, the polynomial or the set of polynomials representing it, if they do not specialize to 0, then after normalization, specializes to the corresponding element of the reduced Groebner basis. In the full representation at least one of the polynomials representing the I-regular function specializes to non-zero. With the default option (\"rep\",0) the representation of the segment is the P-representation. With option (\"rep\",1) the representation of the segment is the C-representation. With option (\"rep\",2) both representations of the segment are given. The P-representation of a segment is of the form [[p_1,[p_11,..,p_1k1]],..,[p_r,[p_r1,..,p_rkr]]] representing the segment Union_i ( V(p_i) \ ( Union_j V(p_ij) ) ), where the p's are prime ideals. The C-representation of a segment is of the form (E,N) representing V(E) \ V(N), and the ideals E and N are radical and N contains E. OPTIONS: An option is a pair of arguments: string, integer. To modify the default options, pairs of arguments -option name, value- of valid options must be added to the call. \"null\",ideal E: The default is (\"null\",ideal(0)). \"nonnull\",ideal N: The default is (\"nonnull\",ideal(1)). When options \"null\" and/or \"nonnull\" are given, then the parameter space is restricted to V(E) \ V(N). \"can\",0-1: The default is (\"can\",1). With the default option the homogenized ideal is computed before obtaining the Groebner Cover, so that the result is the canonical Groebner Cover. Setting (\"can\",0) only homogenizes the basis so the result is not exactly canonical, but the computation is shorter. \"ext\",0-1: The default is (\"ext\",0). With the default (\"ext\",0), only the generic representation of the bases is computed (single polynomials, but not specializing to non-zero for every point of the segment. With option (\"ext\",1) the full representation of the bases is computed (possible sheaves) and sometimes a simpler result is obtained, but the computation is more time consuming. \"rep\",0-1-2: The default is (\"rep\",0) and then the segments are given in canonical P-representation. Option (\"rep\",1) represents the segments in canonical C-representation, and option (\"rep\",2) gives both representations. \"comment\",0-3: The default is (\"comment\",0). Setting \"comment\" higher will provide information about the development of the computation. \"showhom\",0-1: The default is (\"showhom\",0). Setting \"showhom\",1 will output the set of lpp of the homogenized ideal of each segment as last element. One can give none or whatever of these options. NOTE: The basering R, must be of the form Q[a][x], (a=parameters, x=variables), and should be defined previously. The ideal must be defined on R. KEYWORDS: Groebner cover; parametric ideal; canonical; discussion of parametric ideal EXAMPLE: grobcov; shows an example" { list S; int i; int ish=1; list GBR; list BR; int j; int k; ideal idp; ideal idq; int s; ideal ext; list SS; ideal E; ideal N; int canop; int extop; int repop; int comment=0; int m; def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); Rx[1]=0; def D=ring(Rx); def RP=D+P; list L0=#; list Se; int out=0; int showhom=0; int hom; L0[size(L0)+1]="res"; L0[size(L0)+1]=ideal(1); // default options int start=timer; E=ideal(0); N=ideal(1); canop=1; // canop=0 for homogenizing the basis but not the ideal (not canonical) // canop=1 for working with the homogenized ideal repop=0; // repop=0 for representing the segments in Prep // repop=1 for representing the segments in Crep // repop=2 for representing the segments in Prep and Crep extop=0; // extop=0 if only generic representation of the bases are to be computed // extop=1 if the full representation of the bases are to be computed for(i=1;i<=size(L0) div 2;i++) { if(L0[2*i-1]=="can"){canop=L0[2*i];} else { if(L0[2*i-1]=="ext"){extop=L0[2*i];} else { if(L0[2*i-1]=="rep"){repop=L0[2*i];} else { if(L0[2*i-1]=="null"){E=L0[2*i];} else { if(L0[2*i-1]=="nonnull"){N=L0[2*i];} else { if (L0[2*i-1]=="comment"){comment=L0[2*i];} else { if (L0[2*i-1]=="showhom"){showhom=L0[2*i];} } } } } } } } if(not((canop==0) or (canop==1))) { string("Option can = ",canop," is not supported. It is changed to can = 1"); canop=1; } for(i=1;i<=size(L0) div 2;i++) { if(L0[2*i-1]=="can"){L0[2*i]=canop;} } if ((printlevel) and (comment==0)){comment=printlevel;} list LL; LL[1]="can"; LL[2]=canop; LL[3]="comment"; LL[4]=comment; LL[5]="out"; LL[6]=0; LL[7]="null"; LL[8]=E; LL[9]="nonnull"; LL[10]=N; LL[11]="ext"; LL[12]=extop; LL[13]="rep"; LL[14]=repop; LL[15]="showhom"; LL[16]=showhom; if (comment>=1) { string("Begin grobcov with options: ",LL); } kill S; def S=gcover(F,LL); // NOW extendGC if(extop) { S=extendGC(S,LL); } // NOW repop and showhom list Si; list nS; for(i=1;i<=size(S);i++) { if(repop==0){Si=list(S[i][1],S[i][2],S[i][3]);} if(repop==1){Si=list(S[i][1],S[i][2],S[i][4]);} if(repop==2){Si=list(S[i][1],S[i][2],S[i][3],S[i][4]);} if(showhom==1){Si[size(Si)+1]=S[i][5];} nS[size(nS)+1]=Si; } S=nS; if (comment>=1) { string("Time in grobcov = ", timer-start); string("Number of segments of grobcov = ", size(S)); } return(S); } example { "EXAMPLE:"; echo = 2; // Casas conjecture for degree 4: if(defined(R)){kill R;} ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp; short=0; ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0), x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1), x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0), x2^2+(2*a3)*x2+(a2), x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0), x3+(a3); grobcov(F); // EXAMPLE: // M. Rychlik robot; // Complexity and Applications of Parametric Algorithms of // Computational Algebraic Geometry.; // In: Dynamics of Algorithms, R. de la Llave, L. Petzold and J. Lorenz eds.; // IMA Volumes in Mathematics and its Applications, // Springer-Verlag 118: 1-29 (2000).; // (18. Mathematical robotics: Problem 4, two-arm robot)." if (defined(R)){kill R;} ring R=(0,a,b,l2,l3),(c3,s3,c1,s1), dp; short=0; ideal S12=a-l3*c3-l2*c1,b-l3*s3-l2*s1,c1^2+s1^2-1,c3^2+s3^2-1; S12; grobcov(S12); } // Auxiliary routine called by extendGC // extendpoly // input: // poly f: a generic polynomial in the basis // ideal idp: such that ideal(S)=idp // ideal idq: such that S=V(idp) \ V(idq) //// NW the list of ((N1,W1),..,(Ns,Ws)) of red-rep of the grouped //// segments in the lpp-segment NO MORE USED // output: proc extendpoly(poly f, ideal idp, ideal idq) "USAGE: extendGC(poly f,ideal p,ideal q); f is a polynomial in Q[a][x] in generic representation of an I-regular function F defined on the locally closed segment S=V(p) \ V(q). p,q are ideals in Q[a], representing the Crep of segment S. RETURN: the extended representation of F in S. It can consist of a single polynomial or a set of polynomials when needed. NOTE: The basering R, must be of the form Q[a][x], (a=parameters,x=variables), and should be defined previously. The ideals must be defined on R. KEYWORDS: Groebner cover; parametric ideal; locally closed set; parametric ideal; generic representation; full representation; I-regular function EXAMPLE: extendpoly; shows an example" { int te=0; def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); Rx[1]=0; def D=ring(Rx); def RP=D+P; matrix CC; poly Q; list NewMonoms; int i; int j; poly fout; ideal idout; list L=monoms(f); int nummonoms=size(L)-1; Q=L[1][1]; if (nummonoms==0){return(f);} for (i=2;i<=size(L);i++) { CC=matrix(extendcoef(L[i][1],Q,idp,idq)); NewMonoms[i-1]=list(CC,L[i][2]); } if (nummonoms==1) { for(j=1;j<=ncols(NewMonoms[1][1]);j++) { fout=NewMonoms[1][1][2,j]*L[1][2]+NewMonoms[1][1][1,j]*NewMonoms[1][2]; //fout=pnormalf(fout,idp,W); if(ncols(NewMonoms[1][1])>1){idout[j]=fout;} } if(ncols(NewMonoms[1][1])==1) { return(fout); } else { return(idout); } } else { list cfi; list coefs; for (i=1;i<=nummonoms;i++) { kill cfi; list cfi; for(j=1;j<=ncols(NewMonoms[i][1]);j++) { cfi[size(cfi)+1]=NewMonoms[i][1][2,j]; } coefs[i]=cfi; } def indexpolys=findindexpolys(coefs); for(i=1;i<=size(indexpolys);i++) { fout=L[1][2]; for(j=1;j<=nummonoms;j++) { fout=fout+(NewMonoms[j][1][1,indexpolys[i][j]])/(NewMonoms[j][1][2,indexpolys[i][j]])*NewMonoms[j][2]; } fout=cleardenom(fout); if(size(indexpolys)>1){idout[i]=fout;} } if (size(indexpolys)==1) { return(fout); } else { return(idout); } } } example { echo = 2; "EXAMPLE:"; if(defined(R)){kill R;} ring R=(0,a1,a2),(x),lp; short=0; poly f=(a1^2-4*a1+a2^2-a2)*x+(a1^4-16*a1+a2^3-4*a2); ideal p=a1*a2; ideal q=a2^2-a2,a1*a2,a1^2-4*a1; extendpoly(f,p,q);" "; // EXAMPLE; if (defined(R)){kill R;} ring R=(0,a0,b0,c0,a1,b1,c1,a2,b2,c2),(x), dp; short=0; poly f=(b1*a2*c2-c1*a2*b2)*x+(-a1*c2^2+b1*b2*c2+c1*a2*c2-c1*b2^2); ideal p= (-a0*b1*c2+a0*c1*b2+b0*a1*c2-b0*c1*a2-c0*a1*b2+c0*b1*a2), (a1^2*c2^2-a1*b1*b2*c2-2*a1*c1*a2*c2+a1*c1*b2^2+b1^2*a2*c2-b1*c1*a2*b2+c1^2*a2^2), (a0*a1*c2^2-a0*b1*b2*c2-a0*c1*a2*c2+a0*c1*b2^2+b0*b1*a2*c2-b0*c1*a2*b2 - c0*a1*a2*c2+c0*c1*a2^2), (a0^2*c2^2-a0*b0*b2*c2-2*a0*c0*a2*c2+a0*c0*b2^2+b0^2*a2*c2-b0*c0*a2*b2+c0^2*a2^2), (a0*a1*c1*c2-a0*b1^2*c2+a0*b1*c1*b2-a0*c1^2*a2+b0*a1*b1*c2-b0*a1*c1*b2 -c0*a1^2*c2+c0*a1*c1*a2), (a0^2*c1*c2-a0*b0*b1*c2-a0*c0*a1*c2+a0*c0*b1*b2-a0*c0*c1*a2+b0^2*a1*c2 -b0*c0*a1*b2+c0^2*a1*a2), (a0^2*c1^2-a0*b0*b1*c1-2*a0*c0*a1*c1+a0*c0*b1^2+b0^2*a1*c1-b0*c0*a1*b1+c0^2*a1^2), (2*a0*a1*b1*c1*c2-a0*a1*c1^2*b2-a0*b1^3*c2+a0*b1^2*c1*b2-a0*b1*c1^2*a2 -b0*a1^2*c1*c2+b0*a1*b1^2*c2-b0*a1*b1*c1*b2+b0*a1*c1^2*a2-c0*a1^2*b1*c2+c0*a1^2*c1*b2); ideal q= (-a1*c2+c1*a2), (-a1*b2+b1*a2), (-a0*c2+c0*a2), (-a0*b2+b0*a2), (-a0*c1+c0*a1), (-a0*b1+b0*a1), (-a1*b1*c2+a1*c1*b2), (-a0*b1*c2+a0*c1*b2), (-a0*b0*c2+a0*c0*b2), (-a0*b0*c1+a0*c0*b1); extendpoly(f,p,q); } // if L is a list(ideal,ideal) return 1 else returns 0; static proc typeofCrep(L) { if(typeof(L)!="list"){return(0);} if(size(L)!=2){return(0);} if((typeof(L[1])!="ideal") or (typeof(L[2])!="ideal")){return(0);} return(1); } // Input. GC the grobcov of an ideal in generic representation of the // bases computed with option option ("rep",2). // Output The grobcov in full representation. // Option ("comment",1) shows the time. // Can be called from the top proc extendGC(list GC) "USAGE: extendGC(list GC); list GC must the grobcov of a parametric ideal computed with option \"rep\",2. It determines the full representation. The default option of grobcov provides the bases in generic representation (the I-regular functions forming the bases are then given by a single polynomial. They can specialize to zero for some points of the segments, but in general, it is sufficient for many purposes. Nevertheless the I-regular functions allow a full representation given by a set of polynomials specializing to the value of the function (after normalization) or to zero, but at least one of the polynomials specializes to non-zero. The full representation can be obtained by computing the grobcov with option \"ext\",1. (The default option there is \"ext\",0). With option \"ext\",1 the computation can be much more time consuming, but the result can be simpler. Alternatively, one can compute the full representation of the bases after computing grobcov with the default option for \"ext\" and the option \"rep\",2, that outputs both the Prep and the Crep of the segments, and then call \"extendGC\" to its output. RETURN: When calling extendGC(grobcov(S,\"rep\",2)) the result is of the form [[[lpp_1,basis_1,segment_1,lpph_1], ... , [lpp_s,basis_s,segment_s,lpph_s]] ], where each function of the basis can be given by an ideal of representants. NOTE: The basering R, must be of the form Q[a][x], (a=parameters, x=variables), and should be defined previously. The ideal must be defined on R. KEYWORDS: Groebner cover; parametric ideal; canonical, discussion of parametric ideal; full representation EXAMPLE: extendGC; shows an example" { int te; def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); Rx[1]=0; def D=ring(Rx); def RP=D+P; list S=GC; ideal idp; ideal idq; int i; int j; int m; int s; int k; m=0; i=1; while((i<=size(S)) and (m==0)) { if(typeof(S[i][2])=="list"){m=1;} i++; } if(m==1) { "Warning! grobcov has already extended bases"; return(S); } if(typeofCrep(S[1][3])){k=3;} else{if(typeofCrep(S[1][4])){k=4;};} if(k==0) { "Warning! extendGC make sense only when grobcov has been called with option 'rep',1 or 'rep',2"; // if(te==0){kill @R; kill @RP; kill @P;} return(S); } poly leadc; poly ext; list SS; // Now extendGC for (i=1;i<=size(S);i++) { m=size(S[i][2]); for (j=1;j<=m;j++) { idp=S[i][k][1]; idq=S[i][k][2]; if (size(idp)>0) { leadc=leadcoef(S[i][2][j]); kill ext; def ext=extendpoly(S[i][2][j],idp,idq); if (typeof(ext)=="poly") { S[i][2][j]=pnormalf(ext,idp,idq); } else { if(size(ext)==1) { S[i][2][j]=ext[1]; } else { kill SS; list SS; for(s=1;s<=size(ext);s++) { ext[s]=pnormalf(ext[s],idp,idq); } for(s=1;s<=size(S[i][2]);s++) { if(s!=j){SS[s]=S[i][2][s];} else{SS[s]=ext;} } S[i][2]=SS; } } } } } return(S); } example { "EXAMPLE:"; echo = 2; if(defined(R)){kill R;} ring R=(0,a0,b0,c0,a1,b1,c1),(x), dp; short=0; ideal S=a0*x^2+b0*x+c0, a1*x^2+b1*x+c1; def GCS=grobcov(S,"rep",2); // grobcov(S) with both P and C representations GCS; def FGC=extendGC(GCS,"rep",1); // Full representation FGC; } // Auxiliary routine // nonzerodivisor // input: // poly g in Q[a], // list P=(p_1,..p_r) representing a minimal prime decomposition // output // poly f such that f notin p_i for all i and // g-f in p_i for all i such that g notin p_i static proc nonzerodivisor(poly gr, list Pr) { def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); setring(P); def g=imap(RR,gr); def P=imap(RR,Pr); int i; int k; list J; ideal F; def f=g; ideal Pi; for (i=1;i<=size(P);i++) { option(redSB); Pi=std(P[i]); //attrib(Pi,"isSB",1); if (reduce(g,Pi,5)==0){J[size(J)+1]=i;} } for (i=1;i<=size(J);i++) { F=ideal(1); for (k=1;k<=size(P);k++) { if (k!=J[i]) { F=idint(F,P[k]); } } f=f+F[1]; } setring(RR); def fr=imap(P,f); return(fr); } //Auxiliary routine // nullin // input: // poly f: a polynomial in Q[a] // ideal P: an ideal in Q[a] // called from ring @R // output: // t: with value 1 if f reduces modulo P, 0 if not. static proc nullin(poly f,ideal P) { int t; def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); setring(P); def f0=imap(RR,f); def P0=imap(RR,P); attrib(P0,"isSB",1); if (reduce(f0,P0,5)==0){t=1;} else{t=0;} setring(RR); return(t); } // Auxiliary routine // monoms // Input: A polynomial f // Output: The list of leading terms static proc monoms(poly f) { list L; poly lm; poly lc; poly lp; poly Q; poly mQ; def p=f; int i=1; while (p!=0) { lm=lead(p); p=p-lm; lc=leadcoef(lm); lp=leadmonom(lm); L[size(L)+1]=list(lc,lp); i++; } return(L); } // Auxiliary routine // findindexpolys // input: // list coefs=( (q11,..,q1r_1),..,(qs1,..,qsr_1) ) // of denominators of the monoms // output: // list ind=(v_1,..,v_t) of intvec // each intvec v=(i_1,..,is) corresponds to a polynomial in the sheaf // that will be built from it in extend procedures. static proc findindexpolys(list coefs) { int i; int j; intvec numdens; for(i=1;i<=size(coefs);i++) { numdens[i]=size(coefs[i]); } // def RR=basering; // def Rx=ringlist(RR); // def P=ring(Rx[1]); // setring(P); // def coefsp=imap(RR,coefs); def coefsp=coefs; ideal cof; list combpolys; intvec v; int te; list mp; for(i=1;i<=size(coefsp);i++) { cof=ideal(0); for(j=1;j<=size(coefsp[i]);j++) { cof[j]=factorize(coefsp[i][j],3); } coefsp[i]=cof; } for(j=1;j<=size(coefsp[1]);j++) { v[1]=j; te=1; for (i=2;i<=size(coefsp);i++) { mp=memberpos(coefsp[1][j],coefsp[i]); if(mp[1]) { v[i]=mp[2]; } else{v[i]=0;} } combpolys[j]=v; } combpolys=reform(combpolys,numdens); //"T_combpolys="; combpolys; //setring(RR); //def combpolysT=imap(P,combpolys); // return(combpolysT); return(combpolys); } // Auxiliary routine // extendcoef: given Q,P in Q[a] where P/Q specializes on an open and dense subset // of the whole V(p1 int...int pr), it returns a basis of the module // of all syzygies equivalent to P/Q, static proc extendcoef(poly fP, poly fQ, ideal idp, ideal idq) { def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); setring(P); def PL=ringlist(P); PL[3][1][1]="dp"; def P1=ring(PL); setring(P1); ideal idp0=imap(RR,idp); option(redSB); qring q=std(idp0); poly P0=imap(RR,fP); poly Q0=imap(RR,fQ); ideal PQ=Q0,-P0; module C=syz(PQ); setring(P); def idp1=imap(RR,idp); def idq1=imap(RR,idq); def C1=matrix(imap(q,C)); def redC=selectregularfun(C1,idp1,idq1); setring(RR); def CC=imap(P,redC); return(CC); } // Auxiliary routine // selectregularfun // input: // list L of the polynomials matrix CC // (we assume that one of them is non-null on V(N) \ V(M)) // ideal N, ideal M: ideals representing the locally closed set V(N) \ V(M) // assume to work in @P static proc selectregularfun(matrix C, ideal N, ideal M) { int numcombused; // def RR=basering; // def Rx=ringlist(RR); // def P=ring(Rx[1]); // setring(P); // def C=imap(RR,CC); // def N=imap(RR,NN); // def M=imap(RR,MM); if (ncols(C)==1){return(C);} int i; int j; int k; list c; intvec ci; intvec c0; intvec c1; list T; list T0; list T1; list LL; ideal N1;ideal M1; int te=0; for(i=1;i<=ncols(C);i++) { if((C[1,i]!=0) and (C[2,i]!=0)) { if(c0==intvec(0)){c0[1]=i;} else{c0[size(c0)+1]=i;} } } def C1=submat(C,1..2,c0); for (i=1;i<=ncols(C1);i++) { c=comb(ncols(C1),i); for(j=1;j<=size(c);j++) { ci=c[j]; numcombused++; if(i==1){N1=N+C1[2,j]; M1=M;} if(i>1) { kill c0; intvec c0 ; kill c1; intvec c1; c1=ci[size(ci)]; for(k=1;k0) { if(notfree[1]==0){notfree[1]=combpolys[j][i];} else{notfree[size(notfree)+1]=combpolys[j][i];} } } kill free1; intvec free1; for(j=1;j<=numdens[i];j++) { if(memberpos(j,notfree)[1]==0) { if(free1[1]==0){free1[1]=j;} else{free1[size(free1)+1]=j;} } free[i]=free1; } } list amplcombp; list aux; for(i=1;i<=size(combp0);i++) { v=combp0[i]; kill amplcombp; list amplcombp; amplcombp[1]=intvec(v[1]); for(j=2;j<=size(v);j++) { if(v[j]!=0) { for(k=1;k<=size(amplcombp);k++) { w=amplcombp[k]; w[size(w)+1]=v[j]; amplcombp[k]=w; } } else { kill aux; list aux; for(k=1;k<=size(amplcombp);k++) { for(l=1;l<=size(free[j]);l++) { w=amplcombp[k]; w[size(w)+1]=free[j][l]; aux[size(aux)+1]=w; } } amplcombp=aux; } } for(j=1;j<=size(amplcombp);j++) { combp1[size(combp1)+1]=amplcombp[j]; } } return(combp1); } // Auxiliary routine // precombint // input: L: list of ideals (works in @P) // output: F0: ideal of polys. F0[i] is a poly in the intersection of // all ideals in L except in the ith one, where it is not. // L=(p1,..,ps); F0=(f1,..,fs); // F0[i] \in intersect_{j#i} p_i static proc precombint(list L) { int i; int j; int tes; def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); setring(P); list L0; list L1; list L2; list L3; ideal F; L0=imap(RR,L); L1[1]=L0[1]; L2[1]=L0[size(L0)]; for (i=2;i<=size(L0)-1;i++) { L1[i]=intersect(L1[i-1],L0[i]); L2[i]=intersect(L2[i-1],L0[size(L0)-i+1]); } L3[1]=L2[size(L2)]; for (i=2;i<=size(L0)-1;i++) { L3[i]=intersect(L1[i-1],L2[size(L0)-i]); } L3[size(L0)]=L1[size(L1)]; for (i=1;i<=size(L3);i++) { option(redSB); L3[i]=std(L3[i]); } for (i=1;i<=size(L3);i++) { tes=1; j=0; while((tes) and (j0) { for (i=1;i<=size(L) div 2;i++) { if (L[2*i-1]=="null"){E=L[2*i];} else { if (L[2*i-1]=="nonnull"){N=L[2*i];} else { if (L[2*i-1]=="comment"){comment=L[2*i];} else { if (L[2*i-1]=="out"){out=L[2*i];} } } } } } if (comment>0){string("Begin KSW with null = ",E," nonnull = ",N);} def CG=KSW0(F,E,N,comment); if (comment>0) { string("Number of segments in KSW (total) = ",size(CG)); string("Time in KSW = ",timer-start); } if(out==0) { CG=KSWtocgsdr(CG); //"T_CG="; CG; if( size(CG)>0) { CG=groupKSWsegments(CG); if (comment>0) { string("Number of lpp segments = ",size(CG)); string("Time in KSW + group + Prep = ",timer-start); } } } return(CG); } // Auxiliary routine // sqf static proc sqf(poly f) { def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); setring(P); def ff=imap(RR,f); poly fff=sqrfree(ff,3); setring(RR); def ffff=imap(P,fff); return(ffff); } // Auxiliary routine // KSW0: Kapur-Sun-Wang algorithm for computing a CGS, called by KSW // Input: // F: parametric ideal to be discussed // Options: // The ideal must be defined on C[parameters][variables] // Output: static proc KSW0(ideal F, ideal E, ideal N, int comment) { def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); Rx[1]=0; def D=ring(Rx); def RP=D+P; int i; int j; list emp; list CGS; ideal N0; for (i=1;i<=size(N);i++) { N0[i]=sqf(N[i]); } ideal E0; for (i=1;i<=size(E);i++) { E0[i]=sqf(leadcoef(E[i])); } setring(P); ideal E1=imap(RR,E0); E1=std(E1); ideal N1=imap(RR,N0); N1=std(N1); setring(RR); E0=imap(P,E1); N0=imap(P,N1); if (inconsistent(E0,N0)==1) { return(emp); } setring(RP); def FRP=imap(RR,F); def ERP=imap(RR,E); FRP=FRP+ERP; option(redSB); def GRP=std(FRP); setring(RR); def G=imap(RP,GRP); if (memberpos(1,G)[1]==1) { if(comment>1){"Basis 1 is found"; E; N;} list KK; KK[1]=list(E0,N0,ideal(1)); return(KK); } ideal Gr; ideal Gm; ideal GM; for (i=1;i<=size(G);i++) { if (variables(G[i])[1]==0){Gr[size(Gr)+1]=G[i];} else{Gm[size(Gm)+1]=G[i];} } ideal Gr0; for (i=1;i<=size(Gr);i++) { Gr0[i]=sqf(Gr[i]); } Gr=elimrepeated(Gr0); ideal GrN; for (i=1;i<=size(Gr);i++) { for (j=1;j<=size(N0);j++) { GrN[size(GrN)+1]=sqf(Gr[i]*N0[j]); } } if (inconsistent(E,GrN)){;} else { if(comment>1){"Basis 1 is found in a branch with arguments"; E; GrN;} CGS[size(CGS)+1]=list(E,GrN,ideal(1)); } if (inconsistent(Gr,N0)){return(CGS);} GM=Gm; Gm=MDBasis(Gm); ideal H; for (i=1;i<=size(Gm);i++) { H[i]=sqf(leadcoef(Gm[i])); } H=facvar(H); poly h=sqf(LCMLC(H)); if(comment>1){"H = "; H; "h = "; h;} ideal Nh=N0; if(size(N0)==0){Nh=h;} else { for (i=1;i<=size(N0);i++) { Nh[i]=sqf(N0[i]*h); } } if (inconsistent(Gr,Nh)){;} else { CGS[size(CGS)+1]=list(Gr,Nh,Gm); } poly hc=1; list KS; ideal GrHi; for (i=1;i<=size(H);i++) { kill GrHi; ideal GrHi; Nh=N0; if (i>1){hc=sqf(hc*H[i-1]);} for (j=1;j<=size(N0);j++){Nh[j]=sqf(N0[j]*hc);} if (equalideals(Gr,ideal(0))==1){GrHi=H[i];} else {GrHi=Gr,H[i];} if(comment>1){"Call to KSW with arguments "; GM; GrHi; Nh;} KS=KSW0(GM,GrHi,Nh,comment); for (j=1;j<=size(KS);j++) { CGS[size(CGS)+1]=KS[j]; } if(comment>1){"CGS after KSW = "; CGS;} } return(CGS); } // Auxiliary routine // KSWtocgsdr static proc KSWtocgsdr(list L) { int i; list CG; ideal B; ideal lpp; int j; list NKrep; for(i=1;i<=size(L);i++) { B=redgbn(L[i][3],L[i][1],L[i][2]); lpp=ideal(0); for(j=1;j<=size(B);j++) { lpp[j]=leadmonom(B[j]); } NKrep=KtoPrep(L[i][1],L[i][2]); CG[i]=list(lpp,B,NKrep); } return(CG); } // Auxiliary routine // KtoPrep // Computes the P-representaion of a K-representation (N,W) of a set // input: // ideal E (null conditions) // ideal N (non-null conditions ideal) // output: // the ((p_1,(p_11,..,p_1k_1)),..,(p_r,(p_r1,..,p_rk_r))); // the Prep of V(N) \ V(W) static proc KtoPrep(ideal N, ideal W) { int i; int j; if (N[1]==1) { L0[1]=list(ideal(1),list(ideal(1))); return(L0); } def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); setring(P); ideal B; int te; poly f; ideal Np=imap(RR,N); ideal Wp=imap(RR,W); list L; list L0; list T0; L0=minGTZ(Np); for(j=1;j<=size(L0);j++) { option(redSB); L0[j]=std(L0[j]); } for(i=1;i<=size(L0);i++) { if(inconsistent(L0[i],Wp)==0) { B=L0[i]+Wp; T0=minGTZ(B); option(redSB); for(j=1;j<=size(T0);j++) { T0[j]=std(T0[j]); } L[size(L)+1]=list(L0[i],T0); } } setring(RR); def LL=imap(P,L); return(LL); } // Auxiliary routine // groupKSWsegments // input: the list of vertices of KSW // output: the same terminal vertices grouped by lpp static proc groupKSWsegments(list T) { int i; int j; list L; list lpp; list lppor; list kk; lpp[1]=T[1][1]; j=1; lppor[1]=intvec(1); for(i=2;i<=size(T);i++) { kk=memberpos(T[i][1],lpp); if(kk[1]==0){j++; lpp[j]=T[i][1]; lppor[j]=intvec(i);} else{lppor[kk[2]][size(lppor[kk[2]])+1]=i;} } list ll; for (j=1;j<=size(lpp);j++) { kill ll; list ll; for(i=1;i<=size(lppor[j]);i++) { ll[size(ll)+1]=list(i,T[lppor[j][i]][2],T[lppor[j][i]][3]); } L[j]=list(lpp[j],ll,string(lpp[j])); } return(L); } //********************* End KapurSunWang ************************* //********************* Begin ConsLevels *************************** static proc zeroone(int n) { list L; list L2; intvec e; intvec e2; intvec e3; int j; if(n==1) { e[1]=0; L[1]=e; e[1]=1; L[2]=e; return(L); } if(n>1) { L=zeroone(n-1); for(j=1;j<=size(L);j++) { e2=L[j]; e3=e2; e3[size(e3)+1]=0; L2[size(L2)+1]=e3; e3=e2; e3[size(e3)+1]=1; L2[size(L2)+1]=e3; } } return(L2); } // Auxiliary routine // subsets: the list of subsets of (1,..n) static proc subsets(int n) { list L; list L1; int i; int j; L=zeroone(n); intvec e; intvec e1; for(i=1;i<=size(L);i++) { e=L[i]; kill e1; intvec e1; for(j=1;j<=n;j++) { if(e[n+1-j]==1) { if(e1==intvec(0)){e1[1]=j;} else{e1[size(e1)+1]=j}; } } L1[i]=e1; } return(L1); } // Input a list A of locally closed sets in C-rep // Output a list B of a simplified list of A static proc SimplifyUnion(list A) { int i; int j; list L=A; int n=size(L); if(n<2){return(A);} intvec w; for(i=1;i<=size(L);i++) { for(j=1;j<=size(L);j++) { if(i != j) { if(equalideals(L[i][2],L[j][1])==1) { L[i][2]=L[j][2]; w[size(w)+1]=j; } } } } if(size(w)>0) { for(i=1; i<=size(w);i++) { j=w[size(w)+1-i]; L=elimfromlist(L, j); } } ideal T=ideal(1); intvec v; for(i=1;i<=size(L);i++) { if(equalideals(L[i][2],ideal(1))) { v[size(v)+1]=i; T=intersect(T,L[i][1]); } } if(size(v)>0) { for(i=1; i<=size(v);i++) { j=v[size(v)+1-i]; L=elimfromlist(L, j); } } if(equalideals(T,ideal(1))==0){L[size(L)+1]=list(std(T),ideal(1))}; return(L); } // input list A=[[p1,q1],...,[pn,qn]] : // the list of segments of a constructible set S, where each [pi,qi] is given in C-representation // output list [topA,C] // where topA is the closure of A // C is the list of segments of the complement of A given in C-representation static proc FirstLevel(list A) { int n=size(A); list T=zeroone(n); ideal P; ideal Q; list Cb; ideal Cc=1; int i; int j; intvec t; ideal topA=1; list C; for(i=1;i<=n;i++) { topA=intersect(topA,A[i][1]); } //topA=std(topA); for(i=2; i<=size(T);i++) { t=T[i]; //"T_t"; t; P=0; Q=1; for(j=1;j<=n;j++) { if(t[n+1-j]==1) { P=P+A[j][2]; } else { Q=intersect(Q,A[j][1]); } } Cb=Crep0(P,Q); //"T_Cb="; Cb; if(size(Cb)!=0) { if( Cb[1][1]<>1) { C[size(C)+1]=Cb; } } } if(size(C)>1){C=SimplifyUnion(C);} return(list(topA,C)); } // Input: // Output: static proc ConstoPrep(list L) { list L1; int i; int j; list aux; for(i=1;i<=size(L);i++) { aux=Prep(L[i][2][1],L[i][2][2]); L1[size(L1)+1]=list(L[i][1],aux); } return(L1); } // Input: // list A = [[P1,Q1], .. [Pn,Qn]] // A constructible set as union of locally closed sets represented by pairs of ideals // Output: // list L =[p1,p2,p3,...,pk] // where pi is the ideal of the closure of level i alternatively of A or its complement // Note: the levels of A are [p1,p2], [p3,p4], [p5,p6],... // the levels of C are [p2,p3],[p4,p5], ... // expressed in C-representation // Assumed to be called in the ring Q[a] proc ConsLevels(list A0) "USAGE: ConsLevels(list L); L=[[P1,Q1],...,[Ps,Qs]] is a list of lists of of pairs of ideals represening the constructible set S=V(P1) \ V(Q1) u ... u V(Ps) \ V(Qs). To be called in a ring Q[a][x] or a ring Q[a]. But the ideals can contain only the parameters in Q[a]. RETURN: The list of ideals [a1,a2,...,at] representing the closures of the canonical levels of S and its complement C wrt to the closure of S. The canonical levels of S are represented by theirs Crep. So we have: Levels of S: [a1,a2],[a3,a4],... Levels of C: [a2,a3],[a4,a5],... S=V(a1) \ V(a2) u V(a3) \ V(a4) u ... C=V(a2 \ V(a3) u V(a4) \ V(a5) u ... The expression of S can be obtained from the output of ConsLevels by the call to Levels. NOTE: The algorithm was described in J.M. Brunat, A. Montes. \"Computing the canonical representation of constructible sets.\" Math. Comput. Sci. (2016) 19: 165-178. KEYWORDS: constructible set; locally closed set; canonical form EXAMPLE: ConsLevels; shows an example" { int te; def RR=basering; def Rx=ringlist(RR); if(size(Rx[1])==4) { te=1; def P=ring(Rx[1]); setring P; list A=imap(RR,A0); } // if(defined(@P)){te=1; setring(@P); list A=imap(RR,A0);} else {te=0; def A=A0;} list L; list C; list B; list T; int i; for(i=1; i<=size(A);i++) { T=Crep0(A[i][1],A[i][2]); B[size(B)+1]=T; } list K; while(size(B)>0) { K=FirstLevel(B); //"T_K="; K; L[size(L)+1]=K[1]; B=K[2]; } L[size(L)+1]=ideal(1); if(te==1) {setring(RR); def LL=imap(P,L);} if(te==0){def LL=L;} return(LL); } example { "EXAMPLE:"; echo = 2; if(defined(R)){kill R;} ring R=0,(x,y,z),lp; short=0; ideal P1=(x^2+y^2+z^2-1); ideal Q1=z,x^2+y^2-1; ideal P2=y,x^2+z^2-1; ideal Q2=z*(z+1),y,x*(x+1); ideal P3=x; ideal Q3=5*z-4,5*y-3,x; list Cr1=Crep(P1,Q1); list Cr2=Crep(P2,Q2); list Cr3=Crep(P3,Q3); list L=list(Cr1,Cr2,Cr3); L; def LL=ConsLevels(L); LL; def LLL=Levels(LL); LLL; } // Converts the output of ConsLevels, given by the set of closures of the Levels of the constructible S // to an expression where the Levels are apparent. // Input: The ouput of ConsLevels of the form // [A1,A2,..,Ak], where the Ai's are the closures of the levels. // Output: An expression of the form // L1=[[1,[A1,A2]],[3,[A3,A4]],..,[2l-1,[A_{2l-1},A_{2l}]]] the list of Levels of S proc Levels(list L) "USAGE: Levels(list L); The input list L must be the output of the call to the routine ConsLevels of a constructible set: L=[a1,a2,..,ak], where the a's are the closures of the levels, determined by ConsLevels. Levels selects the levels of the constructible set. To be called in a ring Q[a][x] or a ring Q[a]. But the ideals can contain only the parameters in Q[a]. RETURN: The levels of the constructible set: Lc=[ [1,[a1,a2]],[3,[a3,a4]],.., [2l-1,[a_{2l-1},a_{2l}]] ] the list of levels of S KEYWORDS: constructible sets; canonical form EXAMPLE: Levels shows an example" { int n=size(L) div 2; int i; list L1; list L2; for(i=1; i<=n;i++) { L1[size(L1)+1]=list(2*i-1,list(L[2*i-1],L[2*i])); } return(L1); } example { "EXAMPLE:"; echo = 2; if(defined(R)){kill R;} ring R=0,(x,y,z),lp; short=0; ideal P1=(x^2+y^2+z^2-1); ideal Q1=z,x^2+y^2-1; ideal P2=y,x^2+z^2-1; ideal Q2=z*(z+1),y,x*(x+1); ideal P3=x; ideal Q3=5*z-4,5*y-3,x; list Cr1=Crep(P1,Q1); list Cr2=Crep(P2,Q2); list Cr3=Crep(P3,Q3); list L=list(Cr1,Cr2,Cr3); L; def LL=ConsLevels(L); LL; def LLL=Levels(LL); LLL; } proc DifConsLCSets(list A, list B) "USAGE: DifConsLCSets(list A,list B); Input: The input lists A and B must be each one the canonical representations of the respective constructible sets, i.e. outputs of the routine ConsLevels for a constructible set, or from the routine Grob1Levels applied to the output of grobcov. A=[a1,a2,..,ak], B=[b1,b2,..,bj], where the a's and the b's are the closures of the levels of the constructible and the complements determined by ConsLevels (or GrobLevels) To be called in a ring Q[a][x] or a ring Q[a]. But the ideals can contain only the parameters in Q[a]. RETURN: A list of locally closed sets equivalent to the difference S= A "\" B. Lc=[ [1][p1,q1]] [[2][p2,q2]]..], For obtaining the canonical representation into levels of the constructible A "\" B one have to apply ConsLevels and then optatively Levels. KEYWORDS: constructible sets; canonical form EXAMPLE: DifConsLCSets shows an example" { int n; int m; int t; int i; int j; int k; ideal ABup; ideal ABdw; if (size(B) mod 2==1){B[size(B)+1]=ideal(1);} if (size(A) mod 2==1){A[size(A)+1]=ideal(1);} //"T_A=";A; // "T_B="; B; n=size(A) div 2; m=(size(B) div 2)-1; //string("T_n=",n); //string("T_m=",m); list L; list M; list ABupC; //list LL; for(i=1;i<=n;i++) { //string("T_i=",i); t=1; j=0; //list L; while (t==1 and j<=m) { //string("T_j=",j); ABdw=intersectpar(list(A[2*i],B[2*j+1])); //"T_ABdw="; ABdw; ABup=A[2*i-1]; //"T_ABup1="; ABup; if(j>0) { for(k=1;k<=size(B[2*j]);k++) { ABup[size(ABup)+1]=B[2*j][k]; } } //"T_ABup2="; ABup; ABupC=Crep(ABup,ideal(1)); //"T_ABupC="; ABupC; ABup=ABupC[1]; //"T_ABup="; ABup; if(ABup==1){t=0;} //if(equalideals(ABup,ideal(1))){t=0;} else { M=Crep(ABup,ABdw); //"T_M="; M; //if(not(equalideals(M[1],ideal(1)))) {L[size(L)+1]=M;} if(not(size(M)==0)) {L[size(L)+1]=M;} } //"L="; L; j++; } //LL[size(LL)+1]=L; } return(L); } example { "EXAMPLE:"; echo = 2; if(defined(R)){kill R;} ring R=(0,x,y,z,t),(x1,y1),lp; ideal a1=x; ideal a2=x,y; ideal a3=x,y,z; ideal a4=x,y,z,t; ideal b1=y; ideal b2=y,z; ideal b3=y,z,t; ideal b4=1; list L1=a1,a2,a3,a4; list L2=b1,b2,b3,b4; L1; L2; def LL=DifConsLCSets(L1,L2); LL; def LLL=ConsLevels(LL); LLL; def LLLL=Levels(LLL); LLLL; } //**************************** End ConstrLevels ****************** //******************** Begin locus and envelop ****************************** // indepparameters // Auxiliary routine to detect 'Special' components of the locus // Input: ideal B // Output: // 1 if the ideal does not depend on the parameters // 0 if they depend static proc indepparameters(ideal B) { def RR=basering; list Rx=ringlist(RR); def P=ring(Rx[1]); Rx[1]=0; def D=ring(Rx); def RP=D+P; // if(defined(@P)){kill @P; kill @RP; kill @R;} // setglobalrings(); ideal v=variables(B); setring RP; def BP=imap(RR,B); def vp=imap(RR,v); ideal varpar=variables(BP); int te; te=equalideals(vp,varpar); setring(RR); // kill @P; kill @RP; kill @R; if(te){return(1);} else{return(0);} } // indepparameterspoly // Auxiliary routine to detect 'Special' components of the locus // Input: ideal B // Output: // 1 if the solutions of the ideal (or poly) does not depend on the parameters // 0 if they depend static proc indepparameterspoly(B) { def RR=basering; list Rx=ringlist(RR); def P=ring(Rx[1]); Rx[1]=0; def D=ring(Rx); def RP=D+P; // if(defined(@P)){kill @P; kill @RP; kill @R;} // setglobalrings(); ideal v=variables(B); setring RP; def BP=imap(RR,B); def vp=imap(RR,v); ideal varpar=variables(BP); int te; te=equalideals(vp,varpar); setring(RR); // kill @P; kill @RP; kill @R; if(te){return(1);} else{return(0);} } // dimP0: Auxiliary routine // if the dimension in @P of an ideal in the parameters has dimension 0 then it returns 0 // else it retuns 1 static proc dimP0(ideal N) { def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); setring P; // if(defined(@P)){ kill @P; kill @RP; kill @R;} // setglobalrings(); // setring(@P); int te=1; def NP=imap(RR,N); attrib(NP,"IsSB",1); int d=dim(std(NP)); //"T_d="; d; if(d==0){te=0;} setring(RR); return(te); } // DimPar // Auxilliary routine to NorSing determining the dimension of a parametric ideal // Does not use @P and define it directly because is assumes that // variables and parameters have been inverted static proc DimPar(ideal E) { def RRH=basering; def RHx=ringlist(RRH); def P=ring(RHx[1]); setring(P); def E2=std(imap(RRH,E)); attrib(E2,"IsSB",1); def d=dim(E2); setring RRH; return(d); } // DimPar // Auxilliary routine to locus2 determining the dimension of a parametric ideal // it is identical to DimPar but adds infromation about the character of the component static proc dimComp(ideal PA) { def RR=basering; list Rx=ringlist(RR); int n=size(Rx[1][2]); def P=ring(Rx[1]); setring(P); list Lout; def CP=imap(RR,PA); attrib(CP,"IsSB",1); int d=dim(std(CP)); if(d==n-1){Lout[1]=d;Lout[2]="Degenerate"; } else {Lout[1]=d; Lout[2]="Accumulation";} setring RR; return(Lout); } // Takes a list of intvec and sorts it and eliminates repeated elements. // Auxiliary routine static proc sortpairs(L) { def L1=sort(L); def L2=elimrepeated(L1[1]); return(L2); } // Eliminates the pairs of L1 that are also in L2. // Auxiliary routine static proc minuselements(list L1,list L2) { int i; list L3; for (i=1;i<=size(L1);i++) { if(not(memberpos(L1[i],L2)[1])){L3[size(L3)+1]=L1[i];} } return(L3); } // NorSing // Input: // ideal B: the basis of a component of the grobcov // ideal E: the top of the component (assumed to be of dimension > 0 (single equation) // ideal N: the holes of the component // Output: // int d: the dimension of B on the variables (anti-image). // if d>0 then the component is 'Normal' // if d==0 then the component is 'Singular' static proc NorSing(ideal B, ideal E, ideal N, list #) { int i; int j; int Fenv=0; int env; int dd; list DD=#; def RR=basering; ideal vmov; int version=0; int nv=nvars(RR); int moverdim=nv; //npars(RR); if(nv<4){version=1;} int d; poly F; for(i=1;i<=(size(DD) div 2);i++) { if(DD[2*i-1]=="moverdim"){moverdim=DD[2*i];} if(DD[2*i-1]=="vmov"){vmov=DD[2*i];} if(DD[2*i-1]=="movdim"){moverdim=DD[2*i];} if(DD[2*i-1]=="version"){version=DD[2*i];} if(DD[2*i-1]=="family"){F=DD[2*i];} } if(F!=0){Fenv=1;} list L0; if(dimP0(E)==0){L0=2,"Normal";} // 2 es fals pero ha de ser >0 encara que sigui 0 else { if(version==0) { // Computing std(B+E,plex(x,y,x1,..xn)) one can detect if there is a first part // independent of parameters giving the variables with dimension 0 dd=indepparameters(B); //"T_dd="; dd; if (dd==1){d=0; L0=d,string(B);} // ,determineF(B,F,E) else{d=1; L0=2,"Normal";} } else { def RH=ringlist(RR); //"T_RH="; RH; def H=RH; H[1]=0; H[2]=RH[1][2]+RH[2]; int n=size(H[2]); intvec ll; for(i=1;i<=n;i++) { ll[i]=1; } H[3][1][1]="lp"; H[3][1][2]=ll; def RRH=ring(H); setring(RRH); ideal BH=imap(RR,B); ideal EH=imap(RR,E); ideal NH=imap(RR,N); if(Fenv==1){poly FH=imap(RR,F);} for(i=1;i<=size(EH);i++){BH[size(BH)+1]=EH[i];} BH=std(BH); // MOLT COSTOS!!! ideal G; ideal r; poly q; for(i=1;i<=size(BH);i++) { r=factorize(BH[i],1); q=1; for(j=1;j<=size(r);j++) { if((pdivi(r[j],NH)[1] != 0) or (equalideals(ideal(NH),ideal(1)))) { q=q*r[j]; } } if(q!=1){G[size(G)+1]=q;} } setring RR; def GG=imap(RRH,G); ideal GGG; if(defined(L0)){kill L0; list L0;} for(i=1;i<=size(GG);i++) { if(indepparameters(GG[i])){GGG[size(GGG)+1]=GG[i];} } GGG=std(GGG); ideal GLM; for(i=1;i<=size(GGG);i++) { GLM[i]=leadmonom(GGG[i]); } attrib(GLM,"IsSB",1); d=dim(std(GLM)); string antiimi=string(GGG); L0=d,antiimi; if(d==0) // d0){G2[size(G2)+1]=G[i];} } } if(size(G1)==0){t1=0;} // normal point components if(size(G2)==0){t2=0;} // non-normal point components setring(RP); if(t1) { list G1RP=imap(RR,G1); } else {list G1RP;} list P1RP; ideal B; for(i=1;i<=size(G1RP);i++) { kill B; ideal B; for(k=1;k<=size(G1RP[i][3]);k++) { attrib(G1RP[i][3][k][1],"IsSB",1); G1RP[i][3][k][1]=std(G1RP[i][3][k][1]); for(j=1;j<=size(G1RP[i][2]);j++) { B[j]=reduce(G1RP[i][2][j],G1RP[i][3][k][1]); } B=std(B); P1RP[size(P1RP)+1]=list(G1RP[i][3][k][1],G1RP[i][3][k][2],B); } } //In P1RP the basis has been reduced wrt the top setring(RR); ideal h; list P1; if(t1) { P1=imap(RP,P1RP); for(i=1;i<=size(P1);i++) { for(j=1;j<=size(P1[i][3]);j++) { h=factorize(P1[i][3][j],1); P1[i][3][j]=h[1]; for(k=2;k<=size(h);k++) { P1[i][3][j]=P1[i][3][j]*h[k]; } } } } //In P1 factors in the basis independent of the parameters have been obtained ideal BB; int dd; list NS; for(i=1;i<=size(P1);i++) { //"T_P1="; P1; NS=NorSing(P1[i][3],P1[i][1],P1[i][2][1],DD); //"T_NS="; NS; dd=NS[1]; if(dd==0){P1[i][3]=NS;} //"Special"; //dd==0 else{P1[i][3]="Normal";} } list P2; for(i=1;i<=size(G2);i++) { for(k=1;k<=size(G2[i][3]);k++) { P2[size(P2)+1]=list(G2[i][3][k][1],G2[i][3][k][2]); } } list l; for(i=1;i<=size(P1);i++){Q1[i]=l; Q1[i][1]=P1[i];} P1=Q1; for(i=1;i<=size(P2);i++){Q2[i]=l; Q2[i][1]=P2[i];} P2=Q2; // if(defined(@P)){te=1; kill @P; kill @RP; kill @R;} // setglobalrings(); setring(P); ideal J; if(t1==1) { def C1=imap(RR,P1); def L1=AddLocus(C1); } else{list C1; list L1; kill P1; list P1;} if(t2==1) { def C2=imap(RR,P2); def L2=AddLocus(C2); } else{list L2; list C2; kill P2; list P2;} for(i=1;i<=size(L2);i++) { J=std(L2[i][2]); d=dim(J); if(d+1==dimpar) { L2[i][4]=string("Degenerate",L2[i][4]); } else{L2[i][4]=string("Accumulation",L2[i][4]);} } list LN; if(t1==1) { for(i=1;i<=size(L1);i++){LN[size(LN)+1]=L1[i];} } if(t2==1) { for(i=1;i<=size(L2);i++){LN[size(LN)+1]=L2[i];} } int tLN=1; if(size(LN)==0){tLN=0;} setring(RR); if(tLN) { def L=imap(P,LN); for(i=1;i<=size(L);i++){if(size(L[i][2])==0){L[i][2]=ideal(1);}} //list LL; for(i=1;i<=size(L);i++) { if(typeof(L[i][4])=="list") {L[i][4][1]="Special";} l[1]=L[i][2]; l[2]=L[i][3]; l[3]=L[i][4]; l[4]=L[i][5]; L[i]=l; } } else{list L;} // if(te==1){kill @P; kill @RP; kill @R;} return(L); } // locus2(G): Private routine used by locus (the public routine), that // builds the diferent components. // input: The output G of the grobcov (in generic representation, which is the default option) // The ideal F defining the locus problem (G is the grobcov of F and has been // already determined by locus. // Options: The algorithm allows the following options as pair of arguments: // "moverdim", d : by default moverdim is the number of variabes but it can // be set to minor values. // "version", v : There are two versions of the algorithm. ('version',1) is // a full algorithm that always distinguishes correctly between 'Normal' // and 'Special' components, whereas ('version',0) can decalre a component // as 'Normal' being really 'Special', but is more effective. By default ('version',1) // is used when the number of variables is less than 4 and 0 if not. // The user can force to use one or other version, but it is not recommended. // "system", ideal F: if the initial systrem is passed as an argument. This is actually // not used. // "comments", c: by default it is 0, but it can be set to 1. // Usually locus problems have mover coordinates, variables and tracer coordinates. // The mover coordinates are to be placed as the last variables, and by default, // its number is equal to the number of parameters of tracer variables, but as option // it can be set to other values. // output: // list, the canonical P-representation of the Normal and Non-Normal locus: // The Normal locus has two kind of components: Normal and Special. // The Non-normal locus has two kind of components: Accumulation and Degenerate. // This routine is compemented by locus that calls it in order to eliminate problems // with degenerate points of the mover. // The output components are given as // ((p1,(p11,..p1s_1),type_1,level_1),..,(pk,(pk1,..pks_k),type_k,level_k) // The components are given in canonical P-representation of the subset. // If all levels of a class of locus are 1, then the set is locally closed. Otherwise the level // gives the depth of the component. static proc locus2(list G, ideal F, list #) { int st=timer; list Snor; list Snonor; int d; int i; int j; def RR=basering; def Rx=ringlist(RR); def RP=ring(Rx[1]); int tax=1; list DD=#; list GG=G; int n=size(Rx[1][2]); int moverdim; for(i=1;i<=(size(DD) div 2);i++) { if(DD[2*i-1]=="moverdim"){moverdim=DD[2*i];} if(DD[2*i-1]=="vmov"){vmov=DD[2*i];} // if(DD[2*i-1]=="vmov"){vmov=DD[2*i];} // if(DD[2*i-1]=="version"){version=DD[2*i];} // if(DD[2*i-1]=="comment"){comment=DD[2*i];} // if(DD[2*i-1]=="grobcov"){GG=DD[2*i];} if(DD[2*i-1]=="anti-image"){tax=DD[2*i];} } for(i=1;i<=size(GG);i++) { attrib(GG[i][1],"IsSB",1); GG[i][1]=std(GG[i][1]); d=dim(GG[i][1]); if(d==0) { for(j=1;j<=size(GG[i][3]);j++) { Snor[size(Snor)+1]=GG[i][3][j]; } } else { if(d>0) { for(j=1;j<=size(GG[i][3]);j++) { Snonor[size(Snonor)+1]=GG[i][3][j]; } } } } //"T_Snor="; Snor; //"T_Snonor="; Snonor; int tnor=size(Snor); int tnonor=size(Snonor); setring RP; list SnorP; list SnonorP; if(tnor) { SnorP=imap(RR,Snor); st=timer; SnorP=LCUnionN(SnorP); } if(tnonor) { SnonorP=imap(RR,Snonor); SnonorP=LCUnionN(SnonorP); } //"T_SnorP after LCUnion="; SnorP; // "T_SnonorP after LCUnion="; SnonorP; setring RR; ideal C; list N; list BAC; list AI; list NSC; list DAC; list L; ideal B; int k; int j0; int k0; int te; poly kkk=1; ideal AI0; int dimP; if(tnor) { Snor=imap(RP,SnorP); for(i=1;i<=size(Snor);i++) { C=Snor[i][1]; N=Snor[i][2]; dimP=DimPar(C); //"T_G="; G; AI=NS(F,G,C,N,DD); //"T_AI="; AI; AI0=AI[2]; if(size(AI)==3){AI[2]="normal"; AI[3]=ideal(0);} else { if(AI[1]>=dimP) { AI[2]="Normal"; if(tax>0){AI[3]=AI0;} } else { //"T_Al="; if(AI[1]=10){LK=n,ideal(0),"normal"; return(LK);} // 8 instead of 10 ? // Reducing basis B modulo C in the ring PR lex(x,u) setring(PR); def buPR=imap(RR,bu); def EE=imap(RR,E); ideal BR; def BB=imap(RR,B); def CC=imap(RR,C); attrib(CC,"IsSB",1); CC=std(CC); for(i=1;i<=size(BB);i++) { BR[i]=reduce(BB[i],CC); } // Computing the GB(B) inthe ring PR lex(x,u) BR=std(BR); ideal BRP; // Eliminating the polynomials in B containing variables v_t for(i=1;i<=size(BR);i++) { if(subset(variables(BR[i]),buPR )){BRP[size(BRP)+1]=BR[i];} } BR=BRP; // Factoring and eliminating the factors in N list H; ideal BP; poly f; int tj; //ideal vv=imap(RR,bu); for(i=1;i<=size(BR);i++) { f=1; H=factorize(BR[i]); if((subset(variables(H[1][2]),buPR)) and (size(H[1]))==2){BP[size(BP)+1]=H[1][2];} else { for(j=1;j<=size(H[1]);j++) { tj=subset(H[1][j],EE); if((subset(variables(H[1][j]),buPR)) and (tj==0)){f=f*H[1][j];} } if(leadcoef(f)!=f){BP[size(BP)+1]=f;} } } // Determining the GB basis of B in Q[u] BP=std(BP); // Determining the dimension of the anti-image // in the ring D setring D; def KK=imap(PR,BP); KK=std(KK); ideal KKM; ideal vmovM=imap(RR,vmov); // selecting the element in Q[w] for(i=1;i<=size(KK);i++) { if(subset(variables(KK[i]),vmovM)) { KKM[size(KKM)+1]=KK[i]; } } // Determining the dimensions list AIM=dimM(KKM,nv,moverdim); //int d2=AIM[1]; def AAIM=AIM[2]; setring(RR); //def BBR=imap(D,KK); def LKK=imap(D,AIM); //LKK=d2,string(BBR); return(LKK); } // Auxilliary algorithm of locus2. // The algorithm searches the basis corresponding to C, in the grobcov. // It reduces the basis modulo the component. // The result is the reduced basis BR. // for each hole of the component // it searches the segment where the hole is included // and select form its basis the polynomials // only dependent on the variables. // These polynomials are non-null in an open set of // the component, and are included in the list NoNul of non-null factors // input: C: The top of the the segment of a non-normal locus // G the grobcov of F // output: B: the anti-image of the input segment static proc AISnonor(ideal C,list G, list #) { // Initializing and defining rings int i; int j; int k; int te; int j0;int k0; int m; list DD=#; def RR=basering; def Rx=ringlist(RR); def Lx=Rx; def P=ring(Rx[1]); // ring Q[bx] Lx[1]=0; def D=ring(Lx); // ring Q[bu] def PR0=P+D; def PRx=ringlist(PR0); PRx[3][2][1]="lp"; def PR=ring(PRx); // ring Q[bx,bu] in lex order int n=size(Rx[1][2]); // number of parameters def nv=nvars(RR); // number of variables int moverdim=nv; // number of mover variables ideal vmov; // mover variables (bw) for(i=1;i<=(size(DD) div 2);i++) { if(DD[2*i-1]=="vmov"){vmov=DD[2*i];} if(DD[2*i-1]=="moverdim"){moverdim=DD[2*i];} // if(DD[2*i-1]=="version"){version=DD[2*i];} // if(DD[2*i-1]=="comment"){comment=DD[2*i];} // if(DD[2*i-1]=="grobcov"){GG=DD[2*i];} // if(DD[2*i-1]=="anti-image"){tax=DD[2*i];} } for(i=1;i<=moverdim;i++) { vmov[size(vmov)+1]=var(i+nv-moverdim); } int ddeg; int dp; list LK; ideal bu; // ideal of all variables (u) for(i=1;i<=nv;i++){bu[i]=var(i);} ideal mv; // ideal of mover varibles for(i=1;i<=nv;i++){mv[size(mv)+1]=var(i);} // Searching the basis associated to Nonor j=2; te=1; while((te) and (j<=size(G))) { k=1; while((te) and (k<=size(G[j][3]))) { if (equalideals(C,G[j][3][k][1])){j0=j; k0=k; te=0;} k++; } j++; } if(te==1){"ERROR";} def B=G[j0][2]; //"T_B=G[j0][2]="; B; //string("T_k0=",k0," G[",j0,"]="); G[j0]; return(B); } // NOT USED // Auxilliary algorithm of locus2. // The algorithm searches the basis corresponding to C, in the grobcov. // It reduces the basis modulo the component. // The result is the reduced basis BR. // for each hole of the component // it searches the segment where the hole is included // and select form its basis the polynomials // only dependent on the variables. // These polynomials are non-null in an open set of // the component, and are included in the list NoNul of non-null factors // input: F: the ideal of the locus problem // G the grobcov of F // C the top of a component of normal points // N the holes of the component // output: (d,tax,a) // where d is the dimension of the anti-image // a is the anti-image of the component and // tax is the taxonomy \"Normal\" if d is equal to the dimension of C // and \"Special\" if it is smaller. static proc NS5(ideal F,list G, ideal C, list N) { // Initializing and defining rings int i; int j; int k; int te; int j0;int k0; int m; def RR=basering; def Rx=ringlist(RR); def Lx=Rx; def P=ring(Rx[1]); Lx[1]=0; def D=ring(Lx); def RP=D+P; def PR0=P+D; def PRx=ringlist(PR0); PRx[3][2][1]="lp"; def PR=ring(PRx); int n=size(Rx[1][2]); // dimension of the tracer space def nv=nvars(RR); // number of variables ideal v; for(i=1;i<=nv;i++){v[i]=var(i);} // variables int nm; // number of mover variables ideal vmov; // mover variables if(nv>=n){nm=n;} else{nm=nv;} for(i=1;i<=nm;i++) { vmov[size(vmov)+1]=var(i+nv-nm); } // Searching the basis associated to C j=2; te=1; while((te) and (j<=size(G))) { k=1; while((te) and (k<=size(G[j][3]))) { if (equalideals(C,G[j][3][k][1])){j0=j; k0=k; te=0;} k++; } j++; } if(te==1){"ERROR";} def B=G[j0][2]; // basis associated to C //"T_B="; B; // Searching the elements in Q[vmov] on bases differents from B that are nul there // and cannot become 0 on the antiimage of C. They are placed on NoNul list NoNul; // elements in Q[vmov] on bases differents from B ideal BNoNul; ideal covertop; // basis of the segment whose top is a hole of our segment int te1; for(i=1;i<=size(N);i++) { j=2; te=1; while(te and j<=size(G)) { if(j!=j0) { k=1; while(te and k<=size(G[j][3])) { covertop=G[j][3][k][1]; if(equalideals(covertop,N[i])) { te=0; te1=1; BNoNul=G[j][2]; } else { if(redPbasis(covertop,N[i])) { te=0; te1=1; m=1; while( te1 and m<=size(G[j][3][k][2]) ) { if(equalideals(G[j][3][k][2][m] ,N[i] )==1){te1=0;} m++; } } if(te1==1){ BNoNul=G[j][2];} } k++; } if((te==0) and (te1==1)) { // Selecting the elements independent of the parameters, // They will be non null on the segment for(m=1;m<=size(BNoNul);m++) { if(indepparameterspoly(BNoNul[m])) { NoNul[size(NoNul)+1]=BNoNul[m]; } } } } j++; } } def NN=NoNul; poly kkk=1; if(size(NN)==0){NN[1]=kkk;} //"T_NN="; NN; // Union of F and B for(i=1;i<=size(F);i++) { B[size(B)+1]=F[i]; } //"T_B+F="; B; // Reducing basis B modulo C in the ring PR setring(PR); def vv=imap(RR,v); ideal BBR; def BB=imap(RR,B); def CC=imap(RR,C); attrib(CC,"IsSB",1); CC=std(CC); for(i=1;i<=size(BB);i++) { BBR[i]=reduce(BB[i],CC); } //"T_BBR="; BBR; // Selecting the elements of the ideal in Q[v] // setring RP; ideal Bv; // def Bv=imap(PR,BBR); //ideal vvv=imap(RR,v); // Bv=std(Bv); ideal BR; for(i=1;i<=size(BBR);i++) { if(subset(variables(BBR[i]),vv)){BR[size(BR)+1]=BBR[i];} } //"T_BR="; BR; // setring PR; // def BR=imap(RP,BRv); def NNN=imap(RR,NN); // Factoring and eliminating the factors in N list H; ideal BP; poly f; int tj; // ideal vv=imap(RR,v); for(i=1;i<=size(BR);i++) { f=1; H=factorize(BR[i]); if((subset(variables(H[1][2]),vv)) and (size(H[1]))==2){BP[size(BP)+1]=H[1][2];} else { for(j=1;j<=size(H[1]);j++) { tj=subset(H[1][j],NNN); if((subset(variables(H[1][j]),vv)) and (tj==0)){f=f*H[1][j];} } if(leadcoef(f)!=f){BP[size(BP)+1]=f;} } } //"T_BP="; BP; // Obtaining Am, the anti-image of C on vmov setring D; ideal A=imap(PR,BP); ideal vm=imap(RR,vmov); ideal Am; for(i=1;i<=size(A);i++) { if (subset(variables(A[i]),vm)){Am[size(Am)+1]=A[i];} } //"T_Am="; Am; list AIM=dimM(Am,nv,nm); setring(RR); def LAM=imap(D,AIM); return(LAM); } static proc dimM(ideal KKM, int nv, int nm) { def RR=basering; list L; int i; def Rx=ringlist(RR); for(i=1;i<=nm;i++) { L[i]=Rx[2][nv-nm+i]; } Rx[2]=L; intvec iv; for(i=1;i<=nm;i++){iv[i]=1;} Rx[3][1][2]=iv; def DM=ring(Rx); //"Rx="; Rx; setring(DM); ideal KKMD=imap(RR,KKM); attrib(KKMD,"IsSB",1); KKMD=std(KKMD); int d=dim(KKMD); setring(RR); def KAIM=imap(DM,KKMD); list LAIM=d,KAIM; return(LAIM); } // Procedure using only standard GB in lex(x,a) order to obtain the // component of the locus. // It is not so fine as locus and cannot evaluate the taxonomy, but // it is much simpler and efficient. // input: ideal S for determining the locus // output: the irreducible components of the locus // Data must be given in Q[a][x] proc stdlocus(ideal F) "USAGE: stdlocus(ideal F) The input ideal must be the set equations defining the locus. Calling sequence: locus(F); The input ring must be a parametrical ideal in Q[x][u], (x=tracer variables, u=remaining variables). (Inverts the concept of parameters and variables of the ring). Special routine for determining the locus of points of a geometrical construction. Given a parametric ideal F representing the system determining the locus of points (x) which verify certain properties, the call to stdlocus(F) determines the different irreducible components of the locus. This is a simple routine, using only standard Groebner basis computation, elimination and prime decomposition instead of using grobcov. It does not determine the taxonomy, nor the holes of the components RETURN:The output is a list of the tops of the components [C_1, .. , C_n] of the locus. Each component is given its top ideal p_i. NOTE: The input must be the locus system. KEYWORDS: geometrical locus; locus EXAMPLE: stdlocus; shows an example" { int i; int te; def RR=basering; list Rx=ringlist(RR); int n=npars(RR); // size(Rx[1][2]); int nv=nvars(RR); ideal vpar; ideal vvar; //"T_n="; n; //"T_nv="; nv; for(i=1;i<=n;i++){vpar[size(vpar)+1]=par(i);} for(i=1;i<=nv;i++){vvar[size(vvar)+1]=var(i);} //string("T_vpar = ", vpar," vvar = ",vvar); def P=ring(Rx[1]); Rx[1]=0; def D=ring(Rx); def RP=D+P; list Lx=ringlist(RP); setring(RP); def FF=imap(RR,F); def vvpar=imap(RR,vpar); //string("T_vvpar = ",vvpar); ideal B=std(FF); //"T_B="; B; ideal Bel; //"T_vvpar="; vvpar; for(i=1;i<=size(B);i++) { if(subset(variables(B[i]),vvpar)) {Bel[size(Bel)+1]=B[i];} } //"T_Bel="; Bel; list H; list FH; H=minAssGTZ(Bel); int t1; if(size(H)==0){t1=1;} setring RR; list empt; if(t1==1){return(empt);} else { def HH=imap(RP,H); return(HH); } } example { "EXAMPLE:"; echo = 2; if(defined(R)){kill R;} ring R=(0,x,y),(x1,y1),dp; short=0; // Concoid ideal S96=x1 ^2+y1 ^2-4,(x-2)*x1 -x*y1 +2*x,(x-x1 )^2+(y-y1 )^2-1; stdlocus(S96); } // locus(F): Special routine for determining the locus of points // of geometrical constructions. // input: The ideal of the locus equations // output: // The output components are given as // ((p1,(p11,..p1s_1),tax_1),..,(pk,(pk1,..pks_k),tax_k) // Elements 1 and 2 represent the P-canonical form of the component. // The third element tax is: // for normal point components, tax=(d,taxonomy,anti-image) // being d=dimension of the anti-image on the mover varaibles, // taxonomy='Normal' or 'Special', and // anti-image=ideal of the anti-image over the mover variables // which are taken to be the last n variables. // for non-normal point components, tax =(d,taxonomy) // being d=dimension of the component and // taxonomy='Accumulation' or 'Degenerate'. // The components are given in canonical P-representation of the subset. l // The normal locus has two kind of components: Normal and Special. // Normal component: // - each point in the component has 0-dimensional anti-image. // - the anti-image in the mover coordinates is equal to the dimension of the component // Special component: // - each point in the component has 0-dimensional anti-image. // - the anti-image in the mover coordinates is smaller than the dimension of the component // The non-normal locus has two kind of components: Accumulation and Degenerate. // Accumulation points: // - each point in the component has anti-image of dimension greater than 0. // - the component has dimension less than n-1. // Degenerate components: // - each point in the component has anti-image of dimension greater than 0. // - the component has dimension n-1. // When a normal point component has degree greater than 9, then the // taxonomy is not determined, and (n,'normal', 0) is returned as third // element of the component. (n is the dimension of the space). proc locus(ideal F, list #) "USAGE: locus(ideal F[,options]) Special routine for determining the locus of points of a geometrical construction. INPUT: The input ideal must be the set equations defining the locus. Calling sequence: locus(F[,options]); The input ring must be a parametrical ideal in Q[x1,..,xn][u1,..,um], (x=tracer variables, u=mover variables). The number of tracer variables is n=dim(space). The mover variables can be restricted by the user to the last k variaboes by adding the option \"moverdim\",k (Inverts the concept of parameters and variables of the ring). OPTIONS: An option is a pair of arguments: string, integer. To modify the default options, pairs of arguments -option name, value- of valid options must be added to the call.The algorithm allows the following options as pair of arguments: \"moverdim\", k to restrict the mover-variables of the anti-image to the last k variables. By defaulat k is equal to the number of u-variables, The restriction to a smaller value can produce an error in the character \"Normal\" or \"Special\" of a locus component. \"grobcov\", list of the previous computed grobcov(F). It is to be used when we modify externally the grobcov, for example to obtain the real grobcov. \"anti-image\",0 if the anti-image is not to be shown for \"Normal\" components. \"comments\", c: by default it is 0, but it can be set to 1. RETURN: The output is a list of the components: ((p1,(p11,..p1s_1),tax_1),..,(pk,(pk1,..pks_k),tax_k) Elements 1 and 2 of a component represent the P-canonical form of the component. The third element tax is: for normal point components, tax=(d,taxonomy,anti-image) being d=dimension of the anti-image on the mover varaibles, taxonomy=\"Normal\" or \"Special\" and anti-image=ideal of the anti-image over the mover variables. for non-normal point components, tax =(d,taxonomy) being d=dimension of the component and taxonomy=\"Accumulation\" or \"Degenerate\". The components are given in canonical P-representation. The normal locus has two kind of components: Normal and Special. Normal component: - each point in the component has 0-dimensional anti-image. - the anti-image in the mover coordinates is equal to the dimension of the component Special component: - each point in the component has 0-dimensional anti-image. - the anti-image in the mover coordinates is smaller than the dimension of the component The non-normal locus has two kind of components: Accumulation and Degenerate. Accumulation component: - each point in the component has anti-image of dimension greater than 0. - the component has dimension less than n-1. Degenerate components: - each point in the component has anti-image of dimension greater than 0. - the component has dimension n-1. When a normal point component has degree greater than 9, then the taxonomy is not determined, and (n,'normal', 0) is returned as third element of the component. (n is the dimension of the space). Given a parametric ideal F representing the system F determining the locus of points (x) which verify certain properties, the call to locus(F) determines the different classes of locus components, following the taxonomy defined in the book: A. Montes. \"The Groebner Cover\" A previous paper gives particular definitions for loci in 2d. M. Abanades, F. Botana, A. Montes, T. Recio, \"An Algebraic Taxonomy for Locus Computation in Dynamic Geometry\", Computer-Aided Design 56 (2014) 22-33. NOTE: The input must be the locus system. KEYWORDS: geometrical locus; locus; dynamic geometry EXAMPLE: locus; shows an example" { int tes=0; int i; int m; int mm; int n; def RR=basering; list GG; //Options list DD=#; ideal vmov; int nv=nvars(RR); // number of variables int moverdim=nv; //size(ringlist(RR)[1][2]); // number of parameters if(moverdim>nv){moverdim=nv;} // for(i=1;i<=nv;i++){vmov[size(vmov)+1]=var(i);} int version=2; // if(nv<4){version=2;} int comment=0; int tax=1; ideal Fm; for(i=1;i<=(size(DD) div 2);i++) { if(DD[2*i-1]=="moverdim"){moverdim=DD[2*i];} if(DD[2*i-1]=="vmov"){vmov=DD[2*i];} if(DD[2*i-1]=="version"){version=DD[2*i];} if(DD[2*i-1]=="comment"){comment=DD[2*i];} if(DD[2*i-1]=="grobcov"){GG=DD[2*i];} if(DD[2*i-1]=="anti-image"){tax=DD[2*i];} } for(i=1;i<=moverdim;i++){vmov[size(vmov)+1]=var(i+nv-moverdim);} if(size(GG)==0){GG=grobcov(F);} //"T_GG=grobcov(F)="; GG; int j; int k; int te; def B0=GG[1][2]; def H0=GG[1][3][1][1]; list nGP; if (equalideals(B0,ideal(1)) ) { if(version==2){return(locus2(GG,F,DD));} else{return(locus0(GG,DD));} } else { n=nvars(R); ideal vB; ideal N; for(i=1;i<=size(B0);i++) { if(subset(variables(B0[i]),vmov)){N[size(N)+1]=B0[i];} } attrib(N,"IsSB",1); N=std(N); if((size(N))>=2) { //def dN=dim(N); te=indepparameters(N); if(te) { string("locus detected that the mover must avoid points (",N,") in order to obtain the correct locus");" "; //eliminates segments of GG where N is contained in the basis nGP[1]=GG[1]; nGP[1][1]=ideal(1); nGP[1][2]=ideal(1); def GP=GG; ideal BP; ideal fBP; for(j=2;j<=size(GP);j++) { te=1; k=1; BP=GP[j][2]; // eliminating multiple factors in the polynomials of BP for(mm=1;mm<=size(BP);mm++) { fBP=factorize(BP[mm],1); BP[mm]=1; for(m=1;m<=size(fBP);m++) { BP[mm]=BP[mm]*fBP[m]; } } // end eliminating multiple factors while((te==1) and (k<=size(N))) { if(pdivi(N[k],BP)[1]!=0){te=0;} k++; } if(te==0){nGP[size(nGP)+1]=GP[j];} } } } else { nGP=GG; " ";string("Unavoidable ",moverdim,"-dimensional locus"); //string("Try option 'vmov',ideal(of mover variables) to avoid some point of the mover"); //" ";"Elements of the basis of the generic segment in mover variables="; N;" "; list L; return(L); } } if(comment>0){"Input for locus2 GB="; nGP; "input for locus F="; F;} if(version==2) { //"T_nGP="; nGP; def LL=locus2(nGP,F,DD); } else{def LL=locus0(nGP,DD);} // kill @RP; kill @P; kill @R; return(LL); } example { "EXAMPLE:"; echo = 2; if(defined(R)){kill R;} ring R=(0,x,y),(x1,y1),dp; short=0; // Concoid ideal S96=x1 ^2+y1 ^2-4,(x-2)*x1 -x*y1 +2*x,(x-x1 )^2+(y-y1 )^2-1; locus(S96); // EXAMPLE // Consider two parallel planes z1=-1 and z1=1, and two orthogonal parabolas on them. // Determine the locus generated by the lines that rely the two parabolas // through the points having parallel tangent vectors. if(defined(R)){kill R;} ring R=(0,x,y,z),(lam,x2,y2,z2,x1,y1,z1), lp; short=0; ideal L=z1+1, x1^2-y1, z2-1, y2^2-x2, 4*x1*y2-1, x-x1-lam*(x2-x1), y-y1-lam*(y2-y1), z-z1-lam*(z2-z1); locus(L); // uses "moverdim",7 // Now observe the effect on the antiimage of option moverdim: // It does not take into account that lam is free. locus(L,"moverdim",3); // Observe that with the option "mpverdim",3 the taxonomy becomes Special because considering only x1,y1,z1 as mover variables does not take into account that lam is free in the global anti-image. } // locusdg(G): Special routine for determining the locus of points // of geometrical constructions in Dynamic Geometry. // It is to be applied to the output of locus and selects // as 'Relevant' the 'Normal' and the 'Accumulation' // components. // input: The output of locus(S); // output: // list, the canonical P-representation of the 'Relevant' components of the locus. // The output components are given as // ((p1,(p11,..p1s_1),type_1,level_1),..,(pk,(pk1,..pks_k),type_k,level_k) // The components are given in canonical P-representation of the subset. // If all levels of a class of locus are 1, then the set is locally closed. Otherwise the level // gives the depth of the component of the constructible set. proc locusdg(list L) "USAGE: locusdg(list L) Calling sequence: locusdg(locus(S)). RETURN: The output is the list of the \"Relevant\" components of the locus in Dynamic Geometry [C1,..,C:m], where C_i= [p_i,[p_i1,..p_is_i], \"Relevant\", level_i] The \"Relevant\" components are \"Normal\" and \"Accumulation\" components of the locus. (See help for locus). KEYWORDS: geometrical locus; locus; dynamic geometry EXAMPLE: locusdg; shows an example" { list LL; int i; for(i=1;i<=size(L);i++) { if(typeof(L[i][3][2])=="string") { if((L[i][3][2]=="Normal") or (L[i][3][2]=="Accumulation")){L[i][3][2]="Relevant"; LL[size(LL)+1]=L[i];} } } return(LL); } example { "EXAMPLE:"; echo = 2; if(defined(R)){kill R;}; ring R=(0,a,b),(x,y),dp; short=0; // Concoid ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1; def L96=locus(S96); L96; locusdg(L96); } // locusto: Transforms the output of locus, locusdg, envelop // into a string that can be reed from different computational systems. // input: // list L: The output of locus or locusdg or envelop. // output: // string s: Converts the input into a string readable by other programs proc locusto(list L) "USAGE: locusto(list L); The argument must be the output of locus or locusdg or envelop. It transforms the output into a string in standard form readable in other languages, not only Singular (Geogebra). RETURN: The locus in string standard form NOTE: It can only be called after computing either - locus(F) -> locusto( locus(F) ) - locusdg(locus(F)) -> locusto( locusdg(locus(F)) ) - envelop(F,C) -> locusto( envelop(F,C) ) KEYWORDS: geometrical locus; locus; envelop EXAMPLE: locusto; shows an example" { int i; int j; int k; string s="["; string sf="]"; string st=s+sf; if(size(L)==0){return(st);} ideal p; ideal q; for(i=1;i<=size(L);i++) { s=string(s,"[["); for (j=1;j<=size(L[i][1]);j++) { s=string(s,L[i][1][j],","); } s[size(s)]="]"; s=string(s,",["); for(j=1;j<=size(L[i][2]);j++) { s=string(s,"["); for(k=1;k<=size(L[i][2][j]);k++) { s=string(s,L[i][2][j][k],","); } s[size(s)]="]"; s=string(s,","); } s[size(s)]="]"; s=string(s,"]"); if(size(L[i])>=3) { s=string(s,",["); if(typeof(L[i][3])=="string") { s=string(s,string(L[i][3]),"]]"); } else { for(k=1;k<=size(L[i][3]);k++) { s=string(s,"[",L[i][3][k],"],"); } s[size(s)]="]"; s=string(s,"]"); } } if(size(L[i])>=4) { s[size(s)]=","; s=string(s,string(L[i][4]),"],"); } s[size(s)]="]"; s=string(s,","); } s[size(s)]="]"; return(s); } example { "EXAMPLE:"; echo = 2; if(defined(R)){kill R;} ring R=(0,x,y),(x1,y1),dp; short=0; ideal S=x1^2+y1^2-4,(y-2)*x1-x*y1+2*x,(x-x1)^2+(y-y1)^2-1; def L=locus(S); locusto(L); locusto(locusdg(L)); } // envelop // Input: // poly F: the polynomial defining the family of hypersurfaces in ring R=0,(x_1,..,x_n),(u_1,..,u_m),lp; // ideal C=g1,..,g_{n-1}: the set of constraints; // options. // Output: the components of the envolvent; proc envelop(poly F, ideal C, list #) "USAGE: envelop(poly F,ideal C[,options]); poly F must represent the family of hyper-surfaces for which on want to compute its envelop. ideal C must be the ideal of restrictions on the variables defining the family, and should contain less polynomials than the number of variables. (x_1,..,x_n) are the variables of the hyper-surfaces of F, that are considered as parameters of the parametric ring. (u_1,..,u_m) are the parameteres of the hyper-surfaces, that are considered as variables of the parametric ring. Calling sequence: ring R=(0,x_1,..,x_n),(u_1,..,u_m),lp; poly F=F(x_1,..,x_n,u_1,..,u_m); ideal C=g_1(u_1,..u_m),..,g_s(u_1,..u_m); envelop(F,C[,options]); where snv){moverdim=nv;} for(i=1;i<=moverdim;i++){vmov[size(vmov)+1]=var(i+nv-moverdim);} //for(i=1;i<=nv;i++) {vmov[size(vmov)+1]=var(i);} int version=2; int comment=0; ideal Fm; int tax=1; int familyinfo; for(i=1;i<=(size(DD) div 2);i++) { if(DD[2*i-1]=="anti-image"){tax=DD[2*i];} if(DD[2*i-1]=="vmov"){vmov=DD[2*i];} if(DD[2*i-1]=="version"){version=2;DD[2*i]=2;} if(DD[2*i-1]=="comment"){comment=DD[2*i];} if(DD[2*i-1]=="familyinfo"){familyinfo=0;DD[2*i]=0;} }; int ng=size(C); ideal S=F; //"T_C="; C; for(i=1;i<=size(C);i++){S[size(S)+1]=C[i];} //"T_S="; S; int s=nv-ng; if(s>0) { matrix M[ng+1][ng+1]; def cc=comb(nv,ng+1); //string("nv=",nv," ng=",ng," s=",s); //"T_cc="; cc; poly J; for(k=1;k<=size(cc);k++) { for(j=1;j<=ng+1;j++) { M[1,j]=diff(F,var(cc[k][j])); } for(i=1;i<=ng;i++) { for(j=1;j<=ng+1;j++) { M[i+1,j]=diff(C[i],var(cc[k][j])); } } J=det(M); S[size(S)+1]=J; } } if(comment>0){"System S before grobcov ="; S;} //def G=grobcov(S,DD); list HHH; DD[size(DD)+1]="moverdim"; // new option DD[size(DD)+1]=moverdim; // string("moverdim=",moverdim); def L=locus(S,DD); // "L="; L; list GL; ideal fam; ideal env; def P=ring(Rx[1]); Rx[1]=0; def D=ring(Rx); def RP=P+D; list LL; list NormalComp; ideal Gi; ideal BBBB; poly B0; return(L); } example { "EXAMPLE:"; echo = 2; if(defined(R)){kill R;} ring R=(0,x,y),(r,s,y1,x1),lp; poly F=(x-x1)^2+(y-y1)^2-r; ideal g=(x1-2*(s+r))^2+(y1-s)^2-s; def E=envelop(F,g); E; E[1][1]; string("Taxonomy=", E[1][3][2]); // EXAMPLE // Steiner Deltoid // 1. Consider the circle x1^2+y1^2-1=0, and a mover point M(x1,y1) on it. // 2. Consider the triangle A(0,1), B(-1,0), C(1,0). // 3. Consider lines passing through M perpendicular to two sides of ABC triangle. // 4. Obtain the envelop of the lines above. if(defined(R)){kill R;} ring R=(0,x,y),(x1,y1,x2,y2),lp; short=0; ideal C=(x1)^2+(y1)^2-1, x2+y2-1, x2-y2-x1+y1; matrix M[3][3]=x,y,1,x2,y2,1,x1,0,1; poly F=det(M); // Curves Family F F; // Conditions C= C; envelop(F,C); } proc AssocTanToEnv(poly F,ideal C, ideal E,list #) "USAGE: AssocTanToEnv(poly F,ideal C,ideal E); poly F must be the family of hyper-surfaces whose envelope is analyzed. It must be defined in the ring R=Q[x_1.,,x_n][u_1,..,u_m], ideal C must be the ideal of restrictions in the variables u1,..um for defining the family. C must contain less polynomials than m. ideal E must be a component of envelop(F,C), previously computed. (x_1,..,x_n) are the variables of the hypersurfaces of F, that are considered as parameters of the parametric ring. (u_1,..,u_m) are the parameteres of the hyper-surfaces, that are considered as variables of the parametric ring. Having computed an envelop component E of a family of hyper-surfaces F, with constraints C, it returns the parameter values of the associated tangent hyper-surface of the family passing at one point of the envelop component E. Calling sequence: (s0) { matrix M[ng+1][ng+1]; def cc=comb(nv,ng+1); poly J; for(k=1;k<=size(cc);k++) { for(j=1;j<=ng+1;j++) { M[1,j]=diff(F,var(cc[k][j])); } for(i=1;i<=ng;i++) { for(j=1;j<=ng+1;j++) { M[i+1,j]=diff(C[i],var(cc[k][j])); } } J=det(M); S[size(S)+1]=J; } } for(i=1;i<=size(EE);i++) { S[size(S)+1]=EE[i]; } //if(comment>0){"System S before grobcov ="; S;} def G=grobcov(S,DD); //"T_G=";G; list GG; for(i=2;i<=size(G);i++) { GG[size(GG)+1]=G[i]; } G=GG; //"T_G=";G; if(moreinfo>0){return(G);} else { int t=0; list HH; i=1; while(t==0 and i<=size(G)) { //string("T_G[",i,"][3][1][1][1]="); G[i][3][1][1][1]; //string("T_EE="); EE; if(G[i][3][1][1][1]==EE) { t=1; HH=G[i]; } i++; } return(HH); } return(G); } example { "EXAMPLE:"; echo = 2; if(defined(R)){kill R;} ring R=(0,x,y),(r,s,y1,x1),lp; poly F=(x-x1)^2+(y-y1)^2-r; ideal g=(x1-2*(s+r))^2+(y1-s)^2-s; def E=envelop(F,g); E; def A=AssocTanToEnv(F,g,E[1][1][1]); A; def M1=coef(A[2][1],x1); def M2=coef(A[2][2],y1); def M3=coef(A[2][3],s); def M4=coef(A[2][4],r); "x1=";-M1[2,2]/M1[2,1]; "y1=";-M2[2,2]/M2[2,1]; "s=";-M3[2,2]/M3[2,1]; "r=";-M4[2,2]/M4[2,1]; } proc FamElemsAtEnvCompPoints(poly F,ideal C, ideal E) "USAGE: FamElemsAtEnvCompPoints(poly F,ideal C,ideal E); poly F must be the family of hyper-surfaces whose envelope is analyzed. It must be defined in the ring R=Q[x_1.,,x_n][u_1,..,u_m], ideal C must be the ideal of restrictions on the variables u1,..um. Must contain less polynomials than m. ideal E must be a component of envelop(F,C), previously computed. After computing the envelop of a family of hyper-surfaces F, with constraints C, Consider a component with top E. The call to FamElemsAtEnvCompPoints(F,C,E) returns the parameter values of the set of all hyper-surfaces of the family passing at one point of the envelop component E. Calling sequence: ring R=(0,x_1,..,x_n),(u_1,..,u_m),lp; poly F=F(x_1,..,x_n,u_1,..,u_m); ideal C=g_1(u_1,..u_m),..,g_s(u_1,..u_m); poly E(x_1,..,x_n); FamElemsAtEnvCompPoints(F,C,E[,options]); RETURN: list [lpp,basis,segment]. The basis determines the parameter values of the of hyper-surfaces that pass at a fixed point of the envelop component E. The lpp determines the dimension of the set. The segment is the component and is given in Prep. Fixing the values of (x_1,..,x_n) inside E, the basis allows to detemine the values of the parameters (u_1,..u_m), of the hyper-surfaces passing at a point of E. See the book A. Montes. \"The Groebner Cover\" (Discussing Parametric Polynomial Systems). NOTE: grobcov is called internally. The basering R, must be of the form Q[a][x] (a=parameters, x=variables). KEYWORDS: geometrical locus; locus; envelop; associated tangent EXAMPLE: FamElemsAtEnvCompPoints; shows an example" { int i; int moreinfo=0; ideal S=C; ideal EE=E; S[size(S)+1]=F; //S[size(S)+1]=E; for(i=1;i<=size(E);i++){S[size(S)+1]=E[i];} def G=grobcov(S); list GG; for(i=2; i<=size(G); i++) { GG[size(GG)+1]=G[i]; } if(moreinfo>0){return(GG);} else { int t=0; list HH; i=1; while(t==0 and i<=size(G)) { //string("T_G[",i,"][3][1][1][1]="); G[i][3][1][1][1]; //string("T_EE="); EE; if(G[i][3][1][1][1]==EE) { t=1; HH=G[i]; } i++; } return(HH); } } example { "EXAMPLE:"; echo = 2; if(defined(R)){kill R;} ring R=(0,x,y),(t),dp; short=0; poly F=(x-5*t)^2+y^2-9*t^2; ideal C; def Env=envelop(F,C); Env; // E is a component of the envelope: ideal E=Env[1][1]; E; def A=AssocTanToEnv(F,C,E); A; // The basis of the parameter values of the associated // tangent component is A[2][1]; // Thus t=-(5/12)*y, and the associated tangent family // element at (x,y) is subst(F,t,-(5/12)*y); def FE=FamElemsAtEnvCompPoints(F,C,E); FE; factorize(FE[2][1]); // Thus the unique family element passing through the envelope point (x,y) // corresponds to the value of t of the Associated Tangent // EXAMPLE: if(defined(R)){kill R;} ring R=(0,x,y),(r,s,y1,x1),lp; poly F=(x-x1)^2+(y-y1)^2-r; ideal g=(x1-2*(s+r))^2+(y1-s)^2-s; def E=envelop(F,g); E; def A=AssocTanToEnv(F,g,E[1][1][1]); A; def M1=coef(A[2][1],x1); def M2=coef(A[2][2],y1); def M3=coef(A[2][3],s); def M4=coef(A[2][4],r); // The parameter values corresponding to the family // element tangent at point (x,y) of the envelope are: "x1=";-M1[2,2]/M1[2,1]; "y1=";-M2[2,2]/M2[2,1]; "s=";-M3[2,2]/M3[2,1]; "r=";-M4[2,2]/M4[2,1]; // Now detect if there are other family elements passing at this point: def FE=FamElemsAtEnvCompPoints(F,g,E[1][1][1]); FE; // FE[1] is the set of lpp. It has dimension 4-2=2. // Thus there are points of the envelope at which // they pass infinitely many circles of the family. // To separe the points of the envelope further analysis must be done. } // discrim proc discrim(poly F0, poly x0) "USAGE: discrim(f,x); poly f: the polynomial in Q[a][x] or Q[x] of degree 2 in x poly x: can be a variable or a parameter of the ring. RETURN: the factorized discriminant of f wrt x for discussing its sign KEYWORDS: second degree; solve EXAMPLE: discrim; shows an example" { def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); Rx[1]=0; def D=ring(Rx); def RP=D+P; int i; int te; int d; int dd; if(size(ringlist(RR)[1])>0) { te=1; // setglobalrings(); setring RP; poly F=imap(RR,F0); poly X=imap(RR,x0); } else {poly F=F0; poly X=x0;} matrix M=coef(F,X); d=deg(M[1,1]); if(d>2){"Degree is higher than 2. No discriminant"; setring RR; return();} poly dis=(M[2,2])^2-4*M[2,1]*M[2,3]; def disp=factorize(dis,0); if(te==0){return(disp);} else { setring RR; def disp0=imap(RP,disp); return(disp0); } } example { "EXAMPLE:"; echo = 2; if(defined(R)){kill R;} ring R=(0,a,b,c),(x,y),dp; short=0; poly f=a*x^2*y+b*x*y+c*y; discrim(f,x); } // AddLocus: auxilliary routine for locus0 that computes the components of the constructible: // Input: the list of locally closed sets to be added, each with its type as third argument // L=[ [LC[11],..,LC[1k_1],.., [LC[r1],..,LC[rk_r] ] where // LC[1]=[p1,[p11,..,p1k],typ] // Output: the list of components of the constructible union of L, with the type of the corresponding top // and the level of the constructible // L4= [[v1,p1,[p11,..,p1l],typ_1,level]_1 ,.. [vs,ps,[ps1,..,psl],typ_s,level_s] static proc AddLocus(list L) { list L1; int i; int j; list L2; list L3; list l1; list l2; intvec v; for(i=1; i<=size(L); i++) { for(j=1;j<=size(L[i]);j++) { l1[1]=L[i][j][1]; l1[2]=L[i][j][2]; l2[1]=l1[1]; if(size(L[i][j])>2){l2[3]=L[i][j][3];} v[1]=i; v[2]=j; l2[2]=v; L1[size(L1)+1]=l1; L2[size(L2)+1]=l2; } } L3=LocusConsLevels(L1); list L4; int level; ideal p1; ideal pp1; int t; int k; int k0; string typ; list l4; for(i=1;i<=size(L3);i++) { level=L3[i][1]; for(j=1;j<=size(L3[i][2]);j++) { p1=L3[i][2][j][1]; t=1; k=1; while((t==1) and (k<=size(L2))) { pp1=L2[k][1]; if(equalideals(p1,pp1)){t=0; k0=k;} k++; } if(t==0) { v=L2[k0][2]; l4[1]=v; l4[2]=p1; l4[3]=L3[i][2][j][2]; l4[5]=level; if(size(L2[k0])>2){l4[4]=L2[k0][3];} L4[size(L4)+1]=l4; } else{"ERROR p1 NOT FOUND";} } } return(L4); } // Input L: list of components in P-rep to be added // [ [[p_1,[p_11,..,p_1,r1]],..[p_k,[p_k1,..,p_kr_k]] ] // Output: // list of lists of levels of the different locally closed sets of // the canonical P-rep of the constructible. // [ [level_1,[ [Comp_11,..Comp_1r_1] ] ], .. , // [level_s,[ [Comp_s1,..Comp_sr_1] ] // ] // where level_i=i, Comp_ij=[ p_i,[p_i1,..,p_it_i] ] is a prime component. // LocusConsLevels: given a set of components of locally closed sets in P-representation, it builds the // canonical P-representation of the corresponding constructible set of its union, // including levels it they are. static proc LocusConsLevels(list L) { list Lc; list Sc; int i; for(i=1;i<=size(L);i++) { Sc=PtoCrep0(list(L[i])); Lc[size(Lc)+1]=Sc; } list S=ConsLevels(Lc); S=Levels(S); list Sout; list Lev; for(i=1;i<=size(S);i++) { Lev=list(S[i][1],Prep(S[i][2][1],S[i][2][2])); Sout[size(Sout)+1]=Lev; } return(Sout); } // used in NS // returns 0 if E does not reduce modulo N // returns 1 if it reduces static proc redPbasis(ideal E, ideal N) { int i; def RR=basering; def Rx=ringlist(RR); def Lx=Rx; def P=ring(Rx[1]); setring P; def EP=imap(RR,E); def NP=imap(RR,N); NP=std(NP); list L; int red=1; i=1; while(red and (i<=size(EP))) { if(reduce(EP[i],NP,5)!=0){red=0;} i++; } setring RR; return(red); } //******************** End locus and envelop ****************************** //********************* Begin WLemma ********************** // input ideal F in @R // ideal a in @R but only depending on parameters // F is a generating ideal in V(a); // output: ideal b in @R but depending only on parameters // ideal G=GBasis(F) in V(a) \ V(b) proc WLemma(ideal F,ideal a, list #) "USAGE: WLemma(F,A[,options]); The first argument ideal F in Q[x_1,..,x_n][u_1,..,u_m]; The second argument ideal A in Q[x_1,..,x_n]. Calling sequence: ring R=(0,x_1,..,x_n),(u_1,..,u_m),lp; ideal F=f_1(x_1,..,x_n,u_1,..,u_m),.., f_s(x_1,..,x_n,u_1,..,u_m); ideal A=g_1(u_1,..u_m),..,g_s(u_1,..u_m); list # : Options Calling sequence: WLemma(F,A[,options]); Given the ideal F and ideal A it returns the list (lpp,B,S) were B is the reduced Groebner basis of the specialized F over the segment S, subset of V(A) with top A, determined by Wibmer's Lemma. S is determined in P-representation (or optionally in C-representation). The basis is given by I-regular functions. OPTIONS: either (\"rep\", 0) or (\"rep\",1) the representation of the resulting segment, by default is 0 =P-representation, (default) but can be set to 1=C-representation. RETURN: list of [lpp,B,S] = [leading power product, basis,segment], being B the reduced Groebner Basis given by I-regular functions in full representation, of the specialized ideal F on the segment S, subset of V(A) with top A. given in P- or C-representation. It is the result of Wibmer's Lemma. See A. Montes , M. Wibmer, \"Groebner Bases for Polynomial Systems with parameters\". JSC 45 (2010) 1391-1425.) or the book A. Montes. \"The Groebner Cover\" (Discussing Parametric Polynomial Systems). NOTE: The basering R, must be of the form Q[a][x] (a=parameters, x=variables). KEYWORDS: Wibmer's Lemma EXAMPLE: WLemma; shows an example" { list L=#; int rep=0; int i; int j; if(size(L)>0) { for(i=1;i<=size(L) div 2;i++) { if(L[2*i-1]=="rep"){rep=L[2*i];} } } def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); Rx[1]=0; def D=ring(Rx); def RP=D+P; setring(RP); ideal FF=imap(RR,F); FF=std(FF); ideal AA=imap(RR,a); AA=std(AA); FF=FF,AA; FF=std(FF); ideal FFa; poly r; for(i=1; i<=size(FF);i++) { r=reduce(FF[i],AA); if(r!=0){FFa[size(FFa)+1]=r;} } // FFa is GB(F+a,>xa) setring RR; ideal Fa=imap(RP,FFa); ideal AAA=imap(RP,AA); ideal lppFa; ideal lcFa; for(i=1;i<=size(Fa);i++) { lppFa[size(lppFa)+1]=leadmonom(Fa[i]); lcFa[size(lcFa)+1]=leadcoef(Fa[i]); } // "T_lppFa="; lppFa; // "T_lcFa="; lcFa; setring RP; ideal lccr=imap(RR,lppFa); lccr=std(lccr); setring RR; ideal lcc=imap(RP,lccr); list J; list Jx; ideal Jci; ideal Jxi; list B; // "T_lcc="; lcc; for(i=1;i<=size(lcc);i++) { kill Jci; ideal Jci; kill Jxi; ideal Jxi; for(j=1;j<=size(Fa);j++) { if(lppFa[j]==lcc[i]) { Jci[size(Jci)+1]=lcFa[j]; Jxi[size(Jxi)+1]=Fa[j]; } } J[size(J)+1]=Jci; B[size(B)+1]=Jxi; } // "T_J="; J; if(size(J)>0) { setring P; list Jp=imap(RR,J); ideal JL=product(Jp); // JL=prod(lc(Fa)) def AAAA=imap(RR,AAA); // "T_AAA="; AAA; // "T_JLA="; JLA; def CPR=Crep(AAAA, JL); def PPR=Prep(AAAA,JL); } setring RR; if(size(J)>0) { def JLA=imap(P,JL); def PR=imap(P,PPR); def CR=imap(P,CPR); // PR=Prep(a,b) // CR=Crep(a,b) for(i=1;i<=size(B);i++) { for(j=1;j<=size(B[i]);j++) { B[i][j]=pnormalf(B[i][j],CR[1],CR[2]); } B[i]=elimrepeated(B[i]); } // B=reduced basis on CR //"T_PR="; PR; //"T_CR="; CR; //"T_B="; B; if(rep==1){return(list(lcc,B,CR));} else{return(list(lcc,B,PR));} } else { "PIP"; lcc=ideal(0); B=ideal(0); list NN; NN[1]=list(AAA,ideal(1)); return(list(lcc,B,NN)); } } example { "EXAMPLE:"; echo = 2; if(defined(RE)){kill RE;} ring RE=(0,a,b,c,d,e,f),(x,y),lp; ideal F=a*x^2+b*x*y+c*y^2,d*x^2+e*x*y+f*y^2; ideal A=a*e-b*d; WLemma(F,A); WLemma(F,A,"rep",1); } // Detect if ideal J is in the list of ideals L // Input ideal J, list L // Output: 1 if J is in L, and 0 if not static proc idinlist(ideal J,list L) { int i=0; int te=0; while(te==0 and i<=size(L)-1) { i++; if(equalideals(J,L[i])){te=1;} } return(te); } // input ideal F in Kˆa][x] // output: a disjoint CGS in full representation of the ideal F using Wibmer's Lemma WLemma proc WLcgs(ideal F) USAGE: WLcgs(ideal F) // WLemma(F,A[,options]); ideal F in Q[x_1,..,x_n][u_1,..,u_m]; ring R=(0,x_1,..,x_n),(u_1,..,u_m),lp; ideal F=f_1(x_1,..,x_n,u_1,..,u_m),.., f_s(x_1,..,x_n,u_1,..,u_m); list # : Options Calling sequence: WLcgs(ideal F) Given the ideal F it returns the list of (lpp,B,S)[i] of the grobcov. B[i] is the reduced Groebner basis in full-representation of the specialized F over the segments S[i], S is determined in P-representation (or optionally in C-representation). The basis is given by I-regular functions in full-representation OPTIONS: either (\"rep\", 0) or (\"rep\",1) the representation of the resulting segment, by default is 0 =P-representation, (default) but can be set to 1=C-representation. RETURN: list of [lpp,B,S][i] = [leading power product, basis,segment], being B[i] the reduced Groebner Basis given by I-regular functions in full representation, of the specialized ideal F on the segment S[i], given in P-representation. It is the result of Wibmer's Lemma. See A. Montes , M. Wibmer, \"Groebner Bases for Polynomial Systems with parameters\". JSC 45 (2010) 1391-1425.) or the book A. Montes. \"The Groebner Cover\" (Discussing Parametric Polynomial Systems). NOTE: The basering R, must be of the form Q[a][x] (a=parameters, x=variables). KEYWORDS: Wibmer's Lemma EXAMPLE: WLcgs; shows an example" { int i,j; list Etot; Etot[1]=ideal(0); list Epend=Etot; list G; list G0; list N; ideal a; while (size(Epend)>0) { a=Epend[1]; Epend=elimidealfromlist(Epend,a); G0=WLemma(F,a); if(size(G0)>0) { G[size(G)+1]=G0; //"T_G0="; G0; N=G0[3][1][2]; //"T_N="; N; for(i=1;i<=size(N);i++) { if(not(equalideals(N[i],ideal(1)) or idinlist(N[i],Etot))) { Etot[size(Etot)+1]=N[i]; Epend[size(Epend)+1]=N[i]; } //"T_i="; i; } //"T_Etot="; Etot; //"T_Epend=";Epend; } } return(G); } example { "EXAMPLE:"; echo = 2; if(defined(RRR)){kill RRR;} ring RRR=(0,b,c,d,e,f),(x,y,t),lp; short=0; ideal S=x^2+2*c*x*y+2*d*x*t+b*y^2+2*e*y*t+f*t^2, x+c*y+d*t,c*x+b*y+e*t; grobcov(S); WLcgs(S); } //********************* End WLemma ************************ // Not used static proc redbasis(ideal B, ideal C) { int i; def RR=basering; def Rx=ringlist(RR); def Lx=Rx; def P=ring(Rx[1]); Lx[1]=0; def D=ring(Lx); def RP=D+P; setring RP; ideal BB=imap(RR,B); ideal CC=imap(RR,C); attrib(CC,"IsSB",1); CC=std(CC); for(i=1;i<=size(BB);i++) { BB[i]=reduce(BB[i],CC); } setring(RR); def BBB=imap(RP,BB); return(BBB); } // not used // Input: ideals E, N // Output: the ideal N without the polynomials in E // Works in any kind of ideal static proc idminusid(ideal E,ideal N) { int i; int j; ideal h; int te=1; for(i=1; i<=size(N);i++) { te=1; for(j=1;j<=size(E);j++) { if(N[i]==E[j]){te=0;} } if(te==1){h[size(h)+1]=N[i];} } return(h); } // not used // eliminar els factors de cada polinomi de F que estiguin a N\ E static proc simpB(ideal F,ideal E,ideal N) { ideal FF; poly ff; int i; int j; ideal J=idminusid(E,N); //"T_J="; J; for(i=1;i<=size(F);i++) { for(j=1;j<=size(J);j++) { ff=elimfacsinP(F[i],J[j]); } FF[size(FF)+1]=ff; } return(FF); } // used in simpB that is not used static proc elimfacsinP(poly f,poly g) { def RR=basering; def Rx=ringlist(RR); int i; int j; int n=size(Rx[1][2]); def Lx=Rx; Lx[1]=0; def D=ring(Lx); def P=ring(Rx[1]); def RP=D+P; setring P; ideal vp; for(i=1;i<=n;i++) { vp[size(vp)+1]=var(i); } setring RP; def gg=imap(RR,g); ideal vpr=imap(P,vp); poly ff=imap(RR,f); def L=factorize(ff); def L1=L[1]; poly p=1; for(i=1;i<=size(L1);i++) { if(L1[i]==gg){;} else{p=p*L1[i];} } setring RR; def pp=imap(RP,p); return(pp); } //****************************** Begin ADGT ************************* // used in ADGT // Given G=grobcov(F,"rep",1) to have the GC in C-representation, // Grob1Levels determines the canonical levels of the constructible subset // of the parameter space for which there exist solutions of F // To be called in Q[a][x] proc Grob1Levels(list G) "USAGE: Grob1Levels(list G); G is the output of grobcov(F,\"rep\",1) for obtaining the segments in C-rep. Then Grob!Levels, selects the set of segments S of G having solutions (i.e. with basis different from 1), and determines the canonical levels of this constructible set. To be called in a ring Q[a][x]. RETURN: The list of ideals [a1,a2,...,at] representing the closures of the canonical levels of S and its complement C wrt to the closure of S. The levels of S and C are Levels of S: [a1,a2],[a3,a4],... Levels of C: [a2,a3],[a4,a5],... S=V(a1) \ V(a2) u V(a3) \ V(a4) u ... C=V(a2 \ V(a3) u V(a4) \ V(a5) u ... The expression of S can be obtained from the output of Grob1Levels by the call to Levels. NOTE: The algorithm was described in J.M. Brunat, A. Montes. \"Computing the canonical representation of constructible sets.\" Math. Comput. Sci. (2016) 19: 165-178. KEYWORDS: constructible set; locally closed set; canonical form EXAMPLE: Grob1Levels; shows an example" { int i; list S; def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); setring(RR); for(i=1;i<=size(G);i++) // select the segments with solutions { if(not(G[i][1][1][1]==1)) { S[size(S)+1]=G[i][3]; } } if(size(S)==0) { list L; L[1]=1; } else { setring P; def SP=imap(RR,S); list LP=ConsLevels(SP); // construct the levels of the constructible setring RR; def L=imap(P,LP); } return(L); } example { "EXAMPLE:"; echo = 2; if (defined(R)) {kill R;} ring R=(0,x,y),(x1,y1,x2,y2),lp; ideal F=-y*x1+(x-1)*y1+y, (x-1)*(x1+1)+y*y1, -y*x2+(x+1)*y2-y, (x+1)*(x2-1)+y*y2, (x1-x)^2+y1^2-(x1-x)^2-y2^2; def G=grobcov(F,"rep",1); G; def L=Grob1Levels(G); L; def LL=Levels(L); LL; } // Auxiliary rutine for intersecting ideal only in the parameters a // when the ideals are defined in Q[a][x] proc intersectpar(list S) "USAGE: interectpar(list of ideals S) list S=ideal I1,...,ideal Ik; RETURN: The intersection of the ideals I1 ... Ik in Q[x,a] NOTE: The routine is called in Q[a][x], The ideals I1,..,Ik can be ideals depending only on [a] or on [x,a] EXAMPLE: intersectpar shows an example" { def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); Rx[1]=0; def D=ring(Rx); def RP=D+P; setring(RP); def SP=imap(RR,S); //"T_SP="; SP; //"T_typeof(SP[1])="; typeof(SP[1]); ideal EP; EP=SP[1]; int i; for(i=2;i<=size(SP);i++) { EP=intersect(EP,SP[i]); } //def EP=intersect(SPL); setring RR; def E=imap(RP,EP); return(E); } example { "EXAMPLE:"; echo = 2; if(defined(R)){kill R;} ring R=(0,x,y,z),(x1,y1,z1),lp; ideal I1=x+y*z*x1; ideal I2=x-y*z*y1; ideal I3=x+y+z*z1; list S=I1,I2,I3; S; intersectpar(S); } proc ADGT(ideal H,ideal T,ideal H1,ideal T1,list #) "USAGE: ADGT(ideal H, ideal T, ideal H1,ideal T1[,options]); H: ideal in Q[a][x] hypothesis T: ideal in Q[a][x] thesis H1: ideal in Q[a][x] negative hypothesis, only dependent on [a] T1: ideal in Q[a][x] negative thesis of the Proposition (H and not H1) => (T and not T1) RETURN: The list [[1,S1],[2,S2],..], S1, S2, .. represent the set of parameter values that must be verified as supplementary conditions for the Proposition to become a Theorem. They are given by default in Prep. If the proposition is generally true, (the proposition is already a theorem), then the generic segment of the internal grobcov called is also returned to provide information about the values of the variables determined for every value of the parameters. If the propsition is false for every values of the parameters, then the emply list is returned. OPTIONS: An option is a pair of arguments: string, integer. To modify the default options, pairs of arguments -option name, value- of valid options must be added to the call. Option \"rep\",0-1: The default is (\"rep\",0) and then the segments are given in canonical P-representation. Option (\"rep\",1) represents the segments in canonical C-representation, Option \"gseg\",0-1: The default is \"gseg\",1 and then when the proposition is generally true, ADGT outputs a second element which is the \"generic segment\" to provide supplementary information. With option \"gseg\",0 this is avoided. Option \"neg\", 0,1: The default is \"neg\",0 With option \"neg\",0 Rabinovitch trick is used for negative hypothesis and thesis With option \"neg\",1 Difference of constructible sets is used instead. NOTE: The basering R, must be of the form Q[a][x], (a=parameters, x=variables), and should be defined previously. The ideals must be defined on R. KEYWORDS: Automatic Deduction; Automatic Demonstration; Geometric Theorems EXAMPLE: ADGT shows an example" { int i; int j; def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); Rx[1]=0; def D=ring(Rx); def RP=D+P; setring RR; list Lopt=#; int start=timer; // options int rep=0; int gseg=1; int neg=0; for(i=1;i<=size(Lopt) div 2;i++) { if(Lopt[2*i-1]=="rep"){rep=Lopt[2*i];} if(Lopt[2*i-1]=="gseg"){gseg=Lopt[2*i];} if(Lopt[2*i-1]=="neg"){neg=Lopt[2*i];} } // begin proc if(neg==0) { // Default option "neg",0 uses Rabinovitch for negative hyothesis and thesis // Option "neg",1 uses Diference of constructive sets for negative hyothesis and thesis if(equalideals(T1,ideal(1))) //nonnullT==0) { def F=H; for(i=1;i<=size(T);i++) { F[size(F)+1]=T[i]; } list G2; if(not(equalideals(H1,ideal(1)))) //nonnullH==1) { G2=grobcov(F,"nonnull",H1,"rep",1); } else { G2=grobcov(F,"rep",1); } } else { def F=H; for(i=1;i<=size(T1);i++) { F[size(F)+1]=T1[i]; } list G3=grobcov(F,"rep",1); def L0=Grob1Levels(G3); setring P; def LP=imap(RR,L0); def H1P=imap(RR,H1); if(not(equalideals(H1P,ideal(1)))) //nonnullH==1) { i=1; while(i<=size(LP)) { if(i mod 2==1) { LP[i]=intersect(LP[i],H1P); } i++; } } setring RR; L0=imap(P,LP); F=H; for(i=1;i<=size(T);i++) { F[size(F)+1]=T[i]; } def G2=grobcov(F,"nonnull",L0[1],"rep",1); list G1; int m=size(L0); int r=m div 2; if(m mod 2==0){r=r-1;} //"L0="; L0; for(i=1;i<=r;i++) { G1=grobcov(F,"null",L0[2*i],"nonnull",L0[2*i+1],"rep",1); for(j=1;j<=size(G1);j++) { if(not(equalideals(G1[j][1],ideal(1)))) { G2[size(G2)+1]=G1[j]; } } } } if(size(G2)==0) { list L; // L[1]=1; } else { list L=Levels(Grob1Levels(G2)); if(rep==0 and size(L)>0) { setring P; def LFP=imap(RR,L); // list LFP1=LFP; for(i=1;i<=size(LFP);i++) { LFP[i][2]=Prep(LFP[i][2][1],LFP[i][2][2]); } setring RR; L=imap(P,LFP); } if(equalideals(G2[1][3][1][1],0) and not(equalideals(G2[1][1],ideal(1)))) { list GL=G2[1]; list LL=GL[3]; if(rep==0) { setring P; def LLP=imap(RR,LL); LLP=Prep(LLP[1],LLP[2]); setring RR; LL=imap(P,LLP); GL[3]=LL; } if(gseg==1) { L[size(L)+1]=list("Generic segment",GL); } } } return(L); } else { if (neg==1) { def LL=ADGTDif(H,T,H1,T1,#); return(LL); } } } example { "EXAMPLE:"; echo = 2; // Determine the supplementary conditions // for the non-degenerate triangle A(-1,0), B(1,0), C(x,y) // to have an ortic non-degenerate isosceles triangle if(defined(R)){kill R;} ring R=(0,x,y),(x2,y2,x1,y1),lp; // Hypothesis H: the triangle A1(x1,y1), B1(x2,y2), C1(x,0), is the // orthic triangle of ABC ideal H=-y*x1+(x-1)*y1+y, (x-1)*(x1+1)+y*y1, -y*x2+(x+1)*y2-y, (x+1)*(x2-1)+y*y2; // Thesis T: the orthic triangle is isosceles ideal T=(x1-x)^2+y1^2-(x2-x)^2-y2^2; // Negative Hypothesis H1: ABC is non-degenerate ideal H1=y; // Negative Thesis T1: the orthic triangle is non-degenerate ideal T1=x*(y1-y2)-y*(x1-x2)+x1*y2-x2*y1; // Complementary conditions for the // Proposition (H and not H1) => (T and not T1) // to be true ADGT(H,T,H1,T1); // Now using diference of constructible sets for negative hypothesis and thesis ADGT(H,T,H1,T1,"neg",1); // The results are identical using both methods for the negative propositions // - Rabinovitch or // - DifConsLCSets // EXAMPLE // Automatic Theorem Proving. // The nine points circle theorem. // Vertices of the triangle: A(-2,0), B(2,0), C(2a,2b) // Heigth foot: A1(x1,y1), // Heigth foot: B1(x2,y2), // Heigth foot: C1(2a,0) // Middle point BC: A2(a+1,b) // Middle point CA: B2 (a-1,b) // Middle point AB: C2(0,0) // Ortocenter: O(2x0,2y0) // Middle point of A and O: A3(x0-1,y0) // Middle point of B and O: B3(x0+1,y0) // Middle point of C and O: C3(x0+a,y0+b) // Nine points circle: c:=(X-x3)^2+(Y-y3)^2-r2 if (defined(R1)){kill R1;} ring R1=(0,a,b),(x1,y1,x2,y2,x0,y0,x3,y3,r2),dp; short=0; ideal H=-x1*b+(a-1)*y1+2*b, (a-1)*x1+b*y1+2*a-2, -x2*b+(a+1)*y2-2*b, (a+1)*x2+b*y2-2*a-2, -x0*y1+(x1+2)*y0-y1, -x0*y2+(x2-2)*y0+y2; ideal T=(x1-2*x3)^2+(y1-2*y3)^2-r2, (a+1-2*x3)^2+(b-2*y3)^2-r2, (x0-1-2*x3)^2+(y0-2*y3)^2-r2, (x2-2*x3)^2+(y2-2*y3)^2-r2, (a-1-2*x3)^2+(b-2*y3)^2-r2, (x0+1-2*x3)^2+(y0-2*y3)^2-r2, (2*a-2*x3)^2+4*y3^2-r2, 4*x3^2+4*y3^2-r2, (x0+a-2*x3)^2+(y0+b-2*y3)^2-r2; ADGT(H,T,b,1); // Thus the nine point circle theorem is true for all real points excepting b=0. } static proc ADGTDif(ideal H,ideal T,ideal H1,ideal T1,list #) { int i; int j; def RR=basering; def Rx=ringlist(RR); def P=ring(Rx[1]); Rx[1]=0; def D=ring(Rx); def RP=D+P; setring RR; list Lopt=#; int start=timer; // options int rep=0; int gseg=1; list LLL; for(i=1;i<=size(Lopt) div 2;i++) { if(Lopt[2*i-1]=="rep"){rep=Lopt[2*i];} else{if(Lopt[2*i-1]=="gseg"){gseg=Lopt[2*i];}} } // begin proc def F=H; for(i=1;i<=size(T);i++) { F[size(F)+1]=T[i]; } list G2=grobcov(F,"rep",1); def L2=Grob1Levels(G2); //"L2="; L2; ideal FN=T1; if(not(equalideals(H1,ideal(1)))) { list FNH1=FN,H1; FN=intersectpar(FNH1); } if(not(equalideals(FN,ideal(1)))) { for(i=1;i<=size(FN);i++) { F[size(F)+1]=FN[i]; } list G3=grobcov(F,"rep",1); def L3=Grob1Levels(G3); //"L3="; L3; } def LL=DifConsLCSets(L2,L3); //"LL="; LL; LL=ConsLevels(LL); LL=Levels(LL); //"LL="; LL; if(rep==1){return(LL);} else { for(i=1;i<=size(LL);i++) { LL[i][2]=Prep(LL[i][2][1],LL[i][2][2]); } } return(LL); } //****************************** End ADGT ************************* ;