with(LinearAlgebra):
with(Groebner):


####################################### 
###This function, returns exponent(f) with respect to n where n is all vairables

power := proc (f, n) 
p := f; 
v := NULL; 
while p <> 0 do 
lm := LeadingMonomial(p, plex(op(n))); 
lt := LeadingTerm(p, plex(op(n)))[1]*LeadingTerm(p, plex(op(n)))[2]; 
u := NULL; 
for i to nops(n) do 
u := u, degree(lm, n[i]) 
end do; 
v := v, Vector([u]); 
p := p-lt; 
end do;
v:=[v];
end proc;

#######################################
###This function, returns monomial defined by a vector

reversePowr := proc (v, n)
local s,i; 
s := 1; 
for i to nops(n) do 
s := s*n[i]^v[i];
end do; 
RETURN(s);
end proc;

#######################################
###These two functions, do updating set of pairs and also suger strategy for pair selection.

degSuger := proc(F,H)
f := F[1];
ltf := F[2];
suf := F[3];
h := H[1];
lth := H[2];
suh := H[3];

mf := lcm(ltf,lth)/ltf;
mh := lcm(ltf,lth)/lth;

degsu := max((degree(mf)+suf),(degree(mh)+suh));
RETURN(degsu);
end proc;

update1 := proc (G, B, h) 
local C, i, D, E, Gnew, Bnew, lth, temp, ltg, f1, hth, htg2, f2, t, Bh, Gh, g; 

C := NULL;
 
for i to nops(G) do 
ds := degSuger(h,G[i]);
C := C, [h, G[i],ds]; 
end do; 

C := [C]; 
D := NULL; 
while nops(C) <> 0 do 
temp := C[1]; 
lth := temp[1][2]; 
ltg := temp[2][2]; 
C := C[2 .. -1]; 
f1 := true; 
f2 := true; 
for i to nops(C) do 
t := C[i]; 
hth := t[1][2]; 
htg2 := t[2][2]; 
if divide(lcm(hth, ltg), lcm(lth, htg2)) = true then 
f1 := false; 
end if; 
end do; 
for i to nops([D]) do 
t := D[i]; 
hth := t[1][2]; 
htg2 := t[2][2]; 
if divide(lcm(hth, ltg), lcm(lth, htg2)) = true then 
f2 := false; 
end if; 
end do; 
if gcd(lth, ltg) = 1 or f1 = true and f2 = true then 
D := D, temp; 
end if; 
end do; 
D := [D]; 
E := NULL; 
while nops(D) <> 0 do 
temp := D[1]; 
D := D[2 .. -1]; 
if gcd(temp[1][2], temp[2][2]) <> 1 then 
E := E, temp; 
end if; 
end do; 
E := [E]; 
Bnew := NULL; 
Bh := B; 
while nops(Bh) <> 0 do 
temp := Bh[1]; 
Bh := Bh[2 .. -1]; 
if divide(lcm(temp[1][2], temp[2][2]), h[2]) <> true or lcm(temp[1][2], h[2]) = lcm(temp[1][2], temp[2][2]) or lcm(h[2], temp[2][2]) = lcm(temp[1][2], temp[2][2]) then 
Bnew := Bnew, temp; 
end if; 
end do; 
Bnew := [Bnew, op(E)]; 
Bnew := sort(Bnew, proc (x, y) options operator, arrow; x[3] <= y[3] end proc);

Gnew := [op(G),h]; 
RETURN([Gnew, Bnew]); 
end proc;

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



removeing:=proc(B)
RETURN(B[2..-1]);
end proc;





#######################################
###This function, returns set of compatible terms of a polynomial with respect to a vertices set

CompatibleLTs := proc (newpoly, boundarySet, ord, n) 
r := newpoly; 
s := boundarySet; 
U := power(r, n); 
V := NULL; 
w := nops(U); 
 
for i to nops(U) do 
Lt := U[i]; 
P := NULL; 
for j to nops(U) do 
check := U[j]; 
q := Lt-check; 
for k to nops(s) do 
w := s[k]; if 0 < w . q then 
P := P, true; 
break; 
end if; 
end do; 
end do; 
P := [P]; 
if nops(P) = nops(U)-1 then 
V := V, reversePowr(Lt, n); 
end if; 
end do; 
RETURN(([V])); 
end proc;

#######################################
###These two functions, used for selecting a good leading term using Hilbert funciton

compare := proc (s1, i1, s2, i2) 
S := simplify(s1-s2); 
if series(S, pp) = 0 then 
RETURN(0); 
end if; 
k1 := [op(series(S, pp, 0))][2]; 
k2 := series(S, pp, k1+1); 
k := [op(k2)][1]; 
if 0 < k then 
RETURN(1); 
else 
RETURN(2); 
end if; 
end proc;


SelectBestHilbert := proc (G,clt, n) 
g := NULL; 
for i to nops(G) do 
g := g, G[i][2]; 
end do; 
L := NULL; 
for i to nops(clt) do 
F := g; 
F := [F, clt[i]]; 
L := L, [HilbertSeries(F, n, pp), i]; 
end do; 
L := [L]; 
minn := L[1]; 
for i from 2 to nops(L) do 
k := compare(minn[1], minn[2], L[i][1], L[i][2]); 
if k = 1 then 
minn := L[i]; 
end if; 
end do; 
RETURN(clt[minn[2]]); 
end proc;

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

Sort := proc (L, n)
l := add(L); 
r := NULL;
for i to nops(L) do
r := r,LeadingTerm(l,tdeg(op(n)))[2];
l := l - LeadingTerm(l,tdeg(op(n)))[1]* LeadingTerm(l,tdeg(op(n)))[2];
end do;
RETURN([r]); 
end proc;


#######################################
###These two functions, used for counting number of terms in the result


power1 := proc (f, n) 
p := f; 
v := NULL; 
while p <> 0 do 
lm := LeadingMonomial(p, plex(op(n))); 
lt := LeadingTerm(p, plex(op(n)))[1]*LeadingTerm(p, plex(op(n)))[2]; 
u := NULL; 
for i to nops(n) do 
u := u, degree(lm, n[i]); 
end do; 
v := v, [u]; 
p := p-lt; 
end do; 
v := [v]; 
RETURN(v); 
end proc;


termCount := proc(m,n)
tt := NULL; 
for i from 1 to nops(m) do
tt := tt,op(power1(m[i],n));
end do;
RETURN(nops({tt}));
end proc;