############################################## # Involutive Divisions # ############################################## ##############Janet Division############# Janet:=proc(u,U,Vars) local n,m,d,L,i,j,dd,V,v,Mult; #option trace; if degree(u)=0 then RETURN(Vars); fi; n:=nops(Vars); m:=nops(U); d:=[seq(max(seq(degree(U[j],Vars[i]),j=1..m)),i=1..n)]; Mult:=NULL; if degree(u,Vars[1])=d[1] then Mult:=Mult,Vars[1]; fi: for j from 2 to n do dd:=[seq(degree(u,Vars[l]),l=1..j-1)]; V:=NULL: for v in U do if [seq(degree(v,Vars[l]),l=1..j-1)]=dd then V:=V,v: fi: od: if degree(u,Vars[j])=max(seq(degree(v,Vars[j]), v in [V])) then Mult:=Mult,Vars[j]; fi: od: RETURN([Mult]); end: ########### Pommaret Division ############# LeftPommaret:=proc(u,U,Vars) local N,Ind,i; N:=NULL: Ind:=indets(u): for i from 1 to nops(Vars) while not Vars[i] in Ind do N:=N,Vars[i]: od: N:=N,Vars[i]: RETURN([N]); end: RightPommaret:=proc(u,U,Vars) local N,Ind,i; N:=NULL: Ind:=indets(u): for i from nops(Vars) by -1 to 1 while not Vars[i] in Ind do N:=N,Vars[i]: od: N:=N,Vars[i]: RETURN([N]); end: ########### Lexicographically Induced Division ############## with(Groebner): LID:=proc(u,U,Vars) local m,n,L,i,j,flag,v,Mult; Mult:=NULL; m,n:=nops(U),nops(Vars): for j from 1 to n do flag:=false: for v in U while not flag do if TestOrder(v,u,plex(op(Vars))) and degree(v,Vars[j])>degree(u,Vars[j]) then flag:=true: fi: od: if not flag then Mult:=Mult,Vars[j]; fi: od: RETURN([Mult]); end: ################## Thomas Division ################## Thomas:=proc(u,U,Vars) local d,L,i,j,Mult; Mult:=NULL; d:=[seq(max(seq(degree(m,Vars[i]),m in U)),i=1..nops(Vars))]; for j from 1 to nops(Vars) do if degree(u,Vars[j])=d[j] then Mult:=Mult,Vars[j]; fi: od: RETURN([Mult]); end: ############################ # Reduction # ############################ ############## Crit:=proc(p,gg,InvDivision) global T,tord,Vars,Q,Polys,c1,c2,LL,LMF; local InvDiv,flag,t,i,NM,y; #option trace; if InvDivision=J then InvDiv:=Janet: elif InvDivision=RP then InvDiv:=RightPommaret: elif InvDivision=LP then InvDiv:=LeftPommaret: elif InvDivision=L then InvDiv:=LID: elif InvDivision=T then InvDiv:=Thomas: fi: N:=NULL; g:=LM(Polys[gg[2]]); #for g in ALM do if LM(Polys[p[2]])*g=p[4]*LM(Polys[p[5]]) then c1:=c1+1; RETURN(true); elif divide(LM(p[4]*Polys[p[5]]),lcm(LM(Polys[p[2]]),g)) and p[4]*LM(Polys[p[5]])<>lcm(LM(Polys[p[2]]),g) then c2:=c2+1; RETURN(true); fi; #od; RETURN(false); end: ################ Normal Form Algorithm ################# with(Groebner): Normalform:=proc(pp,TermOrder,InvDivision) #option trace; local LMF,Vars,p,r,flag,g,u,Mult,InvDiv; global tord,T,Polys; if InvDivision=J then InvDiv:=Janet: elif InvDivision=RP then InvDiv:=RightPommaret: elif InvDivision=LP then InvDiv:=LeftPommaret: elif InvDivision=L then InvDiv:=LID: elif InvDivision=T then InvDiv:=Thomas: fi: LMF:=[seq(LeadingMonomial(Polys[v[1]],tord), v in T)]; Vars:=[op(TermOrder)]; p:=pp[4]*Polys[pp[5]]: r:=0; FF:=[seq(Polys[v[1]], v in T)]; j:=0: while p<>0 do j:=j+1: flag:=false: i:=1; while not flag and i<= nops(FF) do g:=FF[i]: if divide(LeadingMonomial(p,TermOrder),LeadingMonomial(g,TermOrder)) then u:=LeadingMonomial(p,TermOrder)/LeadingMonomial(g,TermOrder); Mult:=InvDiv(LeadingMonomial(g,TermOrder),LMF,Vars); if indets(u) subset {op(Mult)} then flag:=true: if j=1 then if Crit(pp,T[i],InvDivision) then RETURN(c1); fi: fi: p:=simplify(p-(LeadingCoefficient(p,TermOrder)/LeadingCoefficient(g,TermOrder))*u*g); else i:=i+1; fi: else i:=i+1; fi: od: if not flag then j:=j+1; r:=simplify(r+LeadingCoefficient(p,TermOrder)*LeadingMonomial(p,TermOrder)); p:=simplify(p-LeadingCoefficient(p,TermOrder)*LeadingMonomial(p,TermOrder)); fi: od: RETURN(r); end: ################ AutoReduce ################# AutoReduce:=proc(F,TermOrder,InvDivision) local InvDiv,H,i; #option trace; description "AutoReduce procedure Involutive-Reduces the list F of polynomials. For using Janet, LeftPommaret, RightPommaret, Lexicographic induced and Thomas divisions, write J, LP, RP, L and T, respectively for InvDivision."; if InvDivision=J then InvDiv:=Janet: elif InvDivision=RP then InvDiv:=RightPommaret: elif InvDivision=LP then InvDiv:=LeftPommaret: elif InvDivision=L then InvDiv:=LID: elif InvDivision=T then InvDiv:=Thomas: fi: H:=F; for i from 1 to nops(F) do H[i]:=AutoReduction(H[i],[op({op(H)} minus {H[i],0})],TermOrder,InvDivision); od: RETURN([op({op(H)} minus {0})]); end: Describe(AutoReduce); ######################### tidy:=proc(p,q) global Polys,tord; if LM(Polys[p[1]])<>LM(Polys[q[1]]) then RETURN(TestOrder(LM(Polys[p[1]]),LM(Polys[q[1]]),tord)); else RETURN(evalb(p[2]0 then RETURN(LeadingMonomial(f,tord)); fi: RETURN(0); end: ################# LC:=proc(f) global tord; RETURN(LeadingCoefficient(f,tord)); end: ##################### ##################### InvBasis:=proc(Ideal,Tord,InvDivision) #option trace; global tord,T,Q,Polys,Vars,c1,c2,rdz,LMF; local InvDiv,F,p,TT,q,h,NM,np,j,firsttime,firstbytes,secondtime,secondbytes,deg,isr; if InvDivision=J then InvDiv:=Janet: elif InvDivision=RP then InvDiv:=RightPommaret: elif InvDivision=LP then InvDiv:=LeftPommaret: elif InvDivision=L then InvDiv:=LID: elif InvDivision=T then InvDiv:=Thomas: fi: firsttime,firstbytes:=kernelopts(cputime,bytesused); c1:=0:c2:=0: rdz:=0; tord:=Tord; Vars:=[op(Tord)]; F:=AutoReduce(Ideal,tord,InvDiv); F:=sort(F,(a,b)->TestOrder(a,b,tord)); isr:=0: n:=nops(F): deg:=0: Polys:=Array(1..n,i->F[i]); T:=[[1,1,{},1,1]]; Q:=[seq([i,i,{},1,i],i=2..nops(F))]; Q:=sort(Q,tidy); while nops(Q)<>0 do p:=Q[1]; deg:=max(deg,degree(Polys[p[1]])): Q:=Q[2..-1]; h:=Normalform(p,tord,InvDivision); if h=0 then rdz:=rdz+1: fi: if h=0 and p[1]=p[2] then SS:=NULL; S:=Q: for s in S do if s[2]<>p[1] then SS:=SS,s; else isr:=isr+1: fi; od: Q:=[SS]; elif h<>0 and h<>c1 and LM(h)<>LM(Polys[p[1]]) then Polys[p[1]]:=h; T:=[op(T),[p[1],p[1],{},1,p[1]]]; TT:=NULL; for q in T do if divide(LM(Polys[q[1]]), LM(Polys[p[1]])) and LM(Polys[p[1]]) <> LM(Polys[q[1]]) then Q:=[op(Q),q]; else TT:=TT,q; fi: od: T:=[TT]; for i from 1 to nops(T) do q:=T[i]; NM:={op(Vars)} minus {op(InvDiv(LM(Polys[q[1]]),[seq(LM(Polys[w[1]]), w in T)],Vars))}; for y in NM minus q[3] do n:=n+1: Polys(n):=expand(y*Polys[q[1]]); Q:=[op(Q),[n,q[2],{},y,q[1]]]; od: T[i][3]:=(NM minus q[3]) union (NM intersect q[3]); od: elif h<>0 and h<>c1 and LM(h)=LM(Polys[p[1]]) then Polys[p[1]]:=h; T:=[op(T),p]; for i from 1 to nops(T) do q:=T[i]; NM:={op(Vars)} minus {op(InvDiv(LM(Polys[q[1]]),[seq(LM(Polys[w[1]]), w in T)],Vars))}; for y in NM minus q[3] do n:=n+1: Polys(n):=expand(y*Polys[q[1]]); Q:=[op(Q),[n,q[2],{},y,q[1]]]; od: T[i][3]:=(NM minus q[3]) union (NM intersect q[3]); od: fi; Q:=sort(Q,tidy); od: secondtime,secondbytes:=kernelopts(cputime,bytesused); T:=[seq(Polys[tt[1]], tt in T)]; appendto(testVar); printf("%-1s %1s %1s %1s : %3a %3a\n",The, cpu, time, is,secondtime-firsttime,(sec)): printf("%-1s %1s %1s : %3a %3a\n",The,used,memory,secondbytes-firstbytes,(bytes)): printf("%-1s %1s %1s : %3g\n",Num,of,Rds,rdz): printf("%-1s %1s %1s : %3a\n",Num,of,criteria,[c1,c2]): printf("%-1s %1s %1s : %3a\n",Num,of,ISR,isr): printf("%-1s %1s %1s : %3a\n",Max,of,degree,deg): printf("%-1s %1s %1s : %3a\n",Num,of,polynomials,nops(T)): print("IsInvolutive",IsInvolutive(LM(T),tord,InvDivision)); print("IsGrobner",IsGrobner(Ideal,T,Tord)); appendto(terminal); RETURN(T); end: ######################## # IsGrobner # ######################## ######################## IsGrobner:=proc(A,H,T) local s,j,S,p,F,L,LL; F:=H; while member(0, F, 'p') do F:=subsop(p=NULL,F); unassign('p'); od; L:=LeadingMonomial(F, T): LL:=LeadingMonomial(Basis(A,T),T): if evalb(LeadingMonomial(, T)<>LeadingMonomial(, T)) then RETURN(false); fi: for s from 1 to nops(A) do if evalb(Reduce(A[s],F,T)<>0) then RETURN(false); fi; od; RETURN(true); end: ############################ ############################# # Involutive Test # ############################# AutoReduction:=proc(f,F,TermOrder,InvDivision) #option trace; local LMF,Vars,p,r,flag,g,u,Mult,InvDiv; global tord; if InvDivision=J then InvDiv:=Janet: elif InvDivision=RP then InvDiv:=RightPommaret: elif InvDivision=LP then InvDiv:=LeftPommaret: elif InvDivision=L then InvDiv:=LID: elif InvDivision=T then InvDiv:=Thomas: fi: LMF:=[seq(LeadingMonomial(v,tord), v in F)]; Vars:=[op(TermOrder)]; p:=f: r:=0; FF:={op(F)} minus {0}; while p<>0 do flag:=false:i:=1; while not flag and i<= nops(FF) do g:=FF[i]: if divide(LeadingMonomial(p,TermOrder),LeadingMonomial(g,TermOrder)) then u:=LeadingMonomial(p,TermOrder)/LeadingMonomial(g,TermOrder); Mult:=InvDiv(LeadingMonomial(g,TermOrder),LMF,Vars); if indets(u) subset {op(Mult)} then flag:=true: p:=simplify(p-(LeadingCoefficient(p,TermOrder)/LeadingCoefficient(g,TermOrder))*u*g); else i:=i+1; fi: else i:=i+1; fi: od: if not flag then r:=simplify(r+LeadingCoefficient(p,TermOrder)*LeadingMonomial(p,TermOrder)); p:=simplify(p-LeadingCoefficient(p,TermOrder)*LeadingMonomial(p,TermOrder)); fi: od: RETURN(r); end: ############################# IsInvolutive:=proc(F,TermOrder,InvDivision) local InvDiv,Vars,LMF,f,NM,u; #option trace; if InvDivision=J then InvDiv:=Janet: elif InvDivision=RP then InvDiv:=RightPommaret: elif InvDivision=LP then InvDiv:=LeftPommaret: elif InvDivision=L then InvDiv:=LID: elif InvDivision=T then InvDiv:=Thomas: fi: Vars:=[op(TermOrder)]; LMF:=[seq(LeadingMonomial(f,TermOrder), f in F)]; for f in F do NM:={op(Vars)} minus {op(InvDiv(LeadingMonomial(f,TermOrder),LMF,Vars))}; for u in NM do if AutoReduction(u*f,F,TermOrder,InvDivision)<>0 then print(u,f); RETURN(false); fi: od: od: RETURN(true); end: writeto(testVar); "This is the benchmark for VarGerdt algorithm."; writeto(terminal); appendto(testVar); print("####################Example 1"); appendto(terminal); F:=[8*x^2-2*x*y-6*x*z+3*x+3*y^2-7*y*z+10*y+10*z^2-8*z-4,10*x^2-2*x*y+6*x*z-6*x+9*y^2-y*z-4*y-2*z^2+5*z-9 ,5*x^2+8*x*y+4*x*z+8*x+9*y^2-6*y*z+2*y-z^2-7*z+5]: InvBasis(F,tdeg(x,y,z),J): appendto(testVar); print("####################Example 2"); appendto(terminal); F:=[ttt+vvv-aaa,xxx+yyy+zzz+ttt-uuu-www-aaa,xxx*zzz+yyy*zzz+xxx*ttt+zzz*ttt-uuu*www-uuu*aaa-www*aaa]: InvBasis(F,tdeg(xxx,yyy,zzz,ttt,uuu,vvv,www,aaa),J): appendto(testVar); print("####################Example 3"); appendto(terminal); F:=[a*y^2+b*x^3+c, 2*a*y, 3*b*x^2]: InvBasis(F,tdeg(x,y,a,b,c),J): appendto(testVar); print("####################Example 4"); appendto(terminal); F:=[x^2+b*y^2+2*c*x*y+2*d*x+2*e*y+f, x+c*y+d, b*y+c*x+e]: InvBasis(F,tdeg(x,y,a,b,c,e,f),J): appendto(testVar); print("####################Example 5"); appendto(terminal); F:=[a*x^3*y+c*x*y^2, x^2*y+3*d*y, c*x^2+b*x*y]: InvBasis(F,tdeg(x,y,a,b,c,d),J): appendto(testVar); print("####################Example 6"); appendto(terminal); F:=[a*x^2*y+b*x+y^3, a*x^2*y+b*x*y, y^2+b*x^2*y+c*x*y]: InvBasis(F,tdeg(x,y,a,b,c),J): appendto(testVar); print("####################HAAS 3"); appendto(terminal); F := [x^6+a*y^3-y, y^6+b*x^3-x, 36*y^5*x^5-9*x^2*b*a*y^2+3*a*y^2+3*b*x^2-1]: InvBasis(F,tdeg(x,y,a,b),J): appendto(testVar); print("####################Liu"); appendto(terminal); F:=[y*z-y*t0-x*h+a*h, z*t0-z*x-y*h+a*h, t0*x-y*t0-z*h+a*h, x*y-z*x-t0*h+a*h]: InvBasis(F,tdeg(z, y, x, t0, a, h),J): appendto(testVar); print("####################Cyclic 5"); appendto(terminal); F:=[a*b*c*d*e-1, a*b*c*d + a*b*c*e + a*b*d*e + a*c*d*e + b*c*d*e, a*b*c + a*b*e + a*d*e + b*c*d + c*d*e, a*b + a*e + b*c + c*d + d*e, a + b + c + d + e]: InvBasis(F,tdeg(a,b,c,d,e),J): appendto(testVar); print("####################Noon"); appendto(terminal); F := [10*x1^2*x4+10*x2^2*x4+10*x3^2*x4-11*x4*h^2+10*h^3, 10*x1^2*x3+10*x2^2*x3+10*x3*x4^2-11*x3*h^2+10*h^3, 10*x1*x2^2+10*x1*x3^2+10*x1*x4^2-11*x1*h^2+10*h^3, 10*x1^2*x2+10*x2*x3^2+10*x2*x4^2-11*x2*h^2+10*h^3]: InvBasis(F, tdeg(x1, x2, x3, x4, h), J): appendto(testVar); print("####################Weispfenning94"); appendto(terminal); F := [y^4+x*y^2*z+x^2*h^2-2*x*y*h^2+y^2*h^2+z^2*h^2, x*y^4+y*z^4-2*x^2*y*h^2-3*h^5, -y^2*x^3+x*y*z^3+y^4*h+x*y^2*z*h-2*x*y*h^3]: InvBasis(F, tdeg(x,y,z,h), J): appendto(testVar); print("####################Katsura5"); appendto(terminal); F := [2*ax^2+2*ay^2+2*az^2+2*at^2+2*au^2+av^2-av, ax*ay+ay*az+2*az*at+2*at*au+2*au*av-au, 2*ax*az+2*a*y*at+2*a*z*au+au^2+2*at*av-at, 2*ax*at+2*ay*au+2*at*au+2*az*av-az, at^2+2*ax*av+2*ay*av+2*az*av-ay, 2*ax+2*ay+2*az+2*at+2*au+av-1]: InvBasis(F, tdeg(ax, ay, az, at, au, av), J): appendto(testVar); print("####################Eco7"); appendto(terminal); F := [(x1+x1*x2+x2*x3+x3*x4+x4*x5+x5*x6)*x7-1, (x2+x1*x3+x2*x4+x3*x5+x4*x6)*x7-2, (x3+x1*x4+x2*x5+x3*x6)*x7-3, (x4+x1*x5+x2*x6)*x7-4, (x5+x1*x6)*x7-5, x6*x7-6, x1+x2+x3+x4+x5+x6+1]: InvBasis(F, tdeg(x1, x2, x3, x4,x5,x6,x7), J): appendto(testVar); print("####################Lichtblau"); appendto(terminal); F := [x-110*t^2+495*t^3-1320*t^4+2772*t^5-5082*t^6+7590*t^7-8085*t^8+5555*t^9-2189*t^10+374*t^11, y-22*t+110*t^2-330*t^3+1848*t^5-3696*t^6+3300*t^7-1650*t^8+550*t^9-88*t^10-22*t^11]: InvBasis(F, tdeg(x,y,t), J): appendto(testVar); print("####################Cyclic6"); appendto(terminal); F:=[a*b*c*d*e*f -1, a*b*c*d*e + a*b*c*d*f + a*b*c*e*f + a*b*d*e*f + a*c*d*e*f + b*c*d*e*f, a*b*c*d + a*b*c*f + a*b*e*f + a*d*e*f + b*c*d*e + c*d*e*f, a*b*c + a*b*f + a*e*f + b*c*d + c*d*e + d*e*f, a*b + a*f + b*c + c*d + d*e + e*f, a + b + c + d + e + f]: InvBasis(F,tdeg(a,b,c,d,e,f),J): appendto(testVar); print("####################Katsura6"); appendto(terminal); F:=[ 1*x1+2*x2+2*x3+2*x4+2*x5+2*x6+2*x7-1, 2*x4*x3+2*x5*x2+2*x6*x1+2*x7*x2-1*x6, 1*x3^2+2*x4*x2+2*x5*x1+2*x6*x2+2*x7*x3-1*x5, 2*x3*x2+2*x4*x1+2*x5*x2+2*x6*x3+2*x7*x4-1*x4, 1*x2^2+2*x3*x1+2*x4*x2+2*x5*x3+2*x6*x4+2*x7*x5-1*x3, 2*x2*x1+2*x3*x2+2*x4*x3+2*x5*x4+2*x6*x5+2*x7*x6-1*x2, 1*x1^2+2*x2^2+2*x3^2+2*x4^2+2*x5^2+2*x6^2+2*x7^2-1*x1 ]: InvBasis(F, tdeg(x1, x2, x3, x4,x5,x6,x7), J): appendto(testVar); print("####################Sturmfels and Eisenbud"); appendto(terminal); F := [ss*uu+bb*vv, tt*uu+bb*ww, tt*vv+ss*ww, ss*xx+bb*yy, tt*xx+bb*zz, tt*yy+ss*zz, vv*xx+uu*yy, ww*xx+uu*zz, ww*yy+vv*zz]: InvBasis(F, tdeg(bb, xx, yy, zz, ss, tt, uu, vv, ww), J):