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

Update := proc(h, p, t) 
local F, F2, Lmp, r, i, Lcp, Tip, a, g; 
Lcp, Lmp := LeadingTerm(p, t); 
Tip := p - Lcp*Lmp; 
if type(h, monomial) then 
F := [[LeadingTerm(h, t)]]; 
F2 := [F[1][2]]; 
else 
F := [seq([LeadingTerm(a, t)], a in [op(h)])]; 
F2 := [seq(F[i][2], i = 1 .. nops(F))]; 
end if; 
if member(Lmp, F2, 'q') then 
F := subsop(q = [-F[q][1]/Lcp, Tip], F); 
r := expand(add(g[1]*g[2], g in F)); 
else 
r := h; 
end if; 
RETURN(r); 
end proc;

SorT := proc(L, t) 
RETURN(sort(L, (a, b) -> TestOrder(LeadingMonomial(expand(a[1]*a[2]*a[3]), t), LeadingMonomial(expand(b[1]*b[2]*b[3]), t), t))); 
end proc;

SATTraverso := proc(G, B, f, t, n) 
local i, U, p, P, Lmp, h, m, LtG, H, A, Syz, criterion, zero, pair, firsttime, firstbytes, secondtime, secondbytes; 
firsttime, firstbytes := kernelopts(cputime, bytesused); 
criterion := 0; 
zero := 0; 
pair := 0; 
A := []; 
Syz := []; 
U := B; 
H := []; 
LtG := LeadingMonomial(G, t); 
P := [[1, 1, f, f]]; 
pair := pair + 1; 
while P <> [] do 
P := SorT(P, t); 
m := P[1]; 
P := P[2 .. -1]; 
if member(m[1]*m[2], A) = false and NormalForm(m[1]*m[2], Syz, t) <> 0 then 
A := [op(A), m[1]*m[2]]; 
p := NormalForm(expand(m[1]*m[4]), G, t); 
for h in H do 
p := Update(p, h, t); 
end do; 
if p <> 0 then 
Lmp := LeadingMonomial(p, t); 
U := [op({op(U)} minus {Lmp})]; 
for i to n do 
if NormalForm(x[i]*m[1]*m[2], LtG, t) <> 0 then 
P := [op(P), [x[i], m[1]*m[2], f, p]]; 
pair := pair + 1; 
else 
pair := pair + 1; 
criterion := criterion + 1; 
end if; 
end do; 
if evalb(w in indets(Lmp)) = false then 
if NormalForm(w*m[1]*m[2], LtG, t) <> 0 then 
P := [op(P), [w, m[1]*m[2], f, p]]; 
pair := pair + 1; 
else 
pair := pair + 1; 
criterion := criterion + 1; 
end if; 
end if; 
H := [op(H), p]; 
else 
Syz := [op(Syz), m[1]*m[2]]; 
zero := zero + 1; 
end if; 
else
criterion := criterion + 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, pair); 
printf("%-1s %1s %1s %1s:         %3a\n", number, of, zeros, is, zero); 
printf("%-1s %1s %1s:             %3a\n", Two, criteria, is, criterion); 
RETURN([op(G), op(H)]); 
end proc;

lt := proc(f, T) 
RETURN(LeadingTerm(f, T)[1]*LeadingTerm(f, T)[2]); 
end proc;

FL := proc(f, T) 
local L, p; 
L := []; 
p := f; 
while p <> 0 do 
L := [op(L), lt(p, T)]; 
p := expand(p - lt(p, T)); 
end do; 
RETURN(L); 
end proc;

Gauss := proc(FF, t) 
local A, B, C, G, F, u, i; 
F := FF; 
G := [seq(FL(u, t), u = F)]; 
A := [seq(seq(LeadingMonomial(u, t), u = G[i]), i = 1 .. nops(G))]; 
A := [op({op(A)})]; 
A := sort(A, (a, b) -> TestOrder(a, b, t)); 
if evalb(A = [1]) = true then 
RETURN(A); 
elif F = [] then 
RETURN(F); 
end if; 
if A[1] = 1 then 
A := A[2 .. -1]; 
end if; 
if A <> [] then 
for i from nops(A) by -1 to 1 do 
F := subs(A[i] = Y[i], F); 
end do; 
F := Basis(F, plex(seq(Y[i], i = nops(A) .. 1, -1))); 
for i from nops(A) by -1 to 1 do 
F := subs(Y[i] = A[i], F); 
end do; 
end if; 
end proc;

SORT := proc(a, b, t) 
if a[1] < b[1] then 
RETURN(true); 
elif a[1] = b[1] and TestOrder(a[3], b[3], t) = true then RETURN(true); 
else 
RETURN(false); 
end if; 
end proc;

sorT := proc(L, t) 
RETURN(sort(L, (a, b) -> SORT(a, b, t))); 
end proc;

EXTraverso := proc(G, B, F, t, n) 
local i, U, p, h, P, Lmp, k, m, f, g, LtG, H, A, Syz, s, vp; 
global zero, pair, criterion; 
A := [seq([], i = 1 .. nops(F))]; 
Syz := [seq([], i = 1 .. nops(F))]; 
U := B; 
H := []; 
LtG := LeadingMonomial(G, t); 
P := [seq([i, 1, 1, F[i]], i = 1 .. nops(F))]; 
pair := pair + nops(F); 
while P <> [] do 
P := sort(P); 
P := sorT(P, t); 
m := P[1]; 
P := P[2 .. -1]; 
if member(m[2]*m[3], A[m[1]]) = false and NormalForm(m[2]*m[3], Syz[m[1]], t) <> 0 then 
A := subsop(m[1] = [op(A[m[1]]), m[2]*m[3]], A); 
p := NormalForm(expand(m[2]*m[4]), G, t); 
for h in H do 
p := Update(p, h, t); 
end do; 
if p <> 0 then 
Lmp := LeadingMonomial(p, t); 
U := [op({op(U)} minus {Lmp})]; 
for i to n do 
if NormalForm(x[i]*m[2]*m[3], LtG, t) <> 0 then 
P := [op(P), [m[1], x[i], m[2]*m[3], p]]; 
pair := pair + 1; 
else 
criterion := criterion + 1; 
pair := pair + 1; 
end if; 
end do; 
H := [op(H), p]; 
else 
Syz := subsop(m[1] = [op(Syz[m[1]]), m[2]*m[3]], Syz); 
zero := zero + 1; 
end if; 
else 
criterion := criterion + 1; 
end if; 
end do; 
RETURN([op(G), op(H)], U); 
end proc;

Saturation := proc(G, B, t, n) 
local f, i, F, H, R, W, U, TF, satnum, h, CH, flag, firsttime, firstbytes, secondtime, secondbytes, b, g; 
global zero, pair, criterion; 
firsttime, firstbytes := kernelopts(cputime, bytesused); 
zero := 0; 
criterion := 0; 
pair := 0; 
satnum := 0; 
f := w*mul(x[i], i = 1 .. n) - 1; 
R := G; 
U := B; 
F := NormalForm([seq(expand(b*f), b in U)], R, t); 
pair := pair + nops(F); 
flag := false; 
while flag = false do 
H := Gauss(F, t); 
W, CH := selectremove(has, H, w); 
if CH <> [] then 
TF := EXTraverso(R, U, CH, t, n); 
R := TF[1]; 
U := TF[2]; 
satnum := satnum + 1; 
F := [seq(NormalForm(g, R, t), g in W)]; 
pair := pair + nops(F); 
else 
flag := true; 
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, pair); 
printf("%-1s %1s %1s %1s:         %3a\n", number, of, zeros, is, zero); 
printf("%-1s %1s %1s:            %3a\n", Two, criteria, is, criterion); 
RETURN(R); 
end proc;

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 pair, criterion; 
P := Array([]); 
u := LeadingMonomial(g, t); 
n := ArrayNumElems(LtG); 
pair := pair + n; 
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)})]); 
criterion := criterion + n - ArrayNumElems(A1); 
for a in A1 do 
member(a, A, 'q'); 
P := Appendd(P, Array([g, G[q], lcm(u, LtG[q])])); 
end do; 
RETURN(P); 
end proc;

BSH := proc(F, f, t) 
local G, LtG, P, zero, i, p, s, h, g, n, firsttime, firstbytes, secondtime, secondbytes; 
global pair, criterion; 
firsttime, firstbytes := kernelopts(cputime, bytesused); 
criterion := 0; 
pair := 0; 
zero := 0; 
G := Array(F); 
LtG := Array(LeadingMonomial(convert(G, list), t)); 
P := Pair(G, LtG, f, t); 
G := Appendd(G, f); 
LtG := Appendd(LtG, LeadingMonomial(f, t)); 
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 
zero := zero + 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, pair); 
printf("%-1s %1s %1s %1s:         %3a\n", number, of, zeros, is, zero); 
printf("%-1s %1s %1s:           %3a\n", Buch, criteria, is, criterion); 
RETURN(convert(G, list)); 
end proc;


Eco7 := [x[1]*x[2]*x[7] + x[2]*x[3]*x[7] + x[3]*x[4]*x[7] + x[4]*x[5]*x[7] + x[5]*x[6]*x[7] + x[1]*x[7] - 1, x[1]*x[3]*x[7] + x[2]*x[4]*x[7] + x[3]*x[5]*x[7] + x[4]*x[6]*x[7] + x[2]*x[7] - 2, x[1]*x[4]*x[7] + x[2]*x[5]*x[7] + x[3]*x[6]*x[7] + x[3]*x[7] - 3, x[1]*x[5]*x[7] + x[2]*x[6]*x[7] + x[4]*x[7] - 4, x[1]*x[6]*x[7] + x[5]*x[7] - 5, x[6]*x[7] - 6, x[1] + x[2] + x[3] + x[4] + x[5] + x[6] + 1];
t := tdeg(seq(x[i], i = 1 .. 7));
G := Basis(Eco7, t);
B := NormalSet(G, t)[1];
h := NormalForm(w*mul(x[i], i = 1 .. 7) - 1, G, prod(plex(w), t));
SATTraverso(G, B, h, prod(plex(w), t), 7);
Saturation(G, B, prod(plex(w), t), 7);
BSH(G, w*mul(x[i], i = 1 .. 7) - 1, prod(plex(w), t));



Kat5 := [2*x[5] + 2*x[4] + 2*x[3] + 2*x[2] + 2*x[1] + x[6] - 1, 2*x[1]*x[2] + 2*x[1]*x[4] + 2*x[2]*x[5] + 3*x[3]*x[6] - x[3], 2*x[1]*x[3] + 2*x[1]*x[5] + x[2]^2 + x[2]*x[6] + 2*x[4]*x[6] - x[4], 2*x[1]*x[2] + 2*x[1]*x[6] + 2*x[2]*x[3] + 2*x[3]*x[4] + 2*x[4]*x[5] + x[5]*x[6] - x[1], 2*x[1]^2 + 2*x[2]^2 + 2*x[3]^2 + 2*x[4]^2 + 2*x[5]^2 + x[6]^2 - x[6], x[1]^2 + 2*x[1]*x[3] + 2*x[2]*x[4] + 2*x[2]*x[6] + 2*x[3]*x[5] + x[4]*x[6] - x[2]];
t := tdeg(seq(x[i], i = 1 .. 6));
G := Basis(Kat5, t);
B := NormalSet(G, t)[1];
h := NormalForm(w*mul(x[i], i = 1 .. 6) - 1, G, prod(plex(w), t));
SATTraverso(G, B, h, prod(plex(w), t), 6);
Saturation(G, B, prod(plex(w), t), 6);
BSH(G, w*mul(x[i], i = 1 .. 6) - 1, prod(plex(w), t));


Bellido := [x[1]^2 + x[2]^2 + x[3]^2 - 12*x[1] - 68, x[4]^2 + x[5]^2 + x[6]^2 - 12*x[5] - 68, x[7]^2 + x[8]^2 + x[9]^2 - 24*x[8] - 12*x[9] + 100, x[1]*x[4] + x[2]*x[5] + x[3]*x[6] - 6*x[1] - 6*x[5] - 52, x[1]*x[7] + x[2]*x[8] + x[3]*x[9] - 6*x[1] - 12*x[8] - 6*x[9] + 64, x[4]*x[7] + x[5]*x[8] + x[6]*x[9] - 6*x[5] - 12*x[8] - 6*x[9] + 32, 2*x[2] + 2*x[3] - 2*x[6] - x[4] - x[5] - x[7] - x[9] + 18, x[1] + x[2] + 2*x[3] + 2*x[4] + 2*x[6] - 2*x[7] + x[8] - x[9] - 38, x[1] + x[3] + x[5] - x[6] + 2*x[7] - 2*x[8] - 2*x[4] + 8];
t := tdeg(seq(x[i], i = 1 .. 9));
G := Basis(Bellido, t);
B := NormalSet(G, t)[1];
h := NormalForm(w*mul(x[i], i = 1 .. 9) - 1, G, prod(plex(w), t));
SATTraverso(G, B, h, prod(plex(w), t), 9);
Saturation(G, B, prod(plex(w), t), 9);
BSH(G, w*mul(x[i], i = 1 .. 9) - 1, prod(plex(w), t));


CR11 := [-8*x[1]*x[2] - 2*x[1] - x[3] + 1, -8*x[1]*x[2] - 2*x[2] - 3*x[4], -8*x[3]*x[4] + x[1] - 2*x[3] - x[5], -8*x[3]*x[4] + x[2] - 2*x[4] - 3*x[6], -8*x[5]*x[6] + x[3] - 2*x[5] - x[7], -8*x[5]*x[6] + x[4] - 2*x[6] - 3*x[8], -8*x[7]*x[8] + x[5] - 2*x[7] - x[9], -8*x[7]*x[8] + x[6] - 2*x[8] - 3*x[10], -8*x[9]*x[10] + x[7] - 2*x[9] - x[11], -8*x[10]*x[11] + x[8] - 2*x[10], -8*x[10]*x[11] + x[9] - 2*x[11] - 3];
t := tdeg(seq(x[i], i = 1 .. 11));
G := Basis(CR11, t);
B := NormalSet(G, t)[1];
h := NormalForm(w*mul(x[i], i = 1 .. 11) - 1, G, prod(plex(w), t));
SATTraverso(G, B, h, prod(plex(w), t), 11);
Saturation(G, B, prod(plex(w), t), 11);
BSH(G, w*mul(x[i], i = 1 .. 11) - 1, prod(plex(w), t));


t := tdeg(seq(x[i], i = 1 .. 2));
V := [[-10, -7], [-10, -2], [-10, 2], [-10, 5], [-9, -6], [-9, 0], [-9, 1], [-9, 2], [-8, -9], [-8, -8], [-8, 0], [-8, 7], [-7, -8], [-7, -7], [-7, -4], [-7, -3], [-7, 7], [-7, 10], [-6, 10], [-5, -10], [-5, -9], [-5, -4], [-5, 0], [-5, 5], [-5, 9], [-4, -8], [-4, 6], [-3, -8], [-3, -6], [-3, -2], [-2, -10], [-2, -5], [-2, 10], [-1, -10], [-1, -8], [-1, 3], [-1, 7], [-1, 10], [0, -9], [0, -3], [0, -2], [0, 3], [0, 8], [1, 1], [1, 3], [1, 5], [1, 7], [2, 2], [2, 4], [2, 5], [2, 8], [2, 9], [3, -5], [3, -3], [3, -2], [3, 2], [3, 5], [4, -8], [4, -4], [4, -1], [4, 0], [4, 1], [4, 3], [4, 6], [4, 10], [5, -10], [5, -5], [6, -9], [6, -7], [6, -6], [6, -3], [6, 4], [6, 8], [7, -5], [7, 2], [7, 4], [7, 6], [8, -8], [8, -7], [8, -5], [8, -4], [8, -3], [8, -2], [8, -1], [8, 7], [10, -10], [10, -6], [10, 8], [10, 10]];
J := VanishingIdeal(V, [x[1], x[2]]);
G := Basis(J, t);
B := NormalSet(G, t)[1];
h := NormalForm(w*mul(x[i], i = 1 .. 2) - 1, G, prod(plex(w), t));
SATTraverso(G, B, h, prod(plex(w), t), 2);
Saturation(G, B, prod(plex(w), t), 2);
BSH(G, w*mul(x[i], i = 1 .. 2) - 1, prod(plex(w), t));


Kat6 := [2*x[6] + 2*x[5] + 2*x[4] + 2*x[3] + 2*x[2] + 2*x[1] + x[7] - 1, 2*x[1]*x[4] + 2*x[1]*x[6] + 2*x[2]*x[3] + x[2]*x[7] + 2*x[5]*x[7] - x[5], 2*x[1]*x[3] + 2*x[1]*x[5] + x[2]^2 + 2*x[2]*x[6] + x[3]*x[7] + 2*x[4]*x[7] - x[4], 2*x[1]*x[2] + 2*x[1]*x[4] + 2*x[2]*x[5] + 2*x[3]*x[6] + 2*x[3]*x[7] + x[4]*x[7] - x[3], 2*x[1]*x[2] + 2*x[1]*x[7] + 2*x[2]*x[3] + 2*x[3]*x[4] + 2*x[4]*x[5] + 2*x[5]*x[6] + x[6]*x[7] - x[1], 2*x[1]^2 + 2*x[2]^2 + 2*x[3]^2 + 2*x[4]^2 + 2*x[5]^2 + 2*x[6]^2 + x[7]^2 - x[7], x[1]^2 + 2*x[1]*x[3] + 2*x[2]*x[4] + 2*x[2]*x[7] + 2*x[3]*x[5] + 2*x[4]*x[6] + x[5]*x[7] - x[2]];
t := tdeg(seq(x[i], i = 1 .. 7));
G := Basis(Kat6, t);
B := NormalSet(G, t)[1];
h := NormalForm(w*mul(x[i], i = 1 .. 7) - 1, G, prod(plex(w), t));
SATTraverso(G, B, h, prod(plex(w), t), 7);
Saturation(G, B, prod(plex(w), t), 7);
BSH(G, w*mul(x[i], i = 1 .. 7) - 1, prod(plex(w), t));