//////////////////////////////////////////////////////////////////// version="$Id$"; category="Algebraic Geometry"; // summary description of the library info=" LIBRARY: JMBTest.lib A library for Singular which performs JM basis test. AUTHOR: Michela Ceria, email: michela.ceria@unito.it SEE ALSO: JMSConst_lib KEYWORDS: J-marked schemes OVERVIEW: The library performs the J-marked basis test, as described in [CR], [BCLR]. Such a test is performed via the criterion explained in [BCLR], concerning Eliahou-Kervaire polynomials (EK from now on). We point out that all the polynomials are homogeneous and they must be arranged by degree. The fundamental steps are the following:@* -construct the Vm polynomials, via the algorithm VConstructor explained in [CR];@* -construct the Eliahou-Kervaire polynomials defined in [BCLR];@* -reduce the Eliahou-Kervaire polynomials using the Vm's;@* -if it exist an Eliahou-Kervaire polynomial such that its reduction mod Vm is different from zero, the given one is not a J-Marked basis. The algorithm terminates only if the ordering is rp. Anyway, the number of reduction steps is bounded. REFERENCES: [CR] Francesca Cioffi, Margherita Roggero,Flat Families by Strongly Stable Ideals and a Generalization of Groebner Bases, J. Symbolic Comput. 46, 1070-1084, (2011).@* [BCLR] Cristina Bertone, Francesca Cioffi, Paolo Lella, Margherita Roggero, Upgraded methods for the effective computation of marked schemes on a strongly stable ideal, Journal of Symbolic Computation (2012), http://dx.doi.org/10.1016/j.jsc.2012.07.006 @* PROCEDURES: Minimus(ideal) minimal variable in an ideal Maximus(ideal) maximal variable in an ideal StartOrderingV(list,list) ordering of polynomials as in [BCLR] TestJMark(list) tests whether we have a J-marked basis "; LIB "qhmoduli.lib"; LIB "monomialideal.lib"; LIB "ring.lib"; //////////////////////////////////////////////////////////////////// proc mod_init() "USAGE: mod_init(); RETURN: struct: jmp EXAMPLE: example mod_init; shows an example" { newstruct("jmp", "poly h, poly t"); } example { "EXAMPLE:"; echo = 2; mod_init(); } //////////////////////////////////////////////////////////////////// proc Terns(list G, int c) "USAGE: Terns(G,c); G list, c int RETURN: list: T NOTE: Input is a list of J-marked polynomials (arranged by degree) and an integer. EXAMPLE: example Terns; shows an example" { list T=list(); int z; for(int k=1; k<=size(G[c]);k=k+1) { //Loop on G[c] making positions of polynomials in G[c] z=size(T); T=insert(T,list(1,c,k) ,size(T)); } return(T); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,z), rp; jmp r1; r1.h=z^3; r1.t=poly(0); jmp r2; r2.h=z^2*y; r2.t=poly(0); jmp r3; r3.h=z*y^2 ; r3.t=-x^2*y; jmp r4; r4.h=y^5; r4.t=poly(0); list G2F=list(list(r1,r2,r3),list(r4)); Terns(G2F, 1); Terns(G2F, 2); } //////////////////////////////////////////////////////////////////// proc VConst(list G, int c) "USAGE: VConst(G, c); G list, c int RETURN: list: V NOTES: this procedure computes the Vm polynomials following the algorithm in [CR],but it only keeps in memory the monomials by which the G's must be multplied and their positions. EXAMPLE: example VConst; shows an example" { jmp f=G[1][1]; int aJ=deg(f.h); // minimal degree of polynomials in G //print(aJ); list V=list(); V[1]=Terns(G,1); // V[1]=G[1] (keeping in memory only [head, position]) //print(c-aJ+1); int i; int j; int m; list OO; jmp p; for(m=2; m<=c-aJ+1; m=m+1) { //print("entro nel form"); if(m>size(G)) {V[m]=list(); //If we have not G[m] we insert a list() //print("vuota prima"); } else {V[m]=Terns(G,m); //print("piena prima"); } for(i=1; imax){max=L[i];} } return(max); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,z), rp; ideal I=y,x,z; Maximus(I); } //////////////////////////////////////////////////////////////////// proc GJmpMins(jmp P, jmp Q) "USAGE: GJmpMins(P,Q); P jmp, Q jmp RETURN: int: d EXAMPLE: example GJmpMins; shows an example" { int d=1; //-1=lower, 0=equal, 1=higher //At the beginning suppose Q is higher if(deg(P.h)B[1]) { //the ordering MUST be rp //print("Maggiore per Lex"); d=1; } } return(d); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,z), rp; jmp r1; r1.h=z^3; r1.t=poly(0); jmp r2; r2.h=z^2*y; r2.t=poly(0); jmp r3; r3.h=z*y^2 ; r3.t=-x^2*y; jmp r4; r4.h=y^5; r4.t=poly(0); list G2F=list(list(r1,r2,r3),list(r4)); TernCompare([1,1,1],[x,1,1],G2F); } //////////////////////////////////////////////////////////////////// proc MinOfV(list V, list G) "USAGE: Minimal(V,G); V list, G list RETURN: int: R NOTE: Input=lista(terne), G. EXAMPLE: example Minimal; shows an example" { //Minimal element for a given degree list R=list(); list MIN=V[1]; int h=1; int i; for(i=2; i<=size(V); i++) { //I consider the first as minimum //If I find something smaller I change minimum if(TernCompare(V[i],MIN,G)<=0) { MIN=V[i]; h=i; } } //Return: [minimum,position of the minimum] R=MIN,h; return(R); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,z), rp; jmp r1; r1.h=z^3; r1.t=poly(0); jmp r2; r2.h=z^2*y; r2.t=poly(0); jmp r3; r3.h=z*y^2 ; r3.t=-x^2*y; jmp r4; r4.h=y^5; r4.t=poly(0); list G2F=list(list(r1,r2,r3),list(r4)); MinOfV(VConst(G2F,4,basering)[1],G2F); } //////////////////////////////////////////////////////////////////// proc OrderingV(list V,list G,list R) "USAGE: OrderingV(V,G,R); V list, G list, R list RETURN: list: R NOTE: Input: Vm,G,emptylist EXAMPLE: example OrderingV; shows an example" { //Order V[m] //R will contain results but at the beginning it is empty list M=list(); if(size(V)==1) { R=insert(R,V[1],size(R)); } else { M=MinOfV(V,G); R=insert(R,M[1],size(R)); V=delete(V,M[2]); //recursive call R=OrderingV(V,G,R); } return(R); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,z), rp; jmp r1; r1.h=z^3; r1.t=poly(0); jmp r2; r2.h=z^2*y; r2.t=poly(0); jmp r3; r3.h=z*y^2; r3.t=-x^2*y; jmp r4; r4.h=y^5; r4.t=poly(0); list G2F=list(list(r1,r2,r3),list(r4)); OrderingV(VConst(G2F,4,basering)[1],G2F,list()); } //////////////////////////////////////////////////////////////////// proc StartOrderingV(list V,list G) "USAGE: StartOrdina(V,G); V list, G list RETURN: list: R NOTE: Input Vm,G. This procedure uses OrderingV to get the ordered polynomials as in [BCLR]. EXAMPLE: example StartOrderingV; shows an example" { return(OrderingV(V,G, list())); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,z), rp; jmp r1; r1.h=z^3; r1.t=poly(0); jmp r2; r2.h=z^2*y; r2.t=poly(0); jmp r3; r3.h=z*y^2; r3.t=-x^2*y; jmp r4; r4.h=y^5; r4.t=poly(0); list G2F=list(list(r1,r2,r3),list(r4)); StartOrderingV(VConst(G2F,4,basering)[1],G2F); } //////////////////////////////////////////////////////////////////// proc Multiply(list L, list G) "USAGE: moltiplica(L,G); L list, G list RETURN: jmp: K NOTE: Input: a 3-ple,G. It performs the product associated to the 3-uple. EXAMPLE: example Multiply; shows an example" { jmp g=G[L[2]][L[3]]; jmp K; K.h=L[1]*g.h; K.t=L[1]*g.t; return(K); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,z), rp; list P=x^2,1,1; jmp r1; r1.h=z^3; r1.t=poly(0); jmp r2; r2.h=z^2*y; r2.t=poly(0); jmp r3; r3.h=z*y^2 ; r3.t=-x^2*y; jmp r4; r4.h=y^5; r4.t=poly(0); list G2F=list(list(r1,r2,r3),list(r4)); Multiply(P,G2F); } //////////////////////////////////////////////////////////////////// proc IdealOfV(list V) "USAGE: IdealOfV(V); V list RETURN: ideal: I NOTES: this procedure takes a list of Vm's of a certain degree and construct their ideal, multiplying the head by the weighted variable t. EXAMPLE: example IdealOfV; shows an example" { ideal I=0; int i; if (size(V)!=0) { list M=list(); jmp g; for(i=1; i<= size(V); i++) { g=V[i]; g.h=t*g.h; M[i]=g.h+g.t; } I=M[1..size(M)]; //print("IdealOfV"); //I=std(I); } return(I); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,z,t), rp; jmp r1; r1.h=z^3; r1.t=poly(0); jmp r2; r2.h=z^2*y; r2.t=poly(0); jmp r3; r3.h=z*y^2 ; r3.t=-x^2*y; jmp r4; r4.h=y^5; r4.t=poly(0); list G2F=list(list(r1,r2,r3),list(r4)); IdealOfV(G2F[1]); } //////////////////////////////////////////////////////////////////// proc NewWeight(int n) "USAGE: NewWeight(n); n int RETURN: intvec: u EXAMPLE: example NewWeight; shows an example" { intvec u=0; u[n]=1; return(u); } example { "EXAMPLE:"; echo = 2; NewWeight(3); } //////////////////////////////////////////////////////////////////// proc FinalVm(list V1 , list G1 , r) "USAGE: FinalVm(V1, G1, r); V1 list, G1 list , r RETURN: intvec: u EXAMPLE: example NewWeight; shows an example" { //multiply and reduce, degree by degree intvec u=NewWeight(nvars(r)+1); list L=ringlist(r); L[2]=insert(L[2],"t",size(L[2])); //print(L[2]); list ordlist="a",u; L[3]=insert(L[3],ordlist,0); def H=ring(L); //print(V1); //print(G1); list M=list(); jmp p; list N; poly q; poly s; int i; int j; for(i=1; i<=size(G1); i++) { N=list(); for(j=1; j<=size(G1[i]); j++) { p=G1[i][j]; q=p.h; s=p.t; N[j]=list(q,s); } M[i]=N; } p.h=poly(0); p.t=poly(0); setring H; list R=list(); list S=list(); //print("anello definito"); def V=imap(r,V1); //def G=imap(r,G1); //print(V); def MM=imap(r,M); list G=list(); list N=list(); for(i=1; i<=size(MM); i++) { for(j=1; j<=size(MM[i]); j++) { p.h=MM[i][j][1]; p.t=MM[i][j][2]; N[j]=p; } G[i]=N; } ideal I=0; jmp LL; jmp UU; for(i=1; i<=size(V);i++) { R[i]=list(); S[i]=list(); I=0; for(j=1;j<=size(V[i]); j++) { LL=Multiply(V[i][j],G); LL.t=reduce(t*LL.t,I); //I only reduce the tail LL.t=subst(LL.t,t,1); S[i]=insert(S[i],LL,size(S[i])); LL.h=t*LL.h; R[i]=insert(R[i],LL,size(R[i])); UU=R[i][j]; I=I+ideal(UU.h+UU.t); attrib(I,"isSB",1); } } list M=list(); poly q; poly s; for(i=1; i<=size(S); i++) { N=list(); for(j=1; j<=size(S[i]); j++) { p=S[i][j]; q=p.h; s=p.t; N[j]=list(q,s); } M[i]=N; } p.h=poly(0); p.t=poly(0); setring r; def MM=imap(H,M); list MMM=list(); for(i=1; i<=size(MM); i++) { N=list(); for(j=1; j<=size(MM[i]); j++) { p.h=MM[i][j][1]; p.t=MM[i][j][2]; N[j]=p; } MMM[i]=N; } return(MMM); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,z), rp; jmp r1; r1.h=z^3; r1.t=poly(0); jmp r2; r2.h=z^2*y; r2.t=poly(0); jmp r3; r3.h=z*y^2 ; r3.t=-x^2*y; jmp r4; r4.h=y^5; r4.t=poly(0); list G2F=list(list(r1,r2,r3),list(r4)); FinalVm(VConst(G2F,6,r) , G2F, r); } //////////////////////////////////////////////////////////////////// proc ConstructorMain(list G, int c,r) "USAGE: Costruttore(G,c); G list, c int RETURN: list: R NOTE: At the end separated by degree. EXAMPLE: example Costruttore; shows an example" { list V=list(); V= VConst(G,c); //print("VConst"); //V non ordered list L=list(); list R=list(); int i; // head, position //order the different degrees for(i=1; i<=size(V); i++) { L[i]=StartOrderingV(V[i], G); } //multiply and reduce //print("Ordinare"); R=FinalVm(L, G, r); //print("FinalVm"); return(R); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,z), rp; jmp r1; r1.h=z^3; r1.t=poly(0); jmp r2; r2.h=z^2*y; r2.t=poly(0); jmp r3; r3.h=z*y^2 ; r3.t=-x^2*y; jmp r4; r4.h=y^5; r4.t=poly(0); list G2F=list(list(r1,r2,r3),list(r4)); ConstructorMain(G2F,6,r); } //////////////////////////////////////////////////////////////////// proc EKCouples(jmp A, jmp B) "USAGE: CoppiaEK(A,B); A list, B list RETURN: list: L NOTE: At the end the monomials involved by EK. EXAMPLE: example EKCouples; shows an example" { poly E; list L=0,0; string s=varstr(basering); list VVV=varstr(basering); //L will contain results poly h=Minimus(variables(A.h)); //print(h); int l=findvars(h,1)[2][1]; if(l!=nvars(basering)) { //print("vero"); //print(l); for(int j=l+1;j<=nvars(basering); j++) { //print("entrata"); //print(var(j)); E=var(j)*A.h/B.h; //Candidate for * product //print(E); if(E!=0) { //print("primo if passato"); if(Minimus(variables(B.h))>=Maximus(variables(E))) { //Does it work with * ? //print("secondo if passato"); L[1]=j; L[2]=E; break; } } } } return (L); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,z), rp; jmp A; A.h=y*z^2; A.t=poly(0); jmp B; B.h=y^2*z; B.t=poly(0); EKCouples(A,B); EKCouples(B,A); } //////////////////////////////////////////////////////////////////// proc EKPolys(list G) "USAGE: PolysEK(G); G list RETURN: list: EK, list: D NOTE: At the end EK polynomials and their degrees EXAMPLE: example PolysEK; shows an example" { list D=list(); list C=list(); list N=0,0; list EK=list(); int i; int j; int k; int l; jmp p; for(i=1; i<=size(G); i++) { for(j=1; j<=size(G[i]); j++) { for(k=1; k<=size(G); k++) { for(l=1; l<=size(G[k]); l++) { if(i!=k||j!=l) { //Loop on polynomials C=EKCouples(G[i][j], G[k][l]); //print("coppia"); if(C[2]!=0) { C=insert(C,list(i,j,k,l),size(C)); EK=insert(EK,C,size(EK)); p=G[k][l]; D=insert(D,deg(C[2]*p.h),size(D)); } } } } } } //Double Return return(EK, D); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,z), rp; jmp r1; r1.h=z^3; r1.t=poly(0); jmp r2; r2.h=z^2*y; r2.t=poly(0); jmp r3; r3.h=z*y^2; r3.t=-x^2*y; jmp r4; r4.h=y^5; r4.t=poly(0); list G2F=list(list(r1,r2,r3),list(r4)); EKPolys(G2F); } //////////////////////////////////////////////////////////////////// proc EKPolynomials(list EK, list G) "USAGE: EKPolynomials(EK,G); EK list, G list RETURN: list: p NOTE: At the end I obtain the EK polynomials and their degrees. EXAMPLE: example SpolyEK; shows an example" { jmp u=G[EK[3][1]][EK[3][2]]; jmp q=G[EK[3][3]][EK[3][4]]; return(var(EK[1])*(u.h+u.t)-EK[2]*(q.h+q.t)); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,z), rp; jmp r1; r1.h=z^3; r1.t=poly(0); jmp r2; r2.h=z^2*y; r2.t=poly(0); jmp r3; r3.h=z*y^2; r3.t=-x^2*y; jmp r4; r4.h=y^5; r4.t=poly(0); list G2F=list(list(r1,r2,r3),list(r4)); list EK,D=EKPolys(G2F); EKPolynomials(EK[1],G2F); } //////////////////////////////////////////////////////////////////// proc TestJMark(list G1,r) "USAGE: TestJMark(G); G list RETURN: int: i NOTE: This procedure performs J-marked basis test.@* The input is a list of J-marked polynomials (jmp) arranged by degree, so G1 is a list of list.@* The output is a boolean evaluation: True=1/False=0 EXAMPLE: example TestJMark; shows an example" {int flag=1; if(size(G1)==1 && size(G1[1])==1) { //Hypersurface print("Only One Polynomial"); flag=1; } else { int d=0; list EK,D=EKPolys(G1); //print("PolysEK"); //I found EK couples int massimo=Max(D); list V1=ConstructorMain(G1,massimo,r); //print("Costruttore"); //print(V1); jmp mi=V1[1][1]; int minimo=Min(deg(mi.h)); intvec u=NewWeight(nvars(r)+1); list L=ringlist(r); L[2]=insert(L[2],"t",size(L[2])); //print(L[2]); list ordlist="a",u; L[3]=insert(L[3],ordlist,0); def H=ring(L); list JJ=list(); jmp pp; jmp qq; int i; int j; list NN; for(i=size(V1);i>0;i--) { NN=list(); for(j=size(V1[i]);j>0;j--) { //print(j); pp=V1[i][j]; NN[j]=list(pp.h,pp.t); } //print(NN); JJ[i]=NN; //print(JJ[i]); //print(i); } //print(JJ); list KK=list(); list UU=list(); //jmp qq; for(i=size(G1);i>0;i--) { for(j=size(G1[i]);j>0;j--) { //print(j); qq=G1[i][j]; UU[j]=list(qq.h,qq.t); } //print(UU); KK[i]=UU; } setring H; //I defined the new ring with the weighted //variable t poly p; //print("anello definito"); def JJJ=imap(r,JJ); def EK=imap(r,EK); //print(flag); //imap(r,D); list V=list(); jmp fp; //int i; //int j; list N; for(i=size(JJJ); i>0; i--) { N=list(); for(j=size(JJJ[i]); j>0; j--) { fp.h=JJJ[i][j][1]; fp.t=JJJ[i][j][2]; N[j]=fp; } V[i]=N; } //print(V); def KKJ=imap(r,KK); list G=list(); list U=list(); for(i=1; i<=size(KKJ); i++) { for(j=1; j<=size(KKJ[i]); j++) { fp.h=KKJ[i][j][1]; fp.t=KKJ[i][j][2]; U[j]=fp; } G[i]=U; } // print(V); //print(G); //I imported in H everithing I need poly q; ideal I; for(j=1; j<=size(EK);j++) { d=D[j]; p=EKPolynomials(EK[j],G); //print("arrivo"); I=IdealOfV(V[d-minimo+1]); attrib(I,"isSB",1); //print(I); q=reduce(t*p,I); //print(I[1]); //print(t*p); q=subst(q,t,1); //I reduce all the EK polynomials // q=RiduzPoly(V[d-minimo+1], p); if(q!=0) { //check whether reduction is 0 print("NOT A BASIS"); flag=0; break; } } } //print(flag); setring r; //typeof(flag); return(flag); } example { "EXAMPLE:"; echo = 2; ring r=0, (x,y,z), rp; jmp r1; r1.h=z^3; r1.t=poly(0); jmp r2; r2.h=z^2*y; r2.t=poly(0); jmp r3; r3.h=z*y^2 ; r3.t=-x^2*y; jmp r4; r4.h=y^5; r4.t=poly(0); list G2F=list(list(r1,r2,r3),list(r4)); TestJMark(G2F,r); }