with(Groebner);
with(ArrayTools);
with(PolynomialIdeals);
with(ListTools);

Extendd := proc(A, B)
RETURN(Array([op(convert(A, list)), op(convert(B, list))])); 
end proc;

Appendd := proc(A, h) 
RETURN(Array([op(convert(A, list)), h])); 
end proc;

Sort := proc(L, t) 
RETURN(sort(L, (b, a) -> TestOrder(a[3], b[3], t))); 
end proc;

Pair := proc(G, LtG, g, t)
local P, A, A1, a, i, j, u, n, C; 
global pount; P := Array([]); 
u := LeadingMonomial(g, t); 
n := ArrayNumElems(LtG); 
A := Array([seq(lcm(LtG[j], u)/u, j = 1 .. n)]); 
A1 := Array([op({op(Basis(convert(A, list), t))} minus {seq(LtG[j], j = 1 .. n)})]); 
for a in A1 do member(a, A, 'q'); 
P := Appendd(P, Array([g, G[q], lcm(u, LtG[q]), m, u, LtG[q]])); 
pount := pount + 1; 
end do; 
RETURN(P); 
end proc;

BSH := proc(F, t)  # Berkesch and Schreyer algorithm.
local G, LtG, P, i, p, s, h, g, f, n, firsttime, firstbytes, secondtime, secondbytes; 
global count, pount, rount, pk, m; 
firsttime, firstbytes := kernelopts(cputime, bytesused); 
m := 0; 
pount := 0; 
rount := 0; 
G := Array(F); 
LtG := Array([seq(LeadingMonomial(g, t), g in G)]); 
P := Array([]); 
for i from 2 to ArrayNumElems(G) do 
P := Extendd(P, Pair(G[1 .. i - 1], LtG[1 .. i - 1], G[i], t)); 
end do; 
while ArrayNumElems(P) <> 0 do 
P := Sort(P, t); 
p := P[-1]; 
P := Array(1 .. ArrayNumElems(P) - 1, P); 
s := SPolynomial(p[1], p[2], t); 
h := NormalForm(s, convert(G, list), t); 
if h <> 0 then 
P := Extendd(P, Pair(G, LtG, h, t)); 
G := Appendd(G, h); 
LtG := Appendd(LtG, LeadingMonomial(h, t)); 
else 
rount := rount + 1; 
end if; end do; 
secondtime, secondbytes := kernelopts(cputime, bytesused); 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 %1s:         %3a\n", number, of, pairs, is, pount); 
printf("%-1s %1s %1s %1s:         %3a\n", number, of, zeros, is, rount); 
RETURN(G); 
end proc;

Criterion23 := proc(LtG, p, t) 
local n, r, u; 
n := p[4]; 
if ArrayNumElems(LtG) <= n + 1 then 
RETURN(true); 
end if; 
r := NormalForm(p[3], convert(LtG[n + 1 .. -1], list), t, 'k'); 
u := SearchArray(Array(k), location = last); 
if r <> 0 or lcm(LtG[u[1] + n], p[5]) = p[3] or lcm(LtG[u[1] + n], p[6]) = p[3] then 
RETURN(true); 
else 
RETURN(false); 
end if; 
end proc;

BSH23 := proc(F, t) #Improved version of Berkesch and Schreyer algorithm by using our new result to detect as much as possible Buchberger's second criterion.
local G, LtG, P, i, p, s, h, g, f, p23, n, firsttime, firstbytes, secondtime, secondbytes; 
global count, pount, rount, pk, m; 
firsttime, firstbytes := kernelopts(cputime, bytesused); 
p23 := 0; 
pount := 0; 
rount := 0; 
G := Array(F); 
m := ArrayNumElems(G); 
LtG := Array([seq(LeadingMonomial(g, t), g in G)]); 
P := Array([]); 
for i from 2 to ArrayNumElems(G) do 
P := Extendd(P, Pair(G[1 .. i - 1], LtG[1 .. i - 1], G[i], t)); 
end do; 
while ArrayNumElems(P) <> 0 do 
P := Sort(P, t); 
p := P[-1]; 
P := Array(1 .. ArrayNumElems(P) - 1, P); 
if Criterion23(LtG, p, t) = true then 
s := SPolynomial(p[1], p[2], t); 
h := NormalForm(s, convert(G, list), t); 
if h <> 0 then 
m := m + 1; 
P := Extendd(P, Pair(G, LtG, h, t)); 
G := Appendd(G, h); 
LtG := Appendd(LtG, LeadingMonomial(h, t)); 
else 
rount := rount + 1; 
end if; 
else 
p23 := p23 + 1; 
end if; 
end do; 
secondtime, secondbytes := kernelopts(cputime, bytesused); 
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 %1s:         %3a\n", number, of, pairs, is, pount); 
printf("%-1s %1s %1s %1s:         %3a\n", number, of, zeros, is, rount); 
printf("%-1s %1s %1s %1s %1s:   %3a\n", two, out, of, tree, criteria, p23); RETURN(G); end proc;

tord:=plex(seq(x[i], i = 1 .. 10));

F:=[-x[1]*x[2]^3*x[4]^3*x[5]^3*x[6]*x[7]^2*x[9]^2*x[10] + x[1]^2*x[3]^3*x[4]^3*x[5]*x[6]^2*x[7], -x[1]^3*x[2]*x[3]*x[4]^3*x[5]^3*x[7]^2*x[8]^2 - x[1]^3*x[2]*x[5]^2*x[6]*x[9]^2, -x[1]^2*x[4]^2*x[5]^3*x[6]^2*x[8]^3*x[9]*x[10] + x[1]^3*x[2]^2*x[4]*x[6]^3*x[10]];
BSH(F,tord);
BSH23(F,tord);

F:=[-x[1]*x[2]^3*x[4]^3*x[5]^3*x[6]*x[7]^2*x[9]^2*x[10] + x[1]^2*x[3]^3*x[4]^3*x[5]*x[6]^2*x[7], -x[1]^3*x[2]*x[3]*x[4]^3*x[5]^3*x[7]^2*x[8]^2 - x[1]^3*x[2]*x[5]^2*x[6]*x[9]^2, -x[1]^2*x[4]^2*x[5]^3*x[6]^2*x[8]^3*x[9]*x[10] + x[1]^3*x[2]^2*x[4]*x[6]^3*x[10]];
BSH(F,tord);
BSH23(F,tord);

F:=[-x[1]^2 x[2]^2 x[3]^7 x[4]^5 x[5]^2+x[1]^2 x[2] x[3]^4 x[4]^4 x[5]^4,x[1]^7 x[2]^6 x[3]^4 x[4]^4-x[1]^7 x[2]^4 x[3]^3 x[4]^2 x[5]^5,x[1]^4 x[2]^3 x[3] x[4]^7 x[5]^6-x[1]^5 x[2] x[3]^4 x[4]^3 x[5]^5,-x[1]^3 x[2]^2 x[3]^7 x[4]^3 x[5]^3+x[1] x[2]^5 x[3]^6];
BSH(F,tord);
BSH23(F,tord);