with(Groebner): with(PolynomialIdeals):

LM := proc (f, tord)
	options operator, arrow; LeadingMonomial(f, tord); 
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:
 
Between := proc (f, T1, T2) 
local L, N, m, i, flag; 
	L := FL(f, T1); 
	L := sort(L, proc (a, b) options operator, arrow; TestOrder(a, b, T2) end proc); 
	N := NULL; 
	flag := false; 
	m := lt(f, T1); 
	for i from nops(L) by -1 while flag = false do 
		if m = L[i] then 
			N := N, L[i]; 
			flag := true 
		else 
			N := N, L[i] 
		end if; 
	end do; 
RETURN(add(u, u = [N])); 
end proc:

SOrt := proc (p1, p2, tord) 
local u, v; 
global G;
	u := p1[1]; 
	v := p2[1]; 
RETURN(TestOrder(u, v, tord)); 
end proc:
#The next algorithm is the Extended Buchberger's algorithm. 
Stag := proc (FF, tord)
#option trace;
local firsttime, firstbytes, count1, count2, count3, PolList, CritPair, F, g, GG, EB, LCMList, L, TripSet, i, j, k, A, h, t, f1, f2, f, cmp, G, GGG, d; 
	count1 := 0; 
	count2 := 0; 
	count3 := 0; 
	cmp := 0; 
	F := [seq([FF[j], {}, LM(FF[j], tord), j], j = 1 .. nops(FF))]; 
	G := seq([FF[i], [seq(0, l = 1 .. i-1), 1]], i = 1 .. nops(FF)); 
	GGG := NULL; 
	PolList := Array(1 .. nops(F), F); 
	CritPair := Array([]); 
	for i from 2 to nops(F) do 
		LCMList := Array([seq(lcm(PolList[j][3], PolList[i][3]), j = 1 .. i-1)]); 
		L := seq(LCMList[j], j = 1 .. ArrayNumElems(LCMList)); 
		L := Generators(LeadingMonomial(`<,>`(L), tord)); 
		count2 := count2+i-1-nops(L); 
		count3 := count3+i-1; 
		for j in L do 
			member(j, LCMList, 'zz'); 
			CritPair := Array(1 .. ArrayNumElems(CritPair)+1, CritPair); 
			CritPair[-1] := [j, PolList[zz], PolList[i]]; 
		end do;
	end do; 
	while ArrayNumElems(CritPair) <> 0 do 
		CritPair := sort(CritPair, proc (a, b) options operator, arrow; SOrt(a, b, tord) end proc); 
		g := CritPair[1]; 
		CritPair := Array(1 .. ArrayNumElems(CritPair)-1, proc (i) options operator, arrow; CritPair[i+1] end proc); 
		f1 := g[2][1]; f2 := g[3][1]; 
		if type(simplify(f1/LeadingCoefficient(f1, tord)), `+`) or type(simplify(f2/LeadingCoefficient(f2, tord)), `+`) then 
			f := simplify(g[1]*f1/lt(f1, tord)-g[1]*f2/lt(f2, tord)); 
			h := NormalForm(f, [seq(t[1], t = PolList)], tord, 'q'); 
			if h <> 0 then 
				i := g[2][-1]; 
				j := g[3][-1]; 
				G := G, [h, [seq(0, l = 1 .. i-1), g[1]/lt(f1, tord), seq(0, l = i+1 .. ArrayNumElems(PolList))]+[-seq(0, l = 1 .. j-1), -g[1]/lt(f2, tord), -seq(0, l = j+1 .. ArrayNumElems(PolList))]-q]; 
				PolList := Array(1 .. ArrayNumElems(PolList)+1, PolList); 
				PolList[-1] := [h, {}, LM(h, tord), ArrayNumElems(PolList)]; 
				LCMList := Array([seq(lcm(PolList[-1][3], PolList[j][3]), j = 1 .. ArrayNumElems(PolList)-1)]); 
				L := seq(LCMList[j], j = 1 .. ArrayNumElems(LCMList)); 
				L := Generators(LeadingMonomial(`<,>`(L), tord)); 
				count2 := count2+ArrayNumElems(LCMList)-nops(L); 
				count3 := count3+ArrayNumElems(LCMList);
				for j in L do 
					member(j, LCMList, 'zz'); 
					CritPair := Array(1 .. ArrayNumElems(CritPair)+1, CritPair); 
					CritPair[-1] := [j, PolList[zz], PolList[-1]]; 
				end do;
			else 
				i := g[2][-1]; 
				j := g[3][-1]; 
				GGG := GGG, [seq(0, l = 1 .. i-1), g[1]/lt(f1, tord), seq(0, l = i+1 .. ArrayNumElems(PolList))]+[-seq(0, l = 1 .. j-1), -g[1]/lt(f2, tord), -seq(0, l = j+1 .. ArrayNumElems(PolList))]-q; 
				count1 := count1+1; 
			end if; 
		else 
			cmp := cmp+1;
		end if; 
	end do: 
	GG := [seq(PolList[i][1], i = 1 .. ArrayNumElems(PolList))]; 
	lprint("The number of reductions to zero", count1); 
	lprint("The number of removed pairs by Buchberger's criteria", count2); 
	lprint("The number of created critical pairs", count3);
	lprint("The number of removed pairs by new criterion presented in the paper", cmp); 
	lprint("The number of constructed polynomials", ArrayNumElems(PolList)-nops(F)); 
	L := GG; 
	G := [G]; 
RETURN([L, [seq(u[2], u = G),GGG]]); 
end proc:
#The next algorithm is our main algorithm to convert Grobner bases.
Conv := proc (G1, T1, T2)
#option trace;
local F1, F2, G2, A, B, Tr, a, b, G, AA, L, C, TT, c, R, S, SS; 
	G := G1; 
	G2 := op(G); 
	F1 := [seq(Between(f, T1, T2), f = G)]; 
	L := Stag(F1, T2); 
	A := L[1]; 
	B := L[2]; 
	S := HilbertSeries(LeadingMonomial(Basis(G1, T2), T2), {op(T2)}, 't'); 
	while true do 
		G2 := op(G); 
		Tr := NULL; 
		TT := NULL; 
		for b from 1+nops(G) to nops(B) do 
			a := expand(add(B[b][i]*[G2][i], i = 1 .. nops(B[b]))); 
			G2 := G2, a; 
		end do; 
		G2 := op(InterReduce([G2], T2)); 
		SS := HilbertSeries(LeadingMonomial([G2], T2), {op(T2)}, 't'); 
		if S = SS then  
			RETURN([G2]); 
		else 
			F1 := [seq(Between(f, T1, T2), f = [G2])]; 
			AA := A; 
			L := Stag(F1, T2, T1); 
			A := L[1]; 
			B := L[2]; 
			G := [G2]; 
		end if; 
	end do; 
end proc:

###Anna Example 1
F := [x*z+y^2+x, z^2+1]; 
t1 := tdeg(x, y, z); 
t2 := plex(x, y, z); 
F := Basis(F, t1); t := time(); 
A := Conv(F, t1, t2); 
print("The time for the new algorithm is"); 
time()-t; 
t := time(); 
C := Walk(F, t1, t2); 
print("The time for the walk algorithm is"); 
time()-t; 
t := time(); 
Stag(F, t2, t2); 
print("The time for the Berkesch-Schreyer algorithm is"); 
time()-t; 
B := Basis(F, t2); 
evalb(Basis(LeadingMonomial(A, t2), t2) = Basis(LeadingMonomial(B, t2), t2));

###Anna Example 2
F := [x1^3*x2^5*x3+x1^3*x3, x1^7*x2+x1^3+x2*x3*x4, x1^7*x3-x1^3*x2^4*x3-x2^5*x3^2*x4, x2^6*x3+x2*x3, x4^3+x1];
t1 := tdeg(x1, x2, x3, x4); 
t2 := plex(x1, x2, x3, x4); 
F := Basis(F, t1); 
t := time(); 
A := Conv(F, t1, t2); 
print("The time for the new algorithm is"); 
time()-t; 
t := time(); 
C := Walk(F, t1, t2);
print("The time for the walk algorithm is"); 
time()-t; 
t := time(); 
Stag(F, t2, t1):
print("The time for the Berkesch-Schreyer algorithm is"); 
time()-t; 
B := Basis(F, t2); 
evalb(Basis(LeadingMonomial(A, t2), t2) = Basis(LeadingMonomial(B, t2), t2));


###Example 1
F := [y^2*z+x^2, x^3+y^2+z^2]; 
t1 := tdeg(x, y, z); 
t2 := plex(x, y, z); 
F := Basis(F, t1); 
t := time(); 
A := Conv(F, t1, t2); 
print("The time for the new algorithm is"); 
time()-t; t := time(); 
C := Walk(F, t1, t2); 
print("The time for the walk algorithm is"); 
time()-t; 
t := time(); 
Stag(F, t2); 
print("The time for the Berkesch-Schreyer algorithm is"); 
time()-t; 
B := Basis(F, t2);
 evalb(Basis(LeadingMonomial(A, t2), t2) = Basis(LeadingMonomial(B, t2), t2)); 
 
 
 ###Example 2
 F := [-x2^12+x3^12+x3^3*x5^8+x4^8-1, x3^18*x5^4-x4^13*x5^9+x2^8+x2^5*x3^3+x2^4*x3^4-x1+2*x2+1, x4*x5^15-x1*x2+x4+x5+1]; 
 t1 := tdeg(x1, x2, x3, x4, x5); 
 t2 := tdeg(x5, x1, x4, x3, x2); 
 F := Basis(F, t1); 
 t := time(); 
 A := Conv(F, t1, t2); 
 print("The time for the new algorithm is"); 
 time()-t; 
 t := time(); 
 C := Walk(F, t1, t2); 
 print("The time for the walk algorithm is"); 
 time()-t; 
 t := time(); 
 Stag(F, t2); 
 print("The time for the Berkesch-Schreyer algorithm is"); 
 time()-t; 
 B := Basis(F, t2); 
 print("We check the correctness of the new approach:"); 
 evalb(Basis(LeadingMonomial(A, t2), t2) = Basis(LeadingMonomial(B, t2), t2)); 
 
 ###Example 3
 F := [x1^22-x2^22-x1*x2+x2-1, x2^22-x3^22-x1*x2-x3+1, -x1^22+x3^22-x2*x3-x3+1]; 
 t1 := tdeg(x1, x2, x3);
  t2 := tdeg(x3, x2, x1); 
  F := Basis(F, t1); 
  t := time(); 
  A := Conv(F, t1, t2); 
  print("The time for the new algorithm is"); 
  time()-t; 
  t := time(); 
  C := Walk(F, t1, t2); 
  print("The time for the walk algorithm is"); 
  time()-t; 
  t := time(); 
  Stag(F, t1, t2); 
  print("The time for the Berkesch-Schreyer algorithm is"); 
  time()-t; 
  B := Basis(F, t2); 
  print("We check the correctness of the new approach:"); 
  evalb(Basis(LeadingMonomial(A, t2), t2) = Basis(LeadingMonomial(B, t2), t2)); 
  
  ###Example 4
  F := [x1^30-x2^30-x4*x5+x2-1, x1^25-x3^25-x1*x3-x4+1, -x1^25+x3^25-x3*x4-x5^2+x1+1]; 
  t1 := tdeg(x1, x2, x3, x4, x5); 
  t2 := tdeg(x3, x2, x5, x4, x1); 
  F := Basis(F, t1); 
  t := time(); 
  A := Conv(F, t1, t2); 
  print("The time for the new algorithm is"); 
  time()-t; 
  t := time(); 
  C := Walk(F, t1, t2); 
  print("The time for the walk algorithm is"); 
  time()-t; 
  t := time(); 
  Stag(F, t1, t2); 
  print("The time for the Berkesch-Schreyer algorithm is"); 
  time()-t; 
  B := Basis(F, t2); 
  print("We check the correctness of the new approach:"); 
  evalb(Basis(LeadingMonomial(A, t2), t2) = Basis(LeadingMonomial(B, t2), t2)); 
  
  ###Example 5
  F := [x1^22-x2^22-x1*x2+x2-1, x2^22-x3^22-x1-x2+1, x3^22-x4^22-x3^2-x2+1, -x1^22+x4^22-x1+x3-1]; 
  t1 := tdeg(x1, x2, x3, x4); 
  t2 := tdeg(x3, x2, x1, x4); 
  F := Basis(F, t1); 
  t := time(); 
  A := Conv(F, t1, t2); 
  print("The time for the new algorithm is"); 
  time()-t; 
  t := time(); 
  C := Walk(F, t1, t2); 
  print("The time for the walk algorithm is"); 
  time()-t; 
  t := time(); 
  Stag(F, t2); 
  print("The time for the Berkesch-Schreyer algorithm is"); 
  time()-t; 
  B := Basis(F, t2); 
  print("We check the correctness of the new approach:"); 
  evalb(Basis(LeadingMonomial(A, t2), t2) = Basis(LeadingMonomial(B, t2), t2)); 
  
  ###Example 6
  F := [-x2^50+x5^50+x2^10-x2^9*x5-x4^10, -x2^50+x5^50-x1^20-x2^20+x4^10, -x3^50+x5^50-x2^40-x3^40+x1^16-x4^16+x3+x4, -x1^50+x5^50-x1^30-x3^30+x5^2+x1+1];
  t1 := tdeg(x1, x2, x3, x4, x5); 
  t2 := tdeg(x3, x2, x1, x4, x5); 
  F := Basis(F, t1); 
  t := time(); 
  A := Conv(F, t1, t2); 
  print("The time for the new algorithm is"); 
  time()-t; 
  t := time(); 
  C := Walk(F, t1, t2); 
  print("The time for the walk algorithm is"); 
  time()-t; 
  t := time(); 
  Stag(F, t2); 
  print("The time for the Berkesch-Schreyer algorithm is"); 
  time()-t; 
  B := Basis(F, t2); 
  print("We check the correctness of the new approach:"); 
  evalb(Basis(LeadingMonomial(A, t2), t2) = Basis(LeadingMonomial(B, t2), t2)); 
  
  ###Example 7
  F := [x1^5+x2+x3-1, x2^5+x1+x3-1, x3^5+x2+x4-1]; 
  t1 := tdeg(x1, x2, x3, x4); 
  t2 := plex(x1, x2, x3, x4); 
  F := Basis(F, t1); 
  t := time(); 
  A := Conv(F, t1, t2); 
  print("The time for the new algorithm is"); 
  time()-t; 
  t := time(); 
  C := Walk(F, t1, t2); 
  print("The time for the walk algorithm is"); 
  time()-t; 
  t := time(); 
  Stag(F, t2); 
  print("The time for the Berkesch-Schreyer algorithm is"); 
  time()-t; 
  B := Basis(F, t2); 
  print("We check the correctness of the new approach:"); 
  evalb(Basis(LeadingMonomial(A, t2), t2) = Basis(LeadingMonomial(B, t2), t2)); 
  
  ###Example 8
  F := [x1^4-x2^4-x1*x2+1, x2^4-x3^4-x3+1, x3^4-x4^4-x4+2, -x1^4+x4^4-x4+2]; 
  t1 := tdeg(x1, x2, x3, x4); 
  t2 := plex(x1, x2, x3, x4); 
  F := Basis(F, t1); 
  t := time(); 
  A := Conv(F, t1, t2); 
  print("The time for the new algorithm is"); 
  time()-t; 
  t := time(); 
  C := Walk(F, t1, t2); 
  print("The time for the walk algorithm is"); 
  time()-t; 
  t := time(); 
  Stag(F, t2); 
  print("The time for the Berkesch-Schreyer algorithm is"); 
  time()-t; 
  B := Basis(F, t2); 
  print("We check the correctness of the new approach:"); 
  evalb(Basis(LeadingMonomial(A, t2), t2) = Basis(LeadingMonomial(B, t2), t2)); 
  
  ###Example 9
  F := [x2^15-x5^15+x2^2-x2-x4+x5, x1^15-x2^15+x1^10-x2^10+x4+x5, x3^15-x4^15+x5^10+x2+1]; 
  t1 := tdeg(x1, x2, x3, x4, x5); 
  t2 := plex(x1, x2, x3, x4, x5); 
  F := Basis(F, t1); 
  t := time(); 
  A := Conv(F, t1, t2); 
  print("The time for the new algorithm is"); 
  time()-t; 
  t := time(); 
  C := Walk(F, t1, t2); 
  print("The time for the walk algorithm is"); 
  time()-t; 
  t := time(); 
  Stag(F, t2); 
  print("The time for the Berkesch-Schreyer algorithm is"); 
  time()-t; 
  B := Basis(F, t2); 
  print("We check the correctness of the new approach:"); 
  evalb(Basis(LeadingMonomial(A, t2), t2) = Basis(LeadingMonomial(B, t2), t2));