with(Groebner);
with(ArrayTools);

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; 
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])])); 
end do; 
RETURN(P); 
end proc;

BSH := proc(F, t) #This is the Berkesch-Schreyer algorithm to compute Grobner bases.
local G, P, p, s, h, g, i, LtG, firsttime, firstbytes, secondtime, secondbytes; 
firsttime, firstbytes := kernelopts(cputime, bytesused); 
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)); 
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); 
RETURN(convert(G, list)); 
end proc;

CommonFactor := proc(f) 
local A, gc, i; 
global lc; 
if f = 0 then 
RETURN(0); 
end if; 
A := [op(expand(f))]; 
if type(f, `+`) then 
gc := A[1]; 
for i from 2 to nops(A) do 
gc := gcd(gc, A[i]); 
end do; 
lc := lc*gc; 
RETURN(simplify(f/gc)); 
else lc := lc*f; 
RETURN(1); 
end if; 
end proc;

ImpBSH := proc(F, t)#This algorithm by using the Laurent setting computes a Grobner basis for an ideal between <F> and its saturation.
local G, P, p, s, h, i, f, LtG, g; 
global lc; 
lc := 1; 
G := Array([seq(CommonFactor(expand(f)), f in 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 := CommonFactor(SPolynomial(p[1], p[2], t)); 
h := CommonFactor(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)); 
end if; 
end do; 
RETURN(convert(G, list), [op(indets(lc))]); 
end proc;

Solving := proc(F, X, S, H, t) #This algorithm uses the Improved version of Berkesch-Schreyer algorithm to solve a polynomial system by computing all the nodes of the graph described in the paper.
local G, Y, s, i, T, g, v; 
global B, A, C; 
G, Y := ImpBSH(F, t); 
if nops(Y) = 0 then 
B := [op(B), [S, G]]; 
A := [op(A), S]; 
RETURN(); 
else 
s := nops(Y); 
for i to s do 
if member(sort([op(S), Y[i]]), A) = false then 
T := [seq(subs(Y[i] = 0, g), g in F)]; 
if member([sort(Y), "nonzero"], C) = false then 
C := [op(C), [sort(Y), "nonzero"]]; 
B := [op(B), [[Y, "nonzero"], G]]; 
end if; 
Solving(T, [op({op(X)} minus {Y[i]})], [op(S), Y[i]], [[[op(S), Y[i]], T]], t); 
end if; 
end do; 
end if; 
end proc;

Timing := proc(F, X, t) 
local firsttime, firstbytes, secondtime, secondbytes; 
global A, B, C; 
firsttime, firstbytes := kernelopts(cputime, bytesused); 
A := []; B := []; C := []; 
Solving(F, X, [], [], t); 
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); 
RETURN(B); 
end proc;

tord := plex(x, y);

F:=[x^12 - 66*x^10*y^2 + 495*x^8*y^4 - 924*x^6*y^6 + 495*x^4*y^8 - 66*x^2*y^10 + y^12 - 1, 12*x^11*y - 220*x^9*y^3 + 792*x^7*y^5 - 792*x^5*y^7 + 220*x^3*y^9 - 12*x*y^11];
Timing(F,[x,y],tord);
BSH(F,tord);

F:[x^13 - 78*x^11*y^2 + 715*x^9*y^4 - 1716*x^7*y^6 + 1287*x^5*y^8 - 286*x^3*y^10 + 13*x*y^12 - 1, 13*x^12*y - 286*x^10*y^3 + 1287*x^8*y^5 - 1716*x^6*y^7 + 715*x^4*y^9 - 78*x^2*y^11 + y^13];
Timing(F,[x,y],tord);
BSH(F,tord);

F:=[x^14 - 91*x^12*y^2 + 1001*x^10*y^4 - 3003*x^8*y^6 + 3003*x^6*y^8 - 1001*x^4*y^10 + 91*x^2*y^12 - y^14 - 1, 14*x^13*y - 364*x^11*y^3 + 2002*x^9*y^5 - 3432*x^7*y^7 + 2002*x^5*y^9 - 364*x^3*y^11 + 14*x*y^13];
Timing(F,[x,y],tord);
BSH(F,tord);


tord := tdeg(seq(x[i], i = 1 .. 8));

Katsura(7):=[2*x[7] + 2*x[6] + 2*x[5] + 2*x[4] + 2*x[3] + 2*x[2] + 2*x[1] + x[8] - 1, 2*x[1]*x[5] + 2*x[1]*x[7] + 2*x[2]*x[4] + x[2]*x[8] + x[3]^2 + 2*x[6]*x[8] - x[6], 2*x[1]*x[4] + 2*x[1]*x[6] + 2*x[2]*x[3] + 2*x[2]*x[7] + x[3]*x[8] + 2*x[5]*x[8] - x[5], 2*x[1]*x[3] + 2*x[1]*x[5] + x[2]^2 + 2*x[2]*x[6] + 2*x[3]*x[7] + 3*x[4]*x[8] - 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[8] + 2*x[4]*x[7] + x[5]*x[8] - x[3], 2*x[1]*x[2] + 2*x[1]*x[8] + 2*x[2]*x[3] + 2*x[3]*x[4] + 2*x[4]*x[5] + 2*x[5]*x[6] + 2*x[6]*x[7] + x[7]*x[8] - 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 + 2*x[7]^2 + x[8]^2 - x[8], x[1]^2 + 2*x[1]*x[3] + 2*x[2]*x[4] + 2*x[2]*x[8] + 2*x[3]*x[5] + 2*x[4]*x[6] + 2*x[5]*x[7] + x[6]*x[8] - x[2]];
Timing(F,[seq(x[i], i = 1 .. 8)],tord);
BSH(F,tord);