#### This proram is the main program of Staggered-linear basis ####` with(Groebner): with(PolynomialIdeals): lm := proc (f, tord) options operator, arrow; LeadingMonomial(f, tord) end proc: lc := proc (f, tord) options operator, arrow; LeadingCoefficient(f, tord) end proc: lt := proc (f, tord) options operator, arrow; LeadingTerm(f, tord)[1]*LeadingTerm(f,tord)[2] end proc: LM := proc(u) local i; for i to nops(u) do if u[i] <> 0 then RETURN([u[i], i]) end if: end do: RETURN(0) end proc: #`#########################` pot := proc (a, b, tord) i := a[2]; j := b[2]; t1 := a[1]; t2 := b[1]; if i < j then RETURN(true): elif j < i then RETURN(false): elif TestOrder(t1, t2, tord) then RETURN(true): else RETURN(false): end if: end proc: #`####################` Pot := proc (a, b, tord) #option trace: global Fsub: t1 := a[1]; i := a[2]; t2 := b[1]; j := b[2]; if Fsub[j][3][2] < Fsub[i][3][2] then RETURN(false): elif Fsub[i][3][2] < Fsub[j][3][2] then RETURN(true): elif TestOrder(t1*Fsub[i][3][1], t2*Fsub[j][3][1], tord) then RETURN(true): else RETURN(false): end if: end proc: #`########## Cover function ############` Cover := proc (p1, p2, tord) local u1, u2, v1, v2; u1 := p1[1]; u2 := p2[1]; if not evalb(u1[2] = u2[2]) then RETURN(false): elif divide(u1[1], u2[1],'q') then t := q*p2[2]; if TestOrder(t, p1[2],tord) and t <> p1[2] then RETURN(true): end if: end if; RETURN(false): end proc: #`########## Regtopred function #########` regtopred := proc (p1,tord) #option trace: local Flag; global Fsub: Flag := true; G:=Fsub: n := op(2, op(2, G)); #Q:=[seq(0,i=1..n)]: Q:=Array(1..n,i->0): r:=[0,0]: p := [p1[1], p1[2]]; while p[2]<> 0 do i := 1; flag := false; u1, v1:= p[1], p[2]; while flag = false and i <= n do u2, v2 := G[i][3], G[i][1]; if divide(lm(v1, tord), G[i][5], 'q') then t := q; c := lc(v1, tord)/G[i][6]; if pot([t*u2[1], u2[2]], u1, tord) then if [t*u2[1], u2[2]]<>u1 or c <> 1 then p := [u1, v1-expand(c*t*v2)]; Q[i]:=Q[i]+c*t: flag := true; if [t*u2[1], u2[2]] = u1 then p[2] := p[2]/(1-c): end if: end if: end if: end if; i := i+1: end do; if flag = false then r:=[u1,simplify(r[2]+lt(p[2],tord))]: p[2]:=simplify(p[2]-lt(p[2],tord)): end if: end do; RETURN(r,Q): end proc: #########Checkpair############## checkpair:=proc(m,A) #option trace: flag:=true: i:=1: while flag and i<= nops(A) do if divide(m,A[i],'q')=true then RETURN(false): fi: i:=i+1: od: RETURN(true); end proc: ######### Staggered linear basis algorithm based on the incrimental process###### Stag:=proc(F1,tord) #option trace: global Fsub: tim1,b1:=kernelopts(cputime,bytesused): NumOfTop:=0: NumOfZero:=0: NumOfCov:=0: NumOfStaggeredRem:=0: csp:=0: F:=[seq([F1[i],{},[1,i],i,lm(F1[i],tord),lc(F1[i],tord)],i=1..nops(F1))]: Fsub:=[F[1]]: Fsub:=Array(1..1,Fsub): F:=F[2..-1]: H:=[]: JP:=Array([]): while op(F)<>NULL do n:=op(2, op(2, Fsub)): f:=F[1]: F:=F[2..-1]: f[2]:={seq(g[5],g in Fsub)}: f[4]:=n+1: LMG:=[seq(Fsub[i][5],i=1..n)]: #print("f",f): #print("Fsub",Fsub): for p in Fsub do Lcm:=lcm(p[5],f[5]): m1:=(Lcm)/(p[5]): m2:=(Lcm)/(f[5]): m1u1:=[m1*p[3][1],p[3][2]]: m2u2:=[m2*f[3][1],f[3][2]]: if m1u1<>m2u2 then Flag:= pot(m1u1,m2u2,tord): Num:=op(2,op(2,JP)): if Flag and checkpair(m2,f[2]) then JP:=Array(1..Num+1,JP): JP[Num+1]:=[m2,n+1]: elif Flag=false and checkpair(m1,p[2]) then JP:=Array(1..Num+1,JP): JP[Num+1]:=[m1,p[4]]: #elif (Flag and checkpair(m2,f[2])=false) or (Flag=false and checkpair(m1,p[2])=false) then #NumOfCov:=NumOfCov+1: else #prnit("Removed by staggered"); NumOfStaggeredRem:=NumOfStaggeredRem+1: fi: elif checkpair(m2,f[2]) then csp:=csp+1: p[2]:=p[2] union {m1}: fi: od: Fsub:=Array(1..n+1,Fsub): Fsub[n+1]:=f: n:=n+1: ##In this loop we try to obtain a Groebner basis for the ideal while ArrayNumElems(JP)<>0 do JP:=sort(JP,(a,b)->Pot(a,b,tord)): t,i:=JP[1][1],JP[1][2]: #print("####################################################################"); #print("t,i", t,i); p:=Fsub[i]: #print("p",p); JP:=Array(1..ArrayNumElems(JP)-1,i->JP[i+1]): jp:=[[t*p[3][1],p[3][2]],expand(t*p[1]),t*p[5]]: flag:=false: if checkpair(t,p[2])=false then flag:=true: NumOfStaggeredRem:=NumOfStaggeredRem+1: #print("Remove by staggered"); ###### elif checkpair`(t,LMG)=false then ###### flag:=true: #####NumOfCov:=NumOfCov+1: #print("Remove by forbidden LM"); fi: if flag=false then r,Q:=regtopred(jp,tord): #print("r",r): #print("Q",Q): NumOfTop:=NumOfTop+1: if r[2]=0 then NumOfZero:=NumOfZero+1: Fsub[i][2]:=Fsub[i][2] union {t}: else qoutiontt:={}: for su in Fsub[i][2] do l:=numer(simplify(su/t)): qoutiontt:={op(qoutiontt),l}: end do: r:=[r[2],qoutiontt,r[1],n+1,lm(r[2],tord),lc(r[2],tord)]: Fsub[i][2]:=Fsub[i][2] union {t}: #print("r",r): JPnew:=Array([]): for p in Fsub do Lcm:=lcm(p[5],r[5]): m1:=Lcm/p[5]: m2:=Lcm/r[5]: t1u1:=[m1*p[3][1],p[3][2]]: t2u2:=[m2*r[3][1],r[3][2]]: if t1u1<>t2u2 then Flag:= pot(t1u1,t2u2,tord): if Flag and checkpair(m2,r[2]) then endd:=op(2,op(2,JPnew)): JPnew:=Array(1..endd+1,JPnew): JPnew[endd+1]:=[m2,n+1]: elif Flag=false and checkpair(m1,p[2]) then endd:=op(2,op(2,JPnew)): JPnew:=Array(1..endd+1,JPnew): JPnew[endd+1]:=[m1,p[4]]: #elif(Flag and checkpair(m2,r[2])=false) or (Flag=false and checkpair(m1,p[4])=false) then #` #NumOfCov:=NumOfCov+1: else NumOfStaggeredRem:=NumOfStaggeredRem+1: fi: elif checkpair(m2,r[2]) then csp:=csp+1: p[2]:=p[2] union {m1}: fi: od: n:=n+1: Fsub:=Array(1..n,Fsub): Fsub[n]:=r: #print("Fsub",Fsub): #print("JPnew",JPnew,JP): endJp:=ArrayNumElems(JP); nnn:=ArrayNumElems(JPnew): JP:=Array(1..endJp+nnn,JP): for k from endJp+1 to endJp+nnn do JP[k]:=JPnew[k-endJp]; od: #JP:=[op(JP),op(JPnew)]: JP:=sort(JP,(a,b)->Pot(a,b,tord)): #print("JP",JP): fi: fi: od: #print("JP is NUll"): end do: tim2,b2:=kernelopts(cputime,bytesused): L:=[seq(f[1],f in Fsub)]: EB:=[seq(Fsub[i][1..2],i=1..ArrayNumElems(Fsub))]; lprint("The required time is", tim2-tim1); lprint("The used memory", evalf((b2-b1)/10^6)); lprint("The number of reductions to zero", NumOfZero); lprint("The number of removed pairs by cover criterion", NumOfCov); lprint("The number of removed pairs by csp criterion", csp+NumOfStaggeredRem); lprint("The number of Groebner basis elements", nops(L)); #lprint("Is Grobner?", IsGrobner(F1,L,tord)); #lprint("Is Staggrred?", staggered(EB,tord)); RETURN(); end proc: ############################## IsGrobner:=proc(A,H,T) #option trace; 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: ############################## staggered:=proc(F,tord) #option trace: local count,n,i,j,Lcm,m1,m2,T,t,flag,l; n:=nops(F): count:=0: for i from 1 to n do for j from i+1 to n do count:=count+1: Lcm:=lcm(lm(F[i][1],tord),lm(F[j][1],tord)): m1:=(Lcm)/(lm(F[i][1],tord)): m2:=(Lcm)/(lm(F[j][1],tord)): if not (m1 in LeadingMonomial(, tord)) and not (m2 in LeadingMonomial(, tord)) then print("B is not an independent set and not Staggered basis",count); RETURN(); fi: od: od: #print("B is independent set",count); for i from 1 to n do T:=F[i][2]: while op(T)<> NULL do t:=T[1]: T:=T minus {t}: flag:=false: l:=1: while flag=false and l<= n do if divide(t*lm(F[i][1],tord), lm(F[l][1],tord), 'q') then if (q in LeadingMonomial(, tord)) then flag:=true: fi: fi: l:=l+1: od: if flag=false then print("B is not a generating set and not Staggered"); RETURN(0): fi: od: od: #print("B is a generating set"); RETURN(true): end proc: ############################# ############################# F:=[a*y^2+b*x^3+c, 2*a*y, 3*b*x^2]; Stag(F,tdeg(x,y,a,b,c)); F:=[x*w+w, y^2*z+x,x*y+z*w,x*y^2+z*w^2]; Stag(F,tdeg(x,y,z,w)); F := [b*y^2+2*c*x*y+2*d*x+2*e*y+x^2+f, c*y+d+x, b*y+c*x+e]; Stag(F, tdeg(x, y, a, b, c, d, e, f)); 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]; Stag(F,tdeg(x,y,a,b,c)); #Basis(F,tdeg(x,y,a,b,c)); #F := [x^6+a*y^3-y, y^6+b*x^3-x, 36*x^5*y^5-9*a*b*x^2*y^2+3*a*y^2+3*b*x^2-1]; #Stag(F, tdeg(x, y, a, b)); F := [a*h-h*x-t0*y+y*z, a*h-h*y+t0*z-x*z, a*h-h*z+t0*x-t0*y, a*h-h*t0+x*y-x*z]; Stag(F, tdeg(z, y, x, t0, a, h)); F := [x[2]*x[3]^2+x[1]^2, x[3]*x[5]+x[2], x[2]*x[5]+x[3]*x[4], x[2]*x[3]+x[4], x[1]*x[2]+x[5]^2]; Stag(F, tdeg(x[1], x[2], x[3], x[4], x[5])); F := [x1^2*x4+x2*x3^2, x3*x5+x2, x2*x5+x3*x4, x2*x3+x4^2, x1*x2+x5^2]; Stag(F, tdeg(x1, x2, x3, x4, x5)); F := [x1^2*x4+x2*x3^2, x2^2+x3*x5, x2*x5+x3*x4, x2*x3+x4^2, x1*x2+x5^2]; Stag(F, tdeg(x1, x2, x3, x4, x5)); F := [x1^2*x4+x2*x3^2, x2^2+x3*x5, x1*x2*x5+x3*x4^2, x1*x5+x2*x3+x4^2, x1*x2+x5^2]; Stag(F, tdeg(x1, x2, x3, x4, x5)); 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]; Stag(F,tdeg(xxx,yyy,zzz,ttt,uuu,vvv,www,aaa)); F := [10*h^3-11*h^2*x4+10*x1^2*x4+10*x2^2*x4+10*x3^2*x4, 10*h^3-11*h^2*x3+10*x1^2*x3+10*x2^2*x3+10*x3*x4^2, 10*h^3-11*h^2*x1+10*x1*x2^2+10*x1*x3^2+10*x1*x4^2, 10*h^3-11*h^2*x2+10*x1^2*x2+10*x2*x3^2+10*x2*x4^2]; Stag(F, tdeg(x1, x2, x3, x4, h)); F := [h^2*x^2-2*h^2*x*y+h^2*y^2+h^2*z^2+x*y^2*z+y^4, -3*h^5-2*h^2*x^2*y+x*y^4+y*z^4, -2*h^3*x*y+h*x*y^2*z+h*y^4-x^3*y^2+x*y*z^3]; Stag(F, tdeg(x, y, z, h)); 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*ay*at+2*az*au+au^2+2*at*av-at, 2*ax*at+2*ay*au+2*at*au+2*az*av-az]; Stag(F, tdeg(ax, ay, az, at, au, av)): 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]; Stag(F, tdeg(a, b, c, d, e)): #######Sturmfles F :=[bb*vv+ss*uu, bb*ww+tt*uu, ss*ww+tt*vv, bb*yy+ss*xx, bb*zzz+tt*xx, ss*zzz+tt*yy, uu*yy+vv*xx, uu*zzz+ww*xx, vv*zzz+ww*yy]; Stag(F, tdeg(bb, xx, yy, zzz, ss, tt, uu, vv, ww));