################ This program computes a minimal H-bases for both homegenous and non homogenoeus ideals

with(Groebner):

LT := proc (f, tord) options operator, arrow; LeadingMonomial(f, tord)*LeadingCoefficient(f, tord) end proc:

LM := proc (f, tord) options operator, arrow; LeadingMonomial(f, tord) end proc:

################## This is a divison algorithm #####################
Div := proc (f, L, T)
#option trace;
local p, r, m, q, flag, i;
p := f;
r := 0;
m := nops(L);
q := [seq(0, i = 1 .. m)];
while p <> 0 do
   i := 1;
   flag := false;
   while not flag and i <= m do
      if divide(LT(p,T),LT(L[i][1],T),'t') then
          q[i]:= simplify(q[i]+t);
          p := simplify(p-t*L[i][1]);
          flag := true;
      else
          i := i+1;
      end if;
   end do;
   if not flag then
      r := simplify(r+LT(p,T));
      p := simplify(p-LT(p,T));
   end if;
end do;
RETURN(r, q);
end proc:

################## Rctified Minimal H-bases ###########################

YilmazBase:=proc(F,tord)
#option trace;
local G,T,k,t,H,H1,B,p,i,j,c,a,s,q,r,d,P,R,n,K,L,m,re,flag,rem,l,relation,Ly,Yr,Yo,ls,t1,countReduc;
firsttime, firstbytes := kernelopts(cputime, bytesused);
L:=NULL:
MinB:=NULL:
syz:=NULL:
rem:=[]:
countReduc:=0:
F1:=[F[1]]:
for i from 1 to nops(F) do
  if member(F[i],F1,'w')=false then
     F1:=[op(F1),F[i]]:
  end if:
end do:
G:=[seq([F1[i],i],i=1..nops(F1))]:
T:=NULL:
H:=NULL:
Ly:=NULL:
Yr:=[seq(Y[i],i=1..nops(F1))]:
B:=[seq(seq([G[i],G[j]],j=i+1..nops(G)),i=1..nops(G))]:
while B<>[] do
    relation:=0:
    p:=B[1]:
    B:=B[2..-1]:
    i:=p[1][2]:
    j:=p[2][2]:
    c[i][j]:=lcm(LM(p[1][1],tord),LM(p[2][1],tord)):
    a[i][j]:=(c[i][j])/(LT(p[1][1],tord)):
    a[j][i]:=(c[i][j])/(LT(p[2][1],tord)):
    s[i][j]:=simplify(a[i][j]*p[1][1]-a[j][i]*p[2][1]):
    q:='q';
    r,q:=Div(s[i][j],G,tord);
    countReduc:=countReduc+1:
    q:=-q;
    q[i]:=q[i]+a[i][j];
    q[j]:=q[j]-a[j][i];
    if r<>0 then
        ind:=nops(Yr);
        Yr:=[op(Yr),Y[ind+1]]:
        P:=seq([f,[r,ind+1]],f=G):
        B:=[op(B),P]:
        G:=[op(G),[r,ind+1]]:
        q:=[op(q),-1]:
    end if;
        d:=max(seq(degree(G[l][1]*q[l]),l=1..nops(G)));
        for n from 1 to nops(q) do
                if type(q[n],constant) and q[n]<>0 and degree(G[n][1])=d then
                     relation:=relation+q[n]*Y[n];
                end if:
        end do;
        if relation<>0 then
           Ly:=Ly,relation:
        end if:
        H:=H,q:
end do:
H:=[H]:
Ly:=[Ly]:
for h in H do
  nG:=nops(G)-nops(h):
  h:=[op(h),seq(0,i=1..nG)]:
  syz:=syz,h:
end do:
syz:=[syz]:
T:=Basis(Ly,plex(op(Yr))):
rem:=[seq(LM(f,plex(op(Yr))),f=T)]:
R:=[op({op(Yr)} minus {op(rem)})]:
for i from 1 to nops(R) do
    member(R[i],Yr,'w');
    t2:=w:
    MinB:=MinB,G[t2][1]:
end do:
secondtime, secondbytes := kernelopts(cputime, bytesused);
time:=secondtime-firsttime;
memory:=secondbytes-firstbytes;
print("time",time);
print("memory",memory);
print("The number of reductions",countReduc);
RETURN(R,nops(G));
end proc:


##################################################################
##################################################################
################ This is a program that computes a minimal H-basis with applying update algorithm to remove redundant S-polies satisfying first-second Buchberger's criteria 


###This program is the non-hom version of the Update alg + Yilmaz and it works correctly for hom and non-hom examples

with(Groebner):

LT := proc (f, tord) options operator, arrow; LeadingMonomial(f, tord)*LeadingCoefficient(f, tord) end proc:

LM := proc (f, tord) options operator, arrow; LeadingMonomial(f, tord) end proc:

########### Division algorithm ########################

Div := proc (f, L, T)
#option trace;
local p, r, m, q, flag, i;
p := f;
r := 0;
m := nops(L);
q := [seq(0, i = 1 .. m)];
while p <> 0 do
   i := 1;
   flag := false;
   while not flag and i <= m do
      if divide(LT(p,T), LT(L[i][1],T), 't') then
          q[i] := simplify(q[i]+t);
          p := simplify(p-t*L[i][1]);
          flag := true;
      else
          i := i+1;
      end if;
   end do;
   if not flag then
      r := simplify(r+LT(p,T));
      p := simplify(p-LT(p,T));
   end if;
end do;
RETURN(r, q);
end proc:

########################## Update algorithm ####################
UPDATE:=proc(G0,B0,h,tord)
#option trace;
D1:=[]:
C:=[seq([h,g],g=G0)]:
while C<>[] do
  p:=C[1]:
  C:=C[2..-1]:
  flag1:=true:
  for c in C do
   if divide(lcm(LM(p[1][1],tord),LM(p[2][1],tord)),lcm(LM(c[1][1],tord),LM(c[2][1],tord)))=true then
      flag1:=false:
      break:
   end if:
  end do:
  flag2:=true:
  if flag1 then
    if D1<>[] then
       for d in D1 do
          if divide(lcm(LM(p[1][1],tord),LM(p[2][1],tord)),lcm(LM(d[1][1],tord),LM(d[2][1],tord)))=true then
              flag2:=false;
              break;
          end if:
       end do:
    end if:
  end if:
  flag3:=false:
  if flag1 and flag2 then
    flag3:=true;
  end if:
  if gcd(LM(p[1][1],tord),LM(p[2][1],tord))=1 or flag3=true then
    D1:=[op(D1),p]:
  end if:
end do:
E:=[]:
while D1<>[] do
   p:=D1[1]:
   D1:=D1[2..-1]:
   if gcd(LM(p[1][1],tord),LM(p[2][1],tord))<>1 then
      E:=[op(E),p]:
   end if:
end do:
B1:=[]:
BOld:=B0:
while BOld<>[] do
   p:=BOld[1]:
   BOld:=BOld[2..-1]:
   if divide(lcm(LM(p[1][1],tord),LM(p[2][1],tord)),LM(h[1],tord))=false or lcm(LM(p[1][1],tord),LM(h[1],tord))=lcm(LM(p[1][1],tord),LM(p[2][1],tord)) or
   lcm(LM(h[1],tord),LM(p[2][1],tord))=lcm(LM(p[1][1],tord),LM(p[2][1],tord)) then
       B1:=[op(B1),p]:
   end if:
end do:
B1:=[op(B1),op(E)]:
Gold:=G0:
Gnew:=NULL:
while Gold<>[] do
   g:=Gold[1]:
   Gold:=Gold[2..-1]:
   if divide(g[1][1],h[1][1])=false then
      Gnew:=Gnew,g:
   end if:
end do:
Gnew:=[Gnew,h]:
RETURN(B1,Gnew);
end proc:

##########This function reduces list L from H#######
minuslists:=proc(H,L)
local x,N,plc:
    N:=H:
    for x in L do
       if x in N then
            plc:='plc':
        member(x,N,'plc');
        N:=subsop(plc=NULL,N):
    fi:
    od:
RETURN(N):
end:

#############This function sorts a list of critical pairs w.r.t Normal strategy################
SOrt:=proc(B,tord)
#option trace:
B0:=B:
BOrder:=NULL:
l:=[seq(lcm(LM(B0[i][1][1],tord),LM(B0[i][2][1],tord)),i=1..nops(B0))]:
deg:=[seq(degree(a),a=l)]:
while l<>[] do
  minim:=min(deg);
  member(minim,deg,'w');
  i:=w:
  BOrder:=BOrder,B0[i];
  l:=minuslists(l,[l[i]]):
  deg:=minuslists(deg,[deg[i]]):
  B0:=minuslists(B0,[B0[i]]):
end do:
RETURN([BOrder]);
end proc:


###############################Minimal H-basis for polynomial ideals ###################################
YilmazBase:=proc(FF,tord)
#option trace;
local G,T,k,t,H,H1,B,p,i,j,c,a,s,q,r,d,P,R,n,K,L,m,re,flag,flag3,flag4,pair,rem,v,l,l2,D1,relation,Ly,Yr,Yo,ls,t1,lt;
firsttime, firstbytes := kernelopts(cputime, bytesused);
L:=NULL:
MinB:=NULL:
syz:=NULL:
rem:=[]:
count:=0:
F1:=[FF[1]]:
for i from 1 to nops(FF) do
  if member(FF[i],F1,'w')=false then
     F1:=[op(F1),FF[i]]:
  end if:
end do:
F2:=subs(0=NULL,InterReduce(F1,tord)):
F:=[seq([F2[i],i],i=1..nops(F2))]:
G:=[]:
G2:=F2:
T:=NULL:
H:=NULL:
Ly:=NULL:
Yr:=[seq(Y[i],i=1..nops(F))]:
B:=[]:
step:=1:
while F<>[] do
    i:=F[1]:
    F:=F[2..-1]:
    B,G:=UPDATE(G,B,i,tord):
end do:
B:=SOrt(B,tord):
while B<>[] do
    relation:=0:
    p:=B[1]:
    B:=B[2..-1]:
    i:=p[1][2]:
    j:=p[2][2]:
    c[i][j]:=lcm(LM(p[1][1],tord),LM(p[2][1],tord)):
    a[i][j]:=(c[i][j])/(LT(p[1][1],tord)):
    a[j][i]:=(c[i][j])/(LT(p[2][1],tord)):
    s[i][j]:=simplify(a[i][j]*p[1][1]-a[j][i]*p[2][1]):
    q:='q';
    r,q:=Div(s[i][j],G,tord);
    q:=-q;
    q[i]:=q[i]+a[i][j];
    q[j]:=q[j]-a[j][i];
    count:=count+1:
    if r<>0 then
        ind:=nops(Yr);
        Yr:=[op(Yr),Y[ind+1]]:
        G2:=[op(G2),r]:
        q:=[op(q),-1]:
        B,G:=UPDATE(G,B,[r,ind+1],tord):
        B:=SOrt(B,tord):
    end if;
    d:=max(seq(degree(G2[l]*q[l]),l=1..nops(G2)));
    for n from 1 to nops(q) do
        if type(q[n],constant) and q[n]<>0 and degree(G2[n])=d then
            relation:=relation+q[n]*Y[n];
        end if:
    end do;
        if relation<>0 then
           Ly:=Ly,relation:
        end if:
        H:=H,q:
end do:
H:=[H]:
Ly:=[Ly]:
for h in H do
  nG:=nops(G2)-nops(h):
  h:=[op(h),seq(0,i=1..nG)]:
  syz:=syz,h:
end do:
syz:=[syz]:
T:=Basis(Ly,plex(op(Yr))):
rem:=[seq(LM(f,plex(op(Yr))),f=T)]:
R:=[op({op(Yr)} minus {op(rem)})]:
for i from 1 to nops(R) do
    member(R[i],Yr,'w');
    t2:=w:
    MinB:=MinB,G2[t2]:
end do:
secondtime, secondbytes := kernelopts(cputime, bytesused);
time:=secondtime-firsttime;
memory:=secondbytes-firstbytes;
print("time",time);
print("Memory",memory);
print("Number of reductions",count);
RETURN(R,nops(G2));
end proc:

##########################################################################################
##########################################################################################

###### This is a GVW algorithm + MinimalH-bases algorithm which we removed the cover criterion for Jpairs and Trivial syzygies ####

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:

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:

############ The order here is Pot in the sence that if j<i then x^a ei< x^b ej #############

pot := proc (a, b, tord)
i := a[2];
j := b[2];
t1 := a[1];
t2 := b[1];
if j < i then
   RETURN(true):
elif i < j then
   RETURN(false):
elif TestOrder(t1, t2, tord) then
   RETURN(true):
else RETURN(false):
end if:
end proc:

##################### This function orders Jpairs based on pot described above ###################
Pot := proc (a, b, tord)
#option trace:
global G:
t1 := a[1];
i := a[2];
t2 := b[1];
j := b[2];
if G[j][5][2] < G[i][5][2] then
   RETURN(true):
elif G[i][5][2] < G[j][5][2] then
   RETURN(false):
elif TestOrder(t1*G[i][5][2], t2*G[j][5][2], tord) then
   RETURN(true):
else RETURN(false):
end if:
end proc:

########### Cover function ############
Cover := proc (p1, p2, tord)
#option trace:
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:
global G:
local Flag;
Flag := true;
n := op(2, op(2, G));
p := [p1[1], p1[2]];
Q := Array(1..n,0);
while Flag do
   i := 1;
   flag := false;
   u1, v1:= p[1], p[2];
   while flag = false and i <= n do
      u2, v2 := G[i][5], G[i][1];
      if divide(lm(v1, tord), G[i][2], 'q') then
         t := q;
         c := lc(v1, tord)/G[i][3];
         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
      Flag := false:
   end if:
end do;
RETURN(p,Q):
end proc:


############## GVW Algorithm based on Gao et.al paper ######################
F5 := proc (F, tord)
#option trace;
global G:
tim1,b1:=kernelopts(cputime,bytesused):
Ly:=NULL:
NumOfTop:=0:
NumOfZero:=0:
NumOfCov:=0:
n := nops(F);
JP := [];
G:=[seq([F[i],lm(F[i],tord),lc(F[i],tord),i,[1,i]],i=1..nops(F))]:
print("G",G);

G:=Array(1..n,G):
H := [seq(seq([G[j][2], i], j = i+1 .. n), i = 1 .. n)];
Yr:=[seq(Y[i],i=1..n)]:


#print("making Jpairs in first step");
for i to n do
   for j from i+1 to n do
      Lcm := lcm(G[i][2], G[j][2]);
      t[i][j] := Lcm/G[i][2];
      t[j][i] := Lcm/G[j][2];
      t1T := [t[i][j]*G[i][5][1], G[i][5][2]];
      t2T := [t[j][i]*G[j][5][1], G[j][5][2]];
      if t1T<>t2T or G[i][3] <> G[j][3] then
         if pot(t1T, t2T, tord) then
            JP := [op(JP), [t[j][i], j]]:
            print("######i,j######",i,j);
         else
            JP := [op(JP), [t[i][j], i]]:
            print("######i,j######",i,j);
         end if:
      else
            print("##### Jpair{",i,j,"}is not constructed #####");
      end if:
   end do:
end do;

#print("Cheking whether Jpairs will be covered by the elements of G or not");
for jp in JP do
    m, f := jp[1], G[jp[2]];
    mf := [[m*f[5][1],f[5][2]],m*f[2]];
    flag := false;

    #to being covered by G
    l:=1:
    while flag =false and l <= n do
         g:=G[l]:
         gs:=[g[5],g[2]]:
         if Cover(mf, gs, tord) then
             plc := 'plc';
             member(jp, JP, 'plc');
             JP := subsop(plc = NULL, JP);
             NumOfCov:=NumOfCov+1:
             flag := true;
             print("t,i",m,f);

         else
             l:=l+1:
         fi:
    end do;
end do;

#print("sorting Jpairs based on Pot ordering");
JP := sort(JP,(a,b)->Pot(a,b,tord)):


#print("treating the Jpairs" );
while op(JP) <> NULL do
    relation:=0:
    t,i:=JP[1][1],JP[1][2]:
    print("############# t,i ############",t,i);
    p:=G[i]:
    JP:=JP[2..-1]:
    jp:=[[t*p[5][1],p[5][2]],expand(t*p[1]),t*p[2]]:


    #print("starting Cover criterion");
    flag:=false:

    #to being covered by G
    l:=1:
    Rrelation:=0:
    while flag=false and l<= n do
        g:= G[l]:
        sg:=[g[5],g[2]]:
        if Cover([jp[1],jp[3]],sg,tord) then
            NumOfCov:=NumOfCov+1:
            flag:=true:
            print("tp is covered by gl",l);
            r,Q :=regtopred(jp,tord):
            Q[i]:=Q[i]-t;
            print("r,Q",r,Q);
        else
            l:=l+1:
        fi:
    end do:
    if flag=true then
       nG:=op(2, op(2, G));
       if r[2]<> 0 then
           nG:=nG+1:
           Q :=Array(1..nG,Q):
           Q[nG]:=1:
       fi:
       d:=max(seq(degree(expand(G[l][1]*Q[l])),l=1..nG-1));
       if r[2]<>0 then
          d:=max(d,degree(r[2]));
       fi:
       print("d",d);
       for i from 1 to nG-1 do
           if type(Q[i],constant) and Q[i]<>0  and degree(G[i][1])=d then
                     Rrelation:=Rrelation+Q[i]*Y[i];
           fi:
       end do:
       if r[2]<>0 and degree(r[2])=d and Rrelation<>0 then
          Rrelation:=Rrelation+Y[nG]:
       fi:
       if Rrelation<>0 then
           print("Rrelation",Rrelation);
       end if:  
    fi:

   if  flag=false then
       r,Q :=regtopred(jp,tord):

       Q[i]:=Q[i]-t;
       nq:=op(2, op(2, G));

       NumOfTop:=NumOfTop+1:
       if r[2]=0 then
           NumOfZero:=NumOfZero+1:
           H:=[op(H),r[1]]:
           flag1:=false:

       else
           r:=[r[2],lm(r[2],tord),lc(r[2],tord),n+1,r[1]]:
           print("r",r);
           ind:=nops(Yr);
           Yr:=[op(Yr),Y[ind+1]]:
           print("Yr",Yr);
           Q :=Array(1..nq+1,Q):
           Q[nq+1]:= 1:
           #print("making trivial Syz and adding them to H");
           for p in G do
               T1:=[r[2]*p[5][1],p[5][2]]:
               T2:=[p[2]*r[5][1],r[5][2]]:
               if T1<>T2 or r[3]<>p[3] then
                   if pot(T1,T2,tord) then
                       H:=[op(H),T2]:
                   else
                       H:=[op(H),T1]:
                   fi:
               fi:
           od:

           JPnew:=[]:

           #print("Making new Jpairs");
           for p in G do
               Lcm:=lcm(p[2],r[2]):
               m1:=Lcm/p[2]:
               m2:=Lcm/r[2]:
               t1u1:=[m1*p[5][1],p[5][2]]:
               t2u2:=[m2*r[5][1],r[5][2]]:
               if t1u1<>t2u2 or p[3]<>r[3] then
                   if pot(t1u1,t2u2,tord) then
                        JPnew:=[op(JPnew),[m2,n+1]]:
                   else
                        JPnew:=[op(JPnew),[m1,p[4]]]:
                   fi:
               fi:
           od:

           n:=n+1:
           G:=Array(1..n,G):
           G[n]:=r:

           #print("Cheking if new Jpairs are covered by G");
           for jpn in JPnew do
               t0,i0:=jpn[1],jpn[2]:
               f0:=G[i0]:
               jpairr:=[[t0*f0[5][1],f0[5][2]],t0*f0[2]]:
               flag1 := false;

               #Cover by G
               l:=1:
               while flag1=false and l<=n-1 do
                    p:=G[l]:
                    sp:=[p[5],p[2]]:
                    if Cover(jpairr, sp, tord) then
                       plc := 'plc';
                       member(jpn, JPnew, 'plc');
                       JPnew := subsop(plc = NULL, JPnew);
                       NumOfCov:=NumOfCov+1:
                       flag1 := true;
                    else
                       l:=l+1:
                    fi:
               end do;

           end do:


           #print("Sorting old Jpairs + new Jpairs");
           JP:=[op(JP),op(JPnew)]:
           JP:=sort(JP,(a,b)->Pot(a,b,tord)):
       fi:
       nG:=op(2, op(2, G));
       GL:=copy(G):
       print("G",convert(GL,list));
       QL:=copy(Q):
       print("Q",convert(QL,list));
       d:=max(seq(degree(expand(G[l][1]*Q[l])),l=1..nG));
       print("d",d);
       for i from 1 to nG do
           if type(Q[i],constant) and Q[i]<>0  and degree(G[i][1])=d then
                     relation:=relation+Q[i]*Y[i];
           end if:
       end do:
       if relation<>0 then
           print("relation",relation);
           Ly:=Ly,relation:
       end if:
   fi:
end do:
Ly:=[Ly]:
T:=Basis(Ly,plex(op(Yr))):
rem:=[seq(LeadingMonomial(f,plex(op(Yr))),f=T)]:
R:=[op({op(Yr)} minus {op(rem)})]:
tim2,b2:=kernelopts(cputime,bytesused):
L:=[seq(g[1],g=G)]:
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("%-1s %1s %1s %1s       : %5a\n",Num, of, TopReduction, is,NumOfTop):
printf("%-1s %1s %1s             : %5a\n",Num,of,poly,nops(L)):
print("IsGrobner",IsGrobner(F,L,tord));
print("Minimal basis",R);
RETURN(L):
end proc:

###########This function cheks that the output of the algorithm is a Grobner basis or not###############
IsGrobner := proc (A, H, T)
local s, j, S, p, F, L, LL;
#option trace;
F := H;
while member(0, F, 'p') do
   F := subsop(p = NULL, F);
   unassign('p')
   end do;
   L := LeadingMonomial(F, T);
   LL := LeadingMonomial(Basis(A, T), T);
   print(`minus`({op(LL)}, {op(L)}));
   if evalb(LeadingMonomial(`<,>`(op(L)), T) <> LeadingMonomial(`<,>`(op(LL)), T)) then
      RETURN(false):
   end if;
   for s to nops(A) do
      if evalb(Reduce(A[s], F, T) <> 0) then
         RETURN(false):
      end if:
   end do;
   RETURN(true):
end proc: