#`### 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)]:
 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######
 F5:=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)):
            JP:=Array(1..Num+1,JP):
            if Flag and checkpair(m2,f[2]) then
               JP[Num+1]:=[m2,n+1]:
            elif Flag=false and checkpair(m1,p[2]) then
               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:
     #print("Fsub",Fsub):
     JP:=sort(JP,(a,b)->Pot(a,b,tord)):
     #print("JP",JP):
     ##In this loop we try to obtain a Groebner basis for the ideal <Fsub>
     while op(JP)<>NULL do
        t,i:=JP[1][1],JP[1][2]:
        #print("####################################################################");
        #print("t,i", t,i);
        p:=Fsub[i]:
        #print("p",p);
        init:=op(1,op(2,JP)):
        JP:=JP[init+1..-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):
                      endd:=op(2,op(2,JPnew)):
                      if Flag and checkpair(m2,r[2]) then
                          JPnew:=Array(1..endd+1):
                          JPnew[endd+1]:=[m2,n+1]:
                      elif Flag=false and checkpair(m1,p[2]) then
                          JPnew[endd]:=[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):
               endJp:=op(2,op(2,JP)):
               nnn:=nops(JPnew):
               JP:=Array(1..endJp+nnn,JP):
               for k from endJp+1 to endJp+nnn
                    JP[k]:=JPnes[1]
               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)]:
printf("%-1s %1s %1s %1s         : %3a %3a\\n",The, cpu, time, is,tim2-tim1,(sec)):  
printf("%-1s %1s %1s         : %5a%3a\\n",The,used,memory,b2-b1,(bytes)):  
printf("%-1s %1s %1s %1s       : %5a\\n",Number, of, zero, reduction,NumOfZero):  
printf("%-1s %1s %1s:           %5g\\n",Num,of,CoverCriteria,NumOfCov):
printf("csp", csp); 
#print("IsGrobner", IsGrobner(F1, L, tord)); 
#staggered(Fsub, tord); 
#RETURN(csp)
end proc: