with(Groebner): with(PolynomialIdeals):
deg := proc (u, L)
    RETURN([seq(degree(u, i), `in`(i, L))])
end proc:

cls := proc (u, L)
 local t, i;
    t := deg(u, L):
    for i from 1 to nops(t) do
        if  t[i]<> 0 then
            RETURN(i);
        end if;
    end do;
    RETURN(1);
end proc:

Pommaret:=proc(u,L,G)
local t,A,i;

t:=cls(u,L);

A:=op(L[1..t]):

for i from t+1 to nops(L) do
    if IdealMembership(L[i]*u, <op(G)>) then
        A:=A,L[i];
    fi:
od:
RETURN([A]);
end:

Janet := proc (u, U, L,G)
 local n, A, V,i,A1, t1;
    n := nops(L);
    A:= deg(u, L);
    V := NULL;
    for i from 1 to n-1 do
        A1 := NULL;
        for t1 in U do
            if deg(t1, L)[i+1 .. n] = A[i+1 .. n] then
                A1 := A1,degree(t1, L[i]);
            end if;
        end do;
        if degree(u, L[i]) = max([A1]) or IdealMembership(L[i]*u, <op(G)>) then
            V := V, L[i];
        end if;
    end do;
    if degree(u, L[n]) = max([seq(degree(i,L[n]), `in`(i, U))]) or IdealMembership(L[n]*u, <op(G)>) then
        V := V, L[n];
    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:

NF := proc (f, L, order, L1,G)
local k, Q, r, p, i, flag, U, t1,t2;

k := nops(L);

Q := [seq(0, i = 1 .. k)]; r := 0;

p := NormalForm(f,G,order);
U:=[seq(LeadingMonomial(i,order),`in`(i, L))]; 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)],LeadingMonomial(G,order)))})) then
                Q[i] := 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 flag = false then
   t2 :=LeadingTerm(p, order);
   r := r+t2[1]*t2[2];
   p := simplify(p-t2[1]*t2[2]);
   end if;
   p := NormalForm(p,G,order);
end do;

RETURN(r);
end proc:

NFS := proc (F, order, L1,G) local N, n, i, f,flag;

N :=F;

flag:=false:

while flag=false do

    flag:=true:

    n :=nops(N);
    for i to n do
    f := NF(N[i], subsop(i = NULL, N), order, L1,G);
        if f=0 then
            N:=subsop(i = NULL, N):
            flag:=false;
            break;
        elif f <> 0 and f<>N[i] then
            N[i]:=f:
            flag:=false:
            break;
        end if;
    od:
end do;
RETURN(N);
end proc:


InvAutoreduction := proc (F, order, L1,G)

local N;

N := NFS(F, order, L1,G);

while {op(N)} <> {op(NFS(N, order, L1,G))} do

    N := NFS(N, order, L1,G);

end do;
RETURN(N);
end proc:

Apoly:=proc(FF, order,GG)

local i,j,G,N,F,Q;

N:=NULL:

G:=LeadingMonomial(GG,order);

F:=LeadingMonomial(FF,order);

for i from 1 to nops(F) do

    Q:=Quotient(<op(G)>,<F[i]>):
    Q:=Generators(Q):
    N:=N,seq((FF[i]*q),q=Q):

od:
RETURN([N]):
end:

InvCompletion := proc (FF, order, L1,GG)


local N, n, L, U, A, B, i, j, p,AA,G,F;

G:=Basis(GG,order);

F:=NormalForm(FF,G,order);

N := InvAutoreduction(F, order,L1,G);


n :=nops(N);
L := (op(Apoly(F, order,G)));

U := [seq(LeadingMonomial(i,order),`in`(i,N))];
for i to n do
    A := L1(U[i], U, [op(order)],LeadingMonomial(G,order));
    B := `minus`({op(order)}, {op(A)});
    B := sort([op(B)], proc (a, b) options operator, arrow; TestOrder(a, b, order) end proc);
    for j in B do

            L := L,(N[i]*j);

    end do;
end do;
L := [L];
while nops(L) <> 0 do
    L := sort(L, proc (a, b) options operator, arrow; TestOrder(a,b, order) end proc);
    p := L[1];
    L := L[2 .. -1];
    p := NF(expand(p), N,order, L1,G);
    if p<>0 then
        N := InvAutoreduction([op(N), p], order, L1,G);

        L :=op(Apoly(N, order,G));
        U := [seq(LeadingMonomial(i, order), `in`(i, N))];
        for i to nops(N) do
            A := L1(U[i], U, [op(order)],LeadingMonomial(G,order));
            B := `minus`({op(order)},{op(A)});
            B := sort([op(B)], proc (a, b) options operator, arrow;TestOrder(a, b, order) end proc);
            for j in B do
                if NF(N[i]*j, N,order, L1,G) <> 0 then
                    L := L, (N[i]*j);
                end if;
            end do;
        end do;
        L := [L];
    fi:
end do;
RETURN(N);
end proc:

Test:=proc (AA,FF, order, L1,GG)

local U,i,j,A,B;

U := [seq(LeadingMonomial(i, order), `in`(i, FF))];

for i to nops(FF) do
            A := L1(U[i], U, [op(order)],LeadingMonomial(G,order));
            B := `minus`({op(order)},{op(A)});
            B := sort([op(B)], proc (a, b) options operator, arrow;TestOrder(a, b, order) end proc);
            for j in B do
                if NF(FF[i]*j, FF,order, L1,GG) <> 0 then
                    RETURN(false);
                fi:
            od:
od:

A:=Apoly(FF, order,G):

for j in A do
                if NF(j, FF,order, L1,GG) <> 0 then
                    RETURN(false);
                fi:
od:

if Basis([op(FF),op(GG)],order)<>Basis([op(AA),op(GG)],order) then
    RETURN(false);
fi:
RETURN(true);
 end:

InvCompletion([x1^2*x2-x2^3+x3^3, x1*x2^2-x1^3+x1*x2*x3],
tdeg(x1,x2, x3), Janet, [x1*x2+x2^2]):

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

QuasiStable:= proc (AA, Vars,G)

local B,Var,i,Deg,polys,ee,dd,ListVar,mm,LL,bb,ii,c,VV;

option trace;
    VV:=[seq(Vars[i],i=nops(Vars)..1,-1)]:

    B := Basis(AA,tdeg(op(Vars)));
    B:=LeadingMonomial(NormalForm(B,G,tdeg(op(Vars))),tdeg(op(Vars))):
    Deg:=max(seq(degree(polys), polys in {op(B)} union {op(G)}));
    for bb in B do
        c:=cls(bb,VV):
	    for ii from c+1 to nops(VV) do
	        if not IdealMembership(bb*VV[ii]^Deg/VV[c], <op(B)>) and not IdealMembership(bb*VV[ii]^Deg, <op(LeadingMonomial(G,tdeg(op(Vars))))>) then
		      RETURN(false,VV[c],VV[ii]);
	       fi;
	   od;
    od:
RETURN(true);
 end:


######################## Pommaert Basis #########################

LinearChange := proc (F, Tord,G)

#option trace;

local L, chen, A, u, LL, LLL,i, B,
K,V,uu,m,sj,LR,ltsj,sjj,KK,LH,LS,LLH,VV,GG,LG,LLG;

global tord,T,Q,Vars;

option trace;

GG:=G:

V:=[seq([op(Tord)][i], i = 1 .. nops([op(Tord)]))];


#LH:=InvCompletion(F, Tord,Janet,G);

LH:=Basis(F, Tord);


LG:=G:

chen := NULL;

A := QuasiStable(LH,V,G);

while A <> true do
   	  uu := rand(-5 .. 5)();
      while uu = 0 do
         uu := rand(-5 .. 5)();
      end do;
      m := A[2] = A[2]+uu*A[3];
      LLH:=expand(subs(m, LH));
      LLG:=expand(subs(m, LG));
      #LLH :=InvCompletion(LLH, Tord,Janet,LLG);
      LLH :=Basis(LLH, Tord);
      B := QuasiStable(LLH,V,LLG);;
      if [B] <> [A] then
         LH:=LLH;
         LG:=LLG:
         A := B;
         chen := chen, m;
      end if;
  end do;



return(LH,chen);

end:

LinearChange([x1*x3], tdeg(x3,x2,x1),[x1*x2+x2^2]);