############################################## # 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 # ############################ ###### This proc checks involutive divisibility ###### InvDivide2:=proc(u,v,L,V,InvDivision) global tord,Vars; local InvDiv,uu,M; #option trace; if v=0 then RETURN(false); fi; 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: if divide(u,v) then uu:=u/v; M:=InvDiv(v,L,Vars); if indets(uu) subset {op(M)} then RETURN(true); else RETURN(false); fi; else RETURN(false); fi; RETURN(false); end: ################ AutoReduced Algorithm ################# with(Groebner): 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:=F:#{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: ################ 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: ##################### Criteria2:=proc(p,gg,InvDivision) global T,tord,Vars,Q,Polys,c1,c2,LL,ALM; 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: g:=LM(Polys[gg[2]]); if LM(Polys[p[2]])*g=LM(Polys[p[1]]) then c1:=c1+1; RETURN(true); elif divide(LM(Polys[p[1]]),lcm(LM(Polys[p[2]]),g)) and LM(Polys[p[1]])<>lcm(LM(Polys[p[2]]),g) then c2:=c2+1; RETURN(true); fi; RETURN(false); end: ##################### ##################### HeadNormalForm:=proc(p,InvDivision) global T,tord,Q,Polys,Vars; local h,G,flag,i,divisor,g; #option trace; h:=Polys[p[1]]; #G:=[seq(Polys[q[1]], q in T)]; G:=[seq(LM(Polys[q[1]]), q in T)]; flag:=false: NN:=NULL; for i from 1 to nops(T) while not flag do if InvDivide2(LM(h),LM(Polys[T[i][1]]),G,Vars,InvDivision) then divisor:=i; flag:=true; fi: od: if not flag then RETURN(Polys[p[1]]); fi; if p[1]<>p[2] then if Criteria2(p,T[divisor],InvDivision) then RETURN(c1c2c3c4); else i:=1; while h<>0 and i<=nops(T) do g:=Polys(T[i][1]); if InvDivide2(LM(h),LM(g),G,Vars,InvDivision) then h:=SPolynomial(h,g,tord); i:=1; else i:=i+1; fi; od; fi: else i:=1; while h<>0 and i<=nops(T) do g:=Polys(T[i][1]); if InvDivide2(LM(h),LM(g),G,Vars,InvDivision) then h:=SPolynomial(h,g,tord); i:=1; else i:=i+1; fi; od; fi; RETURN(h); end: ##################### HeadReduce:=proc(QQ,InvDivision) #option trace; global Vars,T,tord,Polys,rdz,isr,n; local S,Q,p,h,SS,s; S:=QQ: Q:=NULL; while nops(S)<>0 do p:=S[1]; S:=S[2..-1]; h:=HeadNormalForm(p,InvDivision); if h<>0 and h<>c1c2c3c4 then if LM(Polys[p[1]])<>LM(h) then n:=n+1: Polys(n):=h; Q:=Q,[n,n,{}]; else Polys(p[1]):=h; Q:=Q,p; fi; elif h=0 then rdz:=rdz+1: if p[1]=p[2] then SS:=NULL; for s in S do if s[2]<>p[1] then SS:=SS,s; else isr:=isr+1: fi; od: S:=[SS]; fi: fi: od: RETURN([Q]); end: ##################### InvBasis2:=proc(Ideal,Tord,InvDivision) #option trace; global tord,T,Q,Polys,Vars,c1,c2,rdz,ALM,isr,n; local InvDiv,F,p,TT,q,h,NM,np,j,firsttime,firstbytes,secondtime,secondbytes,deg; 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,J); F:=sort(F,(a,b)->TestOrder(a,b,tord)); n:=nops(F): Polys:=Array(1..n,i->F[i]); deg:=0: isr:=0: T:=[[1,1,{}]]; Q:=[seq([i,i,{}],i=2..nops(F))]; Q:=HeadReduce(Q,InvDivision); Q:=sort(Q,tidy); while nops(Q)<>0 do p:=Q[1]; Q:=Q[2..-1]; if p[1]=p[2] then 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]; fi; deg:=max(deg,degree(Polys[p[1]])): h:=AutoReduction(Polys[p[1]],[seq(Polys[q[1]],q in T)],tord,InvDivision); 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],{}]]; od: T[i][3]:=(NM minus q[3]) union (NM intersect q[3]); od: Q:=HeadReduce(Q,InvDivision); Q:=sort(Q,tidy); od: secondtime,secondbytes:=kernelopts(cputime,bytesused); T:=[seq(Polys[tt[1]], tt in T)]; appendto(test); 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,ISRcrit,isr): printf("%-1s %1s %1s : %3a\n",Max,of,degr,deg): printf("%-1s %1s %1s : %3a\n",Num,of,polyno,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 # ############################# 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: #################################"Examples writeto(test); "This is Gerdt test"; writeto(terminal); appendto(test); 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]: InvBasis2(F,tdeg(x,y,z),J): appendto(test); 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]: InvBasis2(F,tdeg(xxx,yyy,zzz,ttt,uuu,vvv,www,aaa),J): appendto(test); print("####################Example 3"); appendto(terminal); F:=[a*y^2+b*x^3+c, 2*a*y, 3*b*x^2]: InvBasis2(F,tdeg(x,y,a,b,c),J): appendto(test); 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]: InvBasis2(F,tdeg(x,y,a,b,c,e,f),J): appendto(test); 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]: InvBasis2(F,tdeg(x,y,a,b,c,d),J): appendto(test); 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]: InvBasis2(F,tdeg(x,y,a,b,c),J): appendto(test); 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]: InvBasis2(F,tdeg(x,y,a,b),J): appendto(test); 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]: InvBasis2(F,tdeg(z, y, x, t0, a, h),J): appendto(test); 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]: InvBasis2(F,tdeg(a,b,c,d,e),J): appendto(test); 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]: InvBasis2(F, tdeg(x1, x2, x3, x4, h), J): appendto(test); 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]: InvBasis2(F, tdeg(x,y,z,h), J): appendto(test); 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]: InvBasis2(F, tdeg(ax, ay, az, at, au, av), J): appendto(test); 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]: InvBasis2(F, tdeg(x1, x2, x3, x4,x5,x6,x7), J): appendto(test); 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]: InvBasis2(F, tdeg(bb, xx, yy, zz, ss, tt, uu, vv, ww), J): appendto(test); 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]: InvBasis2(F,tdeg(a,b,c,d,e,f),J): appendto(test); 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]: InvBasis2(F, tdeg(x,y,t), J): appendto(test); 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 ]: InvBasis2(F, tdeg(x1, x2, x3, x4,x5,x6,x7), J):