with(Groebner):
deg := proc (u, L) 
	RETURN([seq(degree(u, i), i in L)]): 
end proc:
cls := proc (u, L)
	#option trace:
	local t, i: 
	t := nops(L): 
	for i from 0 to t-1 do 
		if deg(u, L)[t-i] <> 0 then
			RETURN(t-i): 
		end if: 
	end do: 
	RETURN(1): 
end proc:
otherjvar := proc (u, U, L)
	#option trace:
	local n, A, V, i, A1, t1:
	n := nops(L):
	A := deg(u, L): 
	V := NULL: 
	for i from 2 to n do
		A1 := NULL:
		for t1 in U do 
			if deg(t1, L)[1 .. i-1] = A[1 .. i-1] then
				A1 := A1, degree(t1, L[i]):
			end if: 
		end do: 
		if degree(u, L[i]) = max([A1]) then
			V := V, L[i]:
		end if
	end do: 
	if degree(u, L[1]) = max([seq(degree(i, L[1]), i in U)])then 
		V := V, L[1]:
	end if: 
	V := sort([V], proc (a, b) options operator, arrow; TestOrder(a, b, plex(op(L))) end proc): 
	V := [seq(V[i1], i1 = nops(V) .. 1, -1)]:
	RETURN(V) 
end proc:
DNoethervar := proc (u, U, L, d)
	local n, V, V1, t, i: 
	#option trace:
	n := nops(L): 
	V := NULL: 
	V1 := otherjvar(u, U, L): 
	t := cls(u, L): 
	for i to n do 
		if t <= n-d and evalb(L[i] in V1) then
			V := V, L[i]: 
		elif n-d < t and n-d < i and evalb(L[i] in V1) then
			V := V, L[i]: 
		end if: 
	end do: 
     RETURN([V]):
end proc:
NF := proc (f, L, order, L1, d)
	local k, Q, r, p, i, flag, U, t1, t2:
	#option trace: 
	k := ArrayNumElems(L): 
	r := 0: 
	p := f: 
	U := Array([]): 
	for i in L do 
		U(ArrayNumElems(U)+1) := LeadingMonomial(i, order): 
	end do: 
	while p <> 0 do
		i := 1: 
		flag := false: 
		while i <= k and flag = false do
			t1 := [LeadingTerm(L[i], order)]: 
			t2 := [LeadingTerm(p, order)]: 
			if divide(t2[1]*t2[2], t1[1]*t1[2],'t') then
				if evalb(`subset`(indets(t), {op(L1(t1[2], U, [op(order)], d))})) then
					p := simplify(p-t*L[i]): 
					flag := true: 
				else
					i := i+1: 
				end if: 
			else 
				i := i+1: 
			end if: 
		end do: 
		if flag = false then 
			t2 := [LeadingTerm(p, order)]: 
			r := r+t2[1]*t2[2]; p := p-t2[1]*t2[2]: 
		end if: 
	end do: 
	RETURN(r): 
end proc:
TailNF := proc (T, p, order, L1, d)
	local h, G, h1, i:
	#option trace: 
	h := p[1]: 
	G := Array([]): 
	for i in T do 
		G(ArrayNumElems(G)+1) := i[1]: 
	end do: 
	RETURN(NF(h, G, order, L1, d)): 
end proc:
HeadNF := proc (T, p, order, L1, d)
	local h, G, G1, n, i, j, t1, i1, i3, i2, t2: 
	h := p[1]: 
	G := Array([]): 
	for i in T do 
		G(ArrayNumElems(G)+1) := i[1]: 
	end do: 
	lm(h) := LeadingMonomial(h, order): 
	j := false: 
	n := ArrayNumElems(G): 
	G1 := Array([]): 
	for i in G do
		G1(ArrayNumElems(G1)+1) := LeadingMonomial(i, order): 
	end do: 
	t1 := [LeadingTerm(G[1], order)]: 
	if divide(lm(h), t1[1]*t1[2], 't') then
		if evalb(`subset`(indets(lm(h)), {op(L1(t1[2], G1, [op(order)], d))})) then
			j := true: 
		end if: 
	end if: 
	if j and lm(h) <> LeadingMonomial(p[2], order) then
		if critria(p, T[1], order) then
			RETURN(0): 
		end if: 
	end if: 
	while h <> 0 and lm(h) <> NF(lm(h), G, order, L1, d) do
		i3 := 1: 
		i2 := 1: 
		while i2 <= n and i3 = 1 do
			t1 := [LeadingTerm(G[i2], order)]: 
			t2 := [LeadingTerm(h, order)]: 
			if divide(t2[1]*t2[2], t1[1]*t1[2], 't') then
				if evalb(`subset`(indets(t), {op(L1(t1[2], G1, [op(order)], d))})) then
					h := simplify(h-t*G[i2]): 
					i3 := 0: 
					lm(h) := LeadingMonomial(h, order): 
				else 
					i2 := i2+1: 
				end if: 
			else 
			     i2 := i2+1: 
			end if: 
		end do: 
	end do: 
	RETURN(h): 
end proc:
critria := proc (p, g, order)
	#option trace:
	lm(anc(p)) := LeadingMonomial(p[2], order): 
	lm(anc(g)) := LeadingMonomial(g[2], order): 
	lm(pol(p)) := LeadingMonomial(p[1], order): 
	if lm(anc(p))*lm(anc(g)) = lm(pol(p)) or NormalForm(lm(pol(p)), [lcm(lm(anc(p)), lm(anc(g)))], order) = 0 and lm(pol(p)) <> lcm(lm(anc(p)), lm(anc(g))) then
		RETURN(true): 
	else 
		RETURN(false): 
	end if: 
end proc:
HeadReduce := proc (Q, T, order, L1)
	local S, Q1, p, h, S1, S2, i: 
	S := Q: 
	Q1 := Array([]): 
	while ArrayNumElems(S) <> 0 do
		p := S[1]: 
		S := S(2 .. -1): 
		h := HeadNF(T, p, order, L1, d): 
		if h <> 0 then
			lm(pol(p)) := LeadingMonomial(p[1], order): 
			if lm(pol(p)) <> LeadingMonomial(h, order) then
				Q1(ArrayNumElems(Q1)+1) := [h, h, []] 
			else 
				Q1(ArrayNumElems(Q1)+1) := p: 
			end if: 
		else 
			if lm(pol(p)) = LeadingMonomial(p[2], order) then
				S1 := Array([]): 
				for i to ArrayNumElems(S) do
					if S[i][2] = p[1] then
						S1(ArrayNumElems(S1)+1) := S[i]: 
					end if: 
				end do: 
				S2 := Array([]): 
				for i in S do
					if member(i, S1) = false then
						S2(ArrayNumElems(S2)+1) := i: 
					end if: 
				end do: 
				S := S2: 
			end if: 
		end if: 
	end do: 
	RETURN(Q1): 
end proc:
InvolutiveBasis2 := proc (F, order, L1, d)
	local F1, T, Q, p, n, T1, T2, h, U, A, A1, j, t, i: 
	#option trace:
	F1 := sort(F, proc (a, b) options operator, arrow; TestOrder(a, b, order) end proc): 
	T := Array([]): 
	T(ArrayNumElems(T)+1) := [F1[1], F1[1], []]: 
	Q := Array([]): 
	for i from 2 to ArrayNumElems(F1) do
		Q(ArrayNumElems(Q)+1) := [F1[i], F1[i], []]: 
	end do: 
	Q := HeadReduce(Q, T, order, L1): 
	while ArrayNumElems(Q) <> 0 do
		Q := sort(Q, proc (a, b) options operator, arrow; TestOrder(a[1], b[1], order) end proc): 
		p := Q[1]: 
		Q := Q(2 .. -1): 
		if p[1] = p[2] then
			n := ArrayNumElems(T): 
			lm(pol(p)) := LeadingMonomial(p[1], order): 
			T1 := Array([]): 
			for i to n do
				if divide(LeadingMonomial(T[i][1], order), lm(pol(p)), 'q') and q <> 1 then
					Q(ArrayNumElems(Q)+1) := T[i]: 
					T1(ArrayNumElems(T1)+1) := T[i]: 
				end if: 
			end do: 
			T2 := Array([]): 
			for i in T do
				if member(i, T1) = false then
					T2(ArrayNumElems(T2)+1) := i: 
				end if: 
			end do: 
			T := T2: 
		end if: 
		h := TailNF(T, p, order, L1, d): 
		if h <> 0 then 
			T(ArrayNumElems(T)+1) := [h, p[2], p[3]]: 
		end if: 
		n := ArrayNumElems(T): 
		U := [seq(LeadingMonomial(T[i][1], order), i = 1 .. n)]: 
		for i to n do
			A := `minus`({op(order)}, {op(L1(U[i], U, [op(order)], d))}): 
			A1 := `minus`(A, {op(T[i][3])}): 
			A1 := [op(A1)]: 
			if nops(A1) <> 0 then 
				t := nops(A1): 
				for j to t do
					Q(ArrayNumElems(Q)+1) := [T[i][1]*A1[t], T[i][2], []]: 
					T[i][3] := [op(`union`(`intersect`({op(T[i][3])}, A), {A1[t]}))]: 
				end do: 
			end if: 
		end do: 
		Q := HeadReduce(Q, T, order, L1): 
	end do: 
	T2 := Array([]): 
	for i in T do
		T2(ArrayNumElems(T2)+1) := i[1]: 
	end do: 
	RETURN(T2): 
end proc:
NoetherBasisAlgo := proc (F, order)
	local S, H, flag, chen, U, A, c, B, p, h, i, NM, j, H1, F1, d, firsttime, firstbytes, secondbytes, secondtime: 
	firsttime, firstbytes := kernelopts(cputime, bytesused): 
	d := HilbertDimension(F, {op(order)}): 
	F1 := Array([]): 
	for i in F do
		F1(ArrayNumElems(F1)+1) := i: 
	end do: 
	H := InvolutiveBasis2(F1, order, otherjvar, d): 
	chen := NULL: 
	c := 1: 
	U := Array([]): 
	for i in H do
		U(ArrayNumElems(U)+1) := LeadingMonomial(i, order): 
	end do: 
	A := [TestAlgo(U, [op(order)], d)]: 
	while A[1] <> true do
		H1 := subs(A[3] = A[3]+c*A[2], H): 
		H1 := InvolutiveBasis2(H1, order, otherjvar, d): 
		U := Array([]): 
		for i in H1 do 
			U(ArrayNumElems(U)+1) := LeadingMonomial(i, order):
		end do: 
	     B := [TestAlgo(U, [op(order)], d)]: 
	     if B <> A then
	     	chen := chen, A[3] = A[3]+c*A[2]: 
	     	A := B: 
	     	c := 1: 
	     	H := H1: 
	     else c := c+1: 
	     end if: 
	end do: 
	secondtime, secondbytes := kernelopts(cputime, bytesused): 
	RETURN(H, [chen], cputime = secondtime-firsttime, bytesused = secondbytes-firstbytes): 
end proc:
TestAlgo := proc (U, order1, d)
	local u, u1, u2, V, n, i, U1:
	#option trace:
	n := nops(order1): 
	U1 := NULL: 
	for u in U do 
		if n-d < cls(u, order1) then
			U1 := U1, u: 
		end if: 
	end do: 
	U1 := [U1]: 
	for u in U1 do
		u1 := {op(DNoethervar(u, U, order1, d))}: 
		u2 := {op(otherjvar(u, U, order1))}: 
		if u1 <> u2 then 
			u1 := {op(DNoethervar(u, U, order1, d))}: 
			V := `minus`(u2, u1): 
			V := [op(V)]: 
			V := sort(V, proc (a, b) options operator, arrow; TestOrder(a, b, plex(op(order1))) end proc): 
			RETURN(false, V[-1], order1[cls(u, order1)]): 
		end if: 
	end do: 
	RETURN(true): 
end proc:
##Examples##
#print("###########Weispfenning94");
#Weispfenning94 := [h^2*x^2-2*h^2*x*y+h^2*y^2+h^2*z^2+x*y^2*z+y^4, -3*h^5-2*h^2*x^2*y+x*y^4+y*z^4, -2*h^3*x*y+h*x*y^2*z+h*y^4-x^3*y^2+x*y*z^3]:
#NoetherBasisAlgo(Weispfenning94, tdeg(x, y, z, h)):
#print("####################Sturmfels and Eisenbud"):
#Sturmfels and Eisenbud := [ss*uu+bb*vv, tt*uu+bb*ww, tt*vv+ss*ww, ss*xx+bb*yy, tt*xx+bb*zz, tt*yy+ss*zz, vv*xx+uu*yy, ww*xx+uu*zz, ww*yy+vv*zz]:
#NoetherBasisAlgo(Sturmfels and Eisenbud, tdeg(xx, yy, zz, ss, tt, uu, vv, ww,bb)):
#print("####################Liu");
#Liu:=[y*z-y*t0-x*h+a*h, z*t0-z*x-y*h+a*h, t0*x-y*t0-z*h+a*h, x*y-z*x-t0*h+a*h]:
#NoetherBasisAlgo(Liu,tdeg(x,z, y, t0, a, h)):
#print("*******Gerdt2******");
#Gerdt2 := [35*y^4-30*x*y^2-210*y^2*z+3*x^2+30*x*z-105*z^2+140*y*w-21*u, 5*x*y^3-140*y^3*z-3*x^2*y+45*x*y*z-420*y*z^2+210*y^2*w-25*x*w+70*z*w+126*y*u]: 
#Gerdt2 := Homogenize(Gerdt2, h):
#NoetherBasisAlgo(Gerdt2, tdeg(x, y, z, w, u, h)):
#print("****************Vermeer***************");
#Vermeer:= [x^2-2*x*u+u^2+y^2-2*y*v+v^2-1,v^2-u^3,2*v*x-2*v*u+3*u^2*y-3*u^2*v,6*w^2*u^2*v-3*w*u^2-2*w*v+1]:
#NoetherBasisAlgo(Homogenize(Vermeer,h), tdeg(x,y,u,v,w,h)):
#print("####################Cyclic5");
#Cyclic5:=[a*b*c*d*e-1, a*b*c*d + a*b*c*e + a*b*d*e + a*c*d*e +     b*c*d*e, a*b*c + a*b*e + a*d*e + b*c*d + c*d*e, a*b + a*e + b*c +     c*d + d*e, a + b + c + d + e]:
#NoetherBasisAlgo(Homogenize(Cyclic5,h),tdeg(a,b,c,d,e,h)):
#print("####################Eco7");
#Eco7:= [(x1+x1*x2+x2*x3+x3*x4+x4*x5+x5*x6)*x7-1, (x2+x1*x3+x2*x4+x3*x5+x4*x6)*x7-2, (x3+x1*x4+x2*x5+x3*x6)*x7-3, (x4+x1*x5+x2*x6)*x7-4, (x5+x1*x6)*x7-5, x6*x7-6, x1+x2+x3+x4+x5+x6+1]:
#NoetherBasisAlgo(Homogenize(Eco7,h), tdeg(x1, x2, x3, x4,x5,x6,x7,h));