> with(Groebner);
> with(PolynomialIdeals);
################## JanetDivision Algorithm #####################
> Janet := proc (u, U, V) 
  local d, l, i, j, h, A, v, M; 
  d := [seq(max(seq(degree(U[j], V[i]), j = 1 .. nops(U))), i = 1 .. nops(V))]; 
  M := NULL; 
  if degree(u, V[1]) = d[1] then 
         M := M, V[1] 
  end if; 
  for j from 2 to nops(V) do 
         h := [seq(degree(u, V[l]), l = 1 .. j-1)]; 
		 A := NULL; 
		 for v in U do 
		        if [seq(degree(v, V[l]), l = 1 .. j-1)] = h then 
				       A := A, v 
				end if 
		 end do; 
		 if degree(u, V[j]) = max(seq(degree(v, V[j]), v in [A])) then 
		        M := M, V[j] 
		 end if 
  end do; 
  return {M} 
  end proc;
  
################## PommaretDivision Algorithm #####################
> Pommaret := proc (u, U, V) 
  local clasu, Mult, k, i; 
  Mult := NULL; 
  k := indets(u); 
  for i from nops(V) by -1 to 1 while not V[i] in k do 
         Mult := Mult, V[i] 
  end do; 
  Mult := Mult, V[i]; 
  return {Mult} 
  end proc;
  
##################################################################
> VARE := proc (T) 
  local i, M, N; 
  N := [op(T)]; 
  M := NULL; 
  for i to nops(N) while y[i] in [op(T)] do 
         M := M, y[i] 
  end do; 
  for i to nops(N) while x[i] in [op(T)] do 
         M := M, x[i] 
  end do; 
  return [M] 
  end proc;
> LM := proc (f, T) 
  if f <> 0 then 
         RETURN(LeadingMonomial(f, T)) 
  end if; 
  RETURN(0) 
  end proc;
> LT := proc (f, T) 
  if f <> 0 then 
         return LeadingMonomial(f, T)*LeadingCoefficient(f, T) 
  end if; 
  return 0 
  end proc;
  
################## NormalForm Algorithm #####################
> Nor := proc (f, L, T, InDiv) 
  local p, r, m, q, V, Mult, flag, k, i, u, t; 
  p := f; 
  r := 0; 
  m := nops(L); 
  q := [seq(0, i = 1 .. m)]; 
  V := VARE(T); 
  while p <> 0 do 
         i := 1; 
		 flag := false; 
		 while not flag and i <= m do 
		        if divide(LT(p, T), LT(L[i], T), 't') then 
				       Mult := InDiv(LM(L[i], T), [seq(LM(u, T), u = L)], V); 
					   if `subset`(indets(t), {op(Mult)}) then 
					          q[i] := simplify(q[i]+t); 
							  p := simplify(p-t*L[i]); 
							  flag := true 
					   else 
					          i := i+1 
					   end if 
				else 
				       i := i+1 
				end if 
		 end do; 
		 if not flag then 
		        r := simplify(r+LT(p, T)); 
				p := simplify(p-LT(p, T)) 
		 end if 
  end do; 
  return r 
  end proc;
  
################## 	HeadNormalForm Algorithm #####################
> HeadNorm := proc (f, L, T, InDiv) 
  local p, r, m, q, V, Mult, flag, FLAG, k, i, u, t; 
  p := f; 
  r := 0; 
  m := nops(L); 
  q := [seq(0, i = 1 .. m)]; 
  V := VARE(T); 
  FLAG := false; 
  while FLAG = false and p <> 0 do 
         FLAG := true; 
		 i := 1; 
		 flag := false; 
		 while not flag and i <= m do 
		        if divide(LT(p, T), LT(L[i], T), 't') then 
				       Mult := InDiv(LM(L[i], T), [seq(LM(u, T), u = L)], V); 
					   if `subset`(indets(t), {op(Mult)}) then 
					          q[i] := simplify(q[i]+t); 
							  p := simplify(p-t*L[i]); 
							  flag := true; 
							  FLAG := false 
					   else 
					          i := i+1 
					   end if 
				else 
				       i := i+1 
				end if 
		 end do 
  end do; 
  return p 
  end proc;
  
################## HeadAutoreduce Algorithm #####################  
> headautoreduce := proc (P, T, InDiv) 
  local N, L, p, h, K, i; 
  N := P; 
  K := NULL; 
  for i to nops(N) do 
         h := HeadNorm(N[i], [K, op(N[i+1 .. -1])], T, InDiv); 
		 if h <> 0 then 
		        K := K, h 
		 end if 
  end do; 
  return [K] 
  end proc;
  
################## InvolutiveBasis Algorithm #####################  
> INV := proc (F, T, InDiv) 
  local S, NM, G, H, p, i, FLAG, a, j, V, K, h; 
  G := {}; 
  K := headautoreduce(F, T, InDiv); 
  FLAG := false; 
  V := VARE(T); 
  while FLAG = false do 
         FLAG := true; 
		 S := {}; 
		 for j to nops(K) do 
		        NM := `minus`({op(T)}, InDiv(LM(K[j], T), [seq(LM(K[i], T), i = 1 .. nops(K))], V)); 
				for a in NM do 
				       S := `union`(S, {expand(a*K[j])}) 
				end do 
		 end do; 
		 if nops(S) <> 0 then 
		        H := sort(S, proc (a, b) options operator, arrow; TestOrder(a, b, T) end proc); 
				while nops(H) <> 0 do 
				       h := 0; 
					   p := H[1]; 
					   H := H[2 .. -1]; 
					   h := Nor(p, K, T, InDiv); 
					   if h <> 0 then 
					          FLAG := false; 
							  K := [op(K), h]; 
							  K := headautoreduce(K, T, InDiv) 
					   end if 
				end do 
		 end if 
  end do; 
  return K 
  end proc;
> SEP := proc (V, T) 
  local i, n, M; 
  n := nops([op(T)]); 
  M[1] := {}; 
  M[2] := {}; 
  for i to n do 
         if x[i] in V then 
		        M[1] := `union`(M[1], {x[i]}) 
		 end if 
  end do; 
  M[2] := `minus`(V, M[1]); 
  return [M[1], M[2]] 
  end proc;
> TEST := proc (u, F, T) 
  local M, V, N, k, B, Ind; 
  B := VARE(T); 
  M := SEP(Janet(u, F, B), T); 
  N := SEP(Pommaret(u, F, B), T); 
  if nops(M[1]) = nops(N[1]) then 
         return true 
  else 
         V := sort(`minus`(M[1], N[1]), proc (a, b) options operator, arrow; TestOrder(a, b, T) end proc); 
		 Ind := indets(u); 
		 for k from nops(B) by -1 to 1 do 
		        if B[k] in Ind then 
				       return false, V[1], B[k] 
				end if 
		 end do 
  end if 
  end proc;
> STEST := proc (F, T) 
  local u, A; 
  for u in F do 
         A := TEST(u, F, T); 
		 if nops([A]) <> 1 then 
		        RETURN(A) 
		 end if 
  end do; 
  RETURN(true) 
  end proc;
################## XPommaretBasis Algorithm #####################  
> XPOM := proc (F, T, InDiv) 
  local L, chen, A, m, LL, i, B, firsttime, firstbytes, secondtime, secondbytes, K; 
  firsttime, firstbytes := kernelopts(cputime, bytesused); 
  L := INV(F, T, InDiv); 
  L := sort(L, proc (a, b) options operator, arrow; TestOrder(a, b, T) end proc); 
  chen := NULL; 
  A := STEST([seq(LM(L[i], T), i = 1 .. nops(L))], T); 
  while A <> true do 
         m := A[3] = A[3]+(rand(-5 .. 5))()*A[2]; 
		 L := expand(subs(m, L)); 
		 LL := INV(L, T, InDiv); 
		 K := sort(LL, proc (a, b) options operator, arrow; TestOrder(a, b, T) end proc); 
		 B := STEST([seq(LM(K[i], T), i = 1 .. nops(K))], T); 
		 if B <> A then 
		        chen := chen, m; 
				L := K; 
				A := B 
		 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 L, chen 
  end proc;
 ################## BiQuasiStable Algorithm ##################### 
> BQS := proc (L, T) 
  local v, V, F, i, j, l, k, n, m, u, g; 
  F := [seq(LM(L[i], T), i = 1 .. nops(L))]; 
  V := SEP(T); 
  m := nops(V[2]); 
  n := nops(V[1]); 
  for u in F do 
         if degree(u, V[1]) <> 0 then 
		        for i from n by -1 to 1 do 
				       g := degree(u, V[1][i]); 
					   if 0 < g then 
					          k[1] := i; 
							  break 
					   end if 
				end do; 
				for j to k[1]-1 do 
				       v := u*V[1][j]/V[1][k[1]]; 
					   if not v in `<,>`(F) then 
					          return [false, V[1][k[1]], V[1][j]] 
					   end if 
				end do 
		 end if; 
		 if degree(u, V[2]) <> 0 then 
		        for i from m by -1 to 1 do 
				       g := degree(u, V[2][i]); 
					   if 0 < g then 
					          k[2] := i; 
							  break 
					   end if 
				end do; 
				for j to k[2] do 
				       v := u*V[2][j]/V[2][k[2]]; 
					   if not v in `<,>`(F) then 
					          return [false, V[2][k[2]], V[2][j]] 
					   end if 
				end do 
		 end if 
  end do; 
  return true 
  end proc;
> LTransform := proc (F, T) 
  local i, L, A, LL, K, B, chen, m; 
  L := Basis(F, T); 
  L := sort(L, proc (a, b) options operator, arrow; TestOrder(b, a, T) end proc); 
  chen := NULL; 
  A := BQS(L, T); 
  while A <> true do 
         m := A[2] = A[2]+(rand(-2 .. 5))()*A[3]; 
		 L := expand(subs(m, L)); 
		 LL := Basis(L, T); 
		 K := sort(LL, proc (a, b) options operator, arrow; TestOrder(b, a, T) end proc); 
		 B := BQS(K, T); 
		 if B <> A then 
		        chen := chen, m; 
				L := K; 
				A := B 
		 end if 
  end do; 
  return chen, L 
  end proc;