with(Groebner); with(LinearAlgebra); with(ListTools); with(ArrayTools);

ord := proc (i, Q, P, E)
local A, M, N, e, j, a, l, L, f, lt, u, leadt, C, V, G, s, v1, o, p, S, c, K, b;
K := Array(Q);
C := K[1];
V := K[2];
G := K[3];
p := nops(P);
c := ArrayNumElems(C);
S := P;
A := [];
M := [];
N := [];
e := [];
L := [];
b := i;
if p < b or b <= 0 then print(" the order is not defined") else
s := `union`(seq({op(S[j])}, j = 1 .. p));
f := AddAlongDimension(ElementMultiply(ElementMultiply(C, V), G), 1);
o := LeadingTerm(f, plex(op(s)));
if o[1]*o[2] = f then leadt := (Q[2][1], Q[3][1], 1) else
for j to c do
a[1] := degree(V(j), S[b]);
a[2] := degree(V(j)); a[3] := seq(degree(V(j), S[t]), t = 1 .. b-1);
a[4] := seq(degree(V(j), S[t]), t = b+1 .. p);
a[5] := seq(degree(V(j), S[b][t]), t = 1 .. nops(S[b]));
M := [seq(op(S[t]), t = 1 .. b-1)]; N := [seq(op(S[t]), t = b+1 .. p)];
e := [op(M), op(N)];
a[6] := seq(degree(V(j), e[t]), t = 1 .. nops(e));
l := [a[1], a[2], a[3], a[4], a[5], a[6]];
A := [op(A), l]
end do;
L := Array([seq(mul(v[t]^A[j][t], t = 1 .. nops(A[j])), j = 1 .. nops(A))]);
f := AddAlongDimension(ElementMultiply(L, G), 1);
lt := LeadingMonomial(f, plex(seq(v[t], t = 1 .. nops(A[1])), op(E)));
u := [seq(degree(lt, v[t]), t = 1 .. nops(A[1]))];
for j to nops(A) do if u = A[j] then leadt := (V(j), G(j), j)
end if;
end do;
end if;
end if;
return [leadt];
end proc;
#//////////////////////////////////////////////////////////////////////////////////////////////////////
dlcm := proc (L, K)
local A, B, t, u, LCM;
A := L;
B := K;
t := A[2];
u := B[2];
if t <> u then LCM := 0 else
LCM := lcm(A[1], B[1])*A[2]
end if;
return LCM
end proc;
#//////////////////////////////////////////////////////////////////////////////////////////////////////////
S := proc (n, L, K, H, E, P)
local A, B, c, k, l, j, z, M, a, t, q, f, N, k1, l1, z1, a1, q1, g, d1, d2, q2, q3, q4, sq1, sq2, sq3, sq4, sq5, sq6, sq7, sq, i, nl, dl, nl1, dl1, p1, p2, F, u, u1, n1, n2,n3,m1,m2,m3;
F := [seq(op(P[i]), i = 1 .. nops(P))];
M := Array(L);
N := Array(K);
A := ord(n, M, P, E);
B := ord(n, N, P, E);
c := dlcm(A, B);
p1 := ArrayNumElems(M[1]);
p2 := ArrayNumElems(N[1]);
if c = 0 then return 0 end if;
k := c/(A[1]*A[2]);
l := M[1][A[3]];
if k <> 1 then
for j to nops(F) do
if divide(k, F[j]) and member(H[j], indets(l)) = true then
z := subs(H[j] = H[j]+degree(k, F[j]), l);
l := z
end if;
end do;
a := [];
u := [];
u := Vector[row](u);
for t to p1 do
u(1, t) := M[1][t];
for j to nops(F) do
if divide(k, F[j]) and member(H[j], indets(u(1, t))) = true then
a := subs(H[j] = H[j]+degree(k, F[j]), u(1, t));
u(1, t) := a
end if;
end do;
end do;
m1 := u;
m2 := ElementMultiply(k, M[2]);
m3 := M[3];
f := AddAlongDimension(ElementMultiply(ElementMultiply(m1, m2), m3), 1) else
f := AddAlongDimension(ElementMultiply(ElementMultiply(M[1], M[2]), M[3]), 1)
end if;
k1 := c/(B[1]*B[2]);
l1 := N[1][B[3]];
if k1 <> 1 then
for j to nops(F) do
if divide(k1, F[j]) and member(H[j], indets(l1)) = true then
z1 := subs(H[j] = H[j]+degree(k1, F[j]), l1);
l1 := z1
end if;
end do;
a1 := [];
u1 := [];
u1 := Vector[row](u1);
for t to p2 do
u1(1, t) := N[1][t];
for j to nops(F) do
if divide(k1, F[j]) and member(H[j], indets(u1(1, t))) = true then
a1 := subs(H[j] = H[j]+degree(k1, F[j]), u1(1, t));
u1(1, t) := a1
end if;
end do;
end do;
n1 := u1;
n2 := ElementMultiply(k1, N[2]);
n3 := N[3];
g := AddAlongDimension(ElementMultiply(ElementMultiply(n1, n2), n3), 1) else
g := AddAlongDimension(ElementMultiply(ElementMultiply(N[1], N[2]), N[3]), 1)
end if;
d1 := [coeffs(f, [op(F), op(E)])];
d2 := [coeffs(g, [op(F), op(E)])];
sq1 := [seq(numer(d1[j]), j = 1 .. nops(d1))];
sq2 := [seq(numer(d2[j]), j = 1 .. nops(d2))];
sq3 := [seq(denom(d1[j]), j = 1 .. nops(d1))];
sq4 := [seq(denom(d2[j]), j = 1 .. nops(d2))];
for i to nops(sq1) do
if irreduc(sq1[i]) = false then
subs(sq1[i] = op(factor(sq1[i])), sq1)
end if;
end do;
for i to nops(sq3) do
if irreduc(sq3[i]) = false then
subs(sq3[i] = op(factor(sq3[i])), sq3)
end if;
end do;
for i to nops(sq2) do
if irreduc(sq2[i]) = false then
subs(sq2[i] = op(factor(sq2[i])), sq2)
end if;
end do;
for i to nops(sq4) do
if irreduc(sq4[i]) = false then
subs(sq4[i] = op(factor(sq4[i])), sq4)
end if;
end do;
nl := numer(l);
dl := denom(l);
nl1 := numer(l1);
dl1 := denom(l1);
if irreduc(nl) = false then
nl := [op(factor(nl))]
end if;
if irreduc(nl1) = false then
nl1 := [op(factor(nl1))]
end if;
if irreduc(dl) = false then
dl := [op(factor(dl))]
end if;
if irreduc(dl1) = false then
dl1 := [op(factor(dl1))]
end if;
sq7 := {op(sq1)} union {op(sq2)} union {op(sq3)} union {op(sq4)} union {op(nl)} union {op(nl1)} union {op(dl)} union {op(dl1)};
sq := [op(sq7)];
for i to nops(sq) do
if type(sq[i], monomial) = true then
sq7 := sq7 minus {sq[i]}
end if;
end do;
q2 := expand(f/l-g/l1, op(sq7));
q3 := [op(collect(q2, [op(F), op(E)], distributed))];
q4 := add(simplify(coeffs(q3[i], [op(F), op(E)]))*q3[i]/coeffs(q3[i], [op(F), op(E)]), i = 1 .. nops(q3));
return q4;
end proc;
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
com := proc (L, K)
local j, M;
j := 1;
M := [];
while j <= nops(L) and L[j] <= K[j] do
M := [op(M), L[j]]; j := j+1
end do;
return M;
end proc;
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
red := proc (f, G, H, E, T, P1)
local O, r, g, j, t, k, l, i, z1, a, w, c, q, f1, f2, f3, f4, f5, H2, H1, p, f6, f7, d, P, H3, H4, r1, r2, F1, Ha, G1, q2, q1, g1, f31, co1, co2, g2, d1, d2, sq1, sq2, sq3, sq4, sq11, sq21, sq31, sq41, sq7, sq, r3, U, sw, z3, o, v1, t3, P2, t4, U1, u, m1, m2, m3, F, q3, f_4;
if f = 0 then return 0 end if;
r2 := nops(P1);
r3 := 0;
d := f;
F := [seq(op(P1[i]), i = 1 .. nops(P1))];
F1 := F;
Ha := H;
U := G;
if r2 < nops(T) then print("This reduction is not defined") else
v1 := [seq(d[1][j]*d[2][j]*d[3][j], j = 1 .. nops(d[1]))];
f4 := add(v1[j], j = 1 .. nops(v1));
p := [seq(ord(T[1], U[j], P1, E), j = 1 .. Size(G)[2])];
r1 := nops(p);
O := [seq(p[j][1]*p[j][2], j = 1 .. r1)];
while f4 <> 0 do
i := 1;
sw := false;
while i <= r1 and sw = false do
U := G;
Ha := H;
P2 := ord(T[1], d, P1, E);
P := P2[1]*P2[2];
if divide(P, O[i]) = true then
z3 := P2[3];
if nops(T) <> 1 then
co1 := [seq(degree(P*ord(T[m], U[i], P1, E)[1]/O[i], {op(P1[T[m]])}), m = 2 .. nops(T))];
co2 := [seq(degree(ord(T[n], d, P1, E)[1], {op(P1[T[n]])}), n = 2 .. nops(T))] else
co1 := [0];
co2 := [0]
end if;
if com(co1, co2) = co1 then k := P/O[i]; l := U[i][1, p[i][3]];
if k <> 1 then for t to nops(F1) do
if divide(k, F1[t]) and member(Ha[t], indets(l)) = true then
z1 := subs(Ha[t] = Ha[t]+degree(k, F1[t]), l); l := z1
end if;
end do;
a := [];
u := [];
u := Vector[row](u);
for w to ArrayNumElems(Array(U[i][1])) do
u(1, w) := U[i][1][w];
for t to nops(F1) do
if divide(k, F1[t]) and member(Ha[t], indets(u(1, w))) = true
then a := subs(Ha[t] = Ha[t]+degree(k, F1[t]), u(1, w));
u(1, w) := a
end if;
end do;
end do;
m1 := u;
m2 := ElementMultiply(k, U[i][2]);
m3 := U[i][3];
q := AddAlongDimension(ElementMultiply(ElementMultiply(d[1][z3]*m1/l, m2), m3), 1) else
q := AddAlongDimension(ElementMultiply(ElementMultiply(d[1][z3]*U[i][1]/l, U[i][2]), U[i][3]), 1)
end if;
g1 := [seq(d[1][j]*d[2][j], j = 1 .. nops(d[2]))];
g := Multiply(convert(g1, Matrix), LinearAlgebra[Transpose](convert(d[3], Matrix)));
g2 := op(convert(g, list));
d1 := [coeffs(q, [op(F), op(E)])];
d2 := [coeffs(g2, [op(F), op(E)])];
sq1 := [op({seq(numer(d1[j]), j = 1 .. nops(d1))})];
sq2 := [op({seq(numer(d2[j]), j = 1 .. nops(d2))})];
sq3 := [op({seq(denom(d1[j]), j = 1 .. nops(d1))})];
sq4 := [op({seq(denom(d2[j]), j = 1 .. nops(d2))})];
sq11 := {};
sq21 := {};
sq31 := {};
sq41 := {};
for i to nops(sq1) do
if irreduc(sq1[i]) = false then
sq11 := {op(factor(sq1[i])), op(sq11)} else
sq11 := {sq1[i], op(sq11)}
end if;
end do;
for i to nops(sq3) do
if irreduc(sq3[i]) = false then sq31 := {op(factor(sq3[i])), op(sq31)} else
sq31 := {sq3[i], op(sq31)}
end if;
end do;
for i to nops(sq2) do
if irreduc(sq2[i]) = false
then sq21 := {op(factor(sq2[i])), op(sq21)} else
sq21 := {sq2[i], op(sq21)}
end if;
end do;
for i to nops(sq4) do
if irreduc(sq4[i]) = false then
sq41 := {op(factor(sq4[i])), op(sq41)} else
sq41 := {sq4[i], op(sq41)}
end if;
end do;
sq7 := sq11 union sq21 union sq31 union sq41;
sq := [op(sq7)];
for i to nops(sq) do
if type(sq[i], monomial) = true then
sq7 := sq7 minus {sq[i]}
end if;
end do;
f4 := expand(g2-q, op(sq7));
if simplify(f4) = 0 then
return collect(r3, [op(F), op(E)], distributed); else
f_4 := simplify(f4);
o := [LeadingTerm(f_4, plex(op(F)))];
if o[1]*o[2] = f_4 then
f5 := [f_4] else
q3 := [op(collect(f4, [op(F), op(E)], distributed))];
f4 := add(simplify(coeffs(q3[i], [op(F), op(E)]))*q3[i]/coeffs(q3[i], [op(F), op(E)]), i = 1 .. nops(q3));
f5 := [op(f4)]
end if;
end if;
H1 := subs(seq(F[i] = 1, i = 1 .. nops(F)), seq(E[i] = 1, i = 1 .. nops(E)), f5);
H2 := [seq(f5[c]/H1[c], c = 1 .. nops(f5))];
H3 := subs(seq(E[i] = 1, i = 1 .. nops(E)), H2);
H4 := subs(seq(F[i] = 1, i = 1 .. nops(F)), H2);
d := [H1, H3, H4]; sw := true else
i := i+1
end if;
else i := i+1
end if;
end do;
if sw = false then
t4 := ord(T[1], d, P1, E);
t3 := t4[3];
t := t4[1]*t4[2]*d[1][t3];
r3 := r3+t;
f4 := f4-t;
if simplify(f4) = 0 then return collect(r3, [op(F), op(E)], distributed) else
f_4 := simplify(f4);
o := [LeadingTerm(f_4, plex(op(F)))]; if o[1]*o[2] = f_4 then
f5 := [f_4] else
q3 := [op(collect(f4, [op(F), op(E)], distributed))];
f4 := add(simplify(coeffs(q3[i], [op(F), op(E)]))*q3[i]/coeffs(q3[i], [op(F), op(E)]), i = 1 .. nops(q3)); f5 := [op(f4)]
end if;
H1 := subs(seq(F[i] = 1, i = 1 .. nops(F)), seq(E[i] = 1, i = 1 .. nops(E)), f5);
H2 := [seq(f5[c]/H1[c], c = 1 .. nops(f5))];
H3 := subs(seq(E[i] = 1, i = 1 .. nops(E)), H2);
H4 := subs(seq(F[i] = 1, i = 1 .. nops(F)), H2);
d := [H1, H3, H4]
end if;
end if;
end do;
return collect(r3, [op(F), op(E)], distributed);
end if;
end proc;
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ 
zarb := proc (G) options operator, arrow; [seq(seq(Array([G[i], G[j]]), j = i+1 .. Size(G)[1]), i = 1 .. Size(G)[1])] end proc;
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
zarb2 := proc (A, B) 
local i, j, c, d; 
description "This algorithm computes product of two sets A and B"; 
c := []; 
d := []; 
for i to nops(A) do 
for j to nops(B) do 
if nops(B) = 1 then c := [op(c), [A[i], B[j]]] 
end if; 
if nops(A) = 1 then c := [op(c), [A[i], B[j]]] 
end if; 
if i <> j then c := [op(c), [A[i], B[j]]] 
end if; 
end do; 
end do; 
c := convert(c, set); 
return c; 
end proc;
  #\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
decom := proc (f, F, E) 
local o, f5, H1, H2, H3, H4, g, q3, d; 
d := []; 
g := f; 
if g <> 0 then o := [LeadingTerm(g, plex(op(F)))]; 
if o[1]*o[2] = g then f5 := [g] else q3 := [op(collect(g, [op(F), op(E)], distributed))]; 
g := add(simplify(coeffs(q3[i], [op(F), op(E)]))*q3[i]/coeffs(q3[i], [op(F), op(E)]), i = 1 .. nops(q3)); 
f5 := [op(g)]; 
end if; 
H1 := subs(seq(F[i] = 1, i = 1 .. nops(F)), seq(E[i] = 1, i = 1 .. nops(E)), f5); 
H2 := [seq(simplify(f5[c]/H1[c]), c = 1 .. nops(f5))]; 
H3 := subs(seq(E[i] = 1, i = 1 .. nops(E)), H2); 
H4 := subs(seq(F[i] = 1, i = 1 .. nops(F)), H2); 
d := [H1, H3, H4] 
else 
d := 0 
end if; 
return d; 
end proc;
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
buchalg_1 := proc (M, H, E, T, P, p) 
local t, G, r, h, i, K, l, q, F, G1; 
F := [seq(op(P[i]), i = 1 .. nops(P))]; 
G1 := Array(M); 
G := Array([seq(G1[j], j = 1 .. nops(M))]); 
K := Array(zarb(G1)); i := 1; 
while IsZero(K) = false do t := K[i]; 
if t <> 0 then K := subs(t = 0, K); 
h := S(p, t[1], t[2], H, E, P); 
h := decom(h, F, E); 
r := red(h, G, H, E, T, P); 
if r <> 0 then r := Array(decom(r, F, E)); 
l := Array([seq(Array([G[j], r]), j = 1 .. Size(G)[2])]); 
K := Array([seq(K[m], m = 1 .. Size(K)[2]), seq(l[j], j = 1 .. Size(l)[2])]); 
G := Array([seq(G[j], j = 1 .. Size(G)[2]), r]) 
end if; 
else i := i+1 
end if; 
end do; 
for i to Size(G)[2] do G[i] := [convert(G[i][1], list), convert(G[i][2], list), convert(G[i][3], list)] 
end do; 
return convert(G, list); 
end proc; 
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
buchalg := proc (M, H, E, T, P) 
local i, G, t1, t2, b1, b2, p; 
option trace; 
t1, b1 := kernelopts(cputime, bytesused); 
p := op(T[1]); 
G := buchalg_1(M, H, E, T[1], P, p); 
i := 1; 
while i <= nops(T)-1 do G := buchalg_1(G, H, E, T[i+1], P, p-i); 
i := i+1 
end do; 
t2, b2 := kernelopts(cputime, bytesused); 
print("CPU time:", t2-t1); 
print(); 
print("Used Memory:", b2-b1); 
return G; 
end proc; 
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\