weight_over_grevlex_matrix:=function(rows,columns,WEIGHT_MATRIX)

      monomial_order_matrix:=[[0 : j in [1..columns]] : i in [1..columns]];
      for i in [1..rows] do
         for j in [1..columns] do
            monomial_order_matrix[i,j]:=WEIGHT_MATRIX[i,j];
         end for;
      end for;
      for i in [rows+1..columns] do
         for j in [1..columns] do
            if i+j le columns+1 then
               monomial_order_matrix[i,j]:=1;
            end if;
         end for;
      end for;
      return monomial_order_matrix;
   end function;

   grevlex_over_weight_matrix:=function(rows,columns,WEIGHT_MATRIX)

      monomial_order_matrix:=[[0 : j in [1..columns]] : i in [1..columns]];
      for i in [columns-rows+1..columns] do
         for j in [1..columns] do
            monomial_order_matrix[i,j]:=WEIGHT_MATRIX[i-columns+rows,j];
         end for;
      end for;
      for i in [1..columns-rows] do
         for j in [1..columns] do
            if i+j le columns-rows+1 then
               monomial_order_matrix[i,j]:=1;
            end if;
         end for;
      end for;
      return monomial_order_matrix;
   end function;

// changes a rows x columns matrix A into a rows*columns sequence 
// as used in wieghted monmial orderings in MAGMA


   matrix_to_sequence:=function(rows,cols,A)
      sequence_version_of_A:=[0: k in [1..rows*cols]];
      for i in [1..rows] do
         for j in [1..cols] do
            sequence_version_of_A[rows*(i-1)+j]:=A[i,j];
         end for;
      end for;
      return sequence_version_of_A;
   end function;
   
   NF:=function(A,t,I)
      prd:=1;
      temp:=NormalForm(A,I);
      repeat
         rem:=t mod 2;
         t div:=2;
         if rem eq 1 then
            prd *:=temp;
         end if;
         if t ne 0 then
            temp:=NormalForm(temp^2,I);
         end if;
      until t eq 0;
      return prd;
   end function; 





integral_closure:=function(F,S,char,WT_MATRIX,Q_ideal_basis,f)

//F is the finite field
//S is the simple integral extesion of the integrally closed extension Q over P
//char is a (small) positive integer, the size of the finite field F over which all the coefficients are defined
//WT_MATRIX is an independent by var_no array of pseudoweights
//Q_ideal_basis is a sequence of defining relations of the integrally closed extension Q over P
//f is the polynomial defining the integral extension S over Q

// WT_MATRIX_icS is the pseudoweight matrix of the integral closure 
// integral_closure_S is the integral closure ring of S
// I is the ideal of relations of integral_closure
// delta in P is the denominator
// phi is the map from integral_closure_S delta to S
// psi is the map from S to integral_closure_S

// changes a rows x columns matrix WEIGHT_MATRIX 
// to a columns x columns monomial order matrix

   weight_over_grevlex_matrix:=function(rows,columns,WEIGHT_MATRIX)

      monomial_order_matrix:=[[0 : j in [1..columns]] : i in [1..columns]];
      for i in [1..rows] do
         for j in [1..columns] do
            monomial_order_matrix[i,j]:=WEIGHT_MATRIX[i,j];
         end for;
      end for;
      for i in [rows+1..columns] do
         for j in [1..columns] do
            if i+j le columns+1 then
               monomial_order_matrix[i,j]:=1;
            end if;
         end for;
      end for;
      return monomial_order_matrix;
   end function;

   grevlex_over_weight_matrix:=function(rows,columns,WEIGHT_MATRIX)

      monomial_order_matrix:=[[0 : j in [1..columns]] : i in [1..columns]];
      for i in [columns-rows+1..columns] do
         for j in [1..columns] do
            monomial_order_matrix[i,j]:=WEIGHT_MATRIX[i-columns+rows,j];
         end for;
      end for;
      for i in [1..columns-rows] do
         for j in [1..columns] do
            if i+j le columns-rows+1 then
               monomial_order_matrix[i,j]:=1;
            end if;
         end for;
      end for;
      return monomial_order_matrix;
   end function;

// changes a rows x columns matrix A into a rows*columns sequence 
// as used in wieghted monmial orderings in MAGMA


   matrix_to_sequence:=function(rows,cols,A)
      sequence_version_of_A:=[0: k in [1..rows*cols]];
      for i in [1..rows] do
         for j in [1..cols] do
            sequence_version_of_A[rows*(i-1)+j]:=A[i,j];
         end for;
      end for;
      return sequence_version_of_A;
   end function;
   
   NF:=function(A,t,I)
      prd:=1;
      temp:=NormalForm(A,I);
      repeat
         rem:=t mod 2;
         t div:=2;
         if rem eq 1 then
            prd *:=temp;
         end if;
         if t ne 0 then
            temp:=NormalForm(temp^2,I);
         end if;
      until t eq 0;
      return prd;
   end function; 


// produces an integral closure of a simple integral extension
// of an integrally closed quotient ring S/<Q_ideal_basis,f>
// with pseudoweight matrix WT_MATRIX and characteristic q>0


   independent:=#WT_MATRIX;
   var_no:=#WT_MATRIX[1];
   dependent:=var_no-independent-1;
   monomial_order_matrix:=weight_over_grevlex_matrix(independent,var_no,WT_MATRIX);

//  computation of a denominator Delta in P

   S_ideal_basis:=Append(Q_ideal_basis,f);
   I1:=ideal<S|S_ideal_basis>;
   I2:=ideal<S|I1,Minors(JacobianMatrix(S_ideal_basis),dependent+1)>;
   I3:=EliminationIdeal(I2,dependent+1);
   G3:=GroebnerBasis(I3);
   Delta:=GCD(G3);
"Delta=",Delta;

// computing a P-module basis for the initial module M_0
 
   M:=[S.(dependent+1-k): k in [-1..dependent-1]]; M[1]:=1;
   z:=S.1;
   initial_module_size:=Degree(f,z)*(dependent+1);
   M0:= [S|0 : i in [1..initial_module_size]]; 
   for  k in [0..initial_module_size-1] do
      M0[k+1]:=
         z^(k div (dependent+1))*M[1+k-(k div (dependent+1))*(dependent+1)];
   end for;

// initialization

   next_module_basis:=[S|];
   C:=[S|];

//////////////////////////////////////////////////////////////
   INDEX:=function(gg)
      LTgg:=LeadingTerm(gg);
      return [Degree(LTgg,S.ii) : ii in [1..var_no]];
   end function;
//////////////////////////////////////////////////////////////
   IND:=function(gg)
      LTgg:=LeadingTerm(gg);
      return [Degree(LTgg,S.ii): ii in [1..dependent+1]];
   end function;
//////////////////////////////////////////////////////////////
   ind:=function(gg)
      LTgg:=LeadingTerm(gg);
      return [Degree(LTgg,S.i): i in [dependent+2..var_no]];
   end function;
//////////////////////////////////////////////////////////////
   maxi:=function(ind1,ind2)
      return [Max(ind1[i],ind2[i]): i in [1..independent]];
   end function;
//////////////////////////////////////////////////////////////
   dif:=function(ind1,ind2)
      return [ind1[i]-ind2[i]: i in [1..independent]];
   end function;
//////////////////////////////////////////////////////////////
   qthroot:=function(ind1,ind2)
      prd:=S!1;
      for i in [1..independent] do
         if ind1[i] lt ind2[i] then 
            prd:=0; 
            break;
         else
            ii:=ind1[i]-ind2[i];
            jj:=ii div char;
            if jj*char eq ii then 
               prd:=prd*S.(dependent+1+i)^(jj);
            else
               prd:=0;
               break; 
            end if;
         end if;
      end for;
      return prd;
   end function;
//////////////////////////////////////////////////////////////
   gt_equal:=function(ind1,ind2);
      for i in [1..#ind1] do
         if ind1[i] lt ind2[i] then 
            return false; 
         end if;
      end for;
      return true;
   end function;
//////////////////////////////////////////////////////////////
   round_up:=function(ind)
      prod:=S!1;
      for j in [1..independent] do
         sub:=dependent+1+j;
         sup:=(ind[j]+char-1) div char;
         prod:=prod*(S.sub)^sup;
      end for;
      return prod;
   end function; 
//////////////////////////////////////////////////////////////
   inner:=function(ind1,ind2)
      sum:=0;
      for i in [1..var_no] do
         sum:=sum+ind1[i]*ind2[i];
      end for;
      return sum;
   end function;
//////////////////////////////////////////////////////////////
   wt:=function(gg);
      LTgg:=LeadingTerm(gg);
      return [inner(monomial_order_matrix[i],INDEX(LTgg)): i in [1..independent]];
   end function; 
//////////////////////////////////////////////////////////////
   gt_than:=function(ind1,ind2)
      for i in [1..independent] do
         if ind1[i] le ind2[i] then 
            return false; 
         end if;
      end for;
      return true;
   end function;
//////////////////////////////////////////////////////////////
   module_reduction:=function(f,g)
      r:=0;
      LTg:=LeadingTerm(g);
      INDg:=IND(g);
      indg:=ind(g);
      while f ne 0 do
         LTf:=LeadingTerm(f);
         if IND(f) eq INDg and gt_equal(ind(f),indg) then
            f+:=-(LTf div LTg)*g;
         else
            r+:=LTf;
            f+:=-LTf;
         end if;
      end while;
      return r;
   end function;  

//initializing g,G,H/////////////////////////////////////////

   g:= [S|0 : i in [1..initial_module_size]]; 
   G:= [S|0 : i in [1..initial_module_size]]; 
   for i in [1..initial_module_size] do
      g[i]:=M0[i];
      G[i]:=NF(g[i],char,I1);
i,g[i],G[i]; 
   end for;

   for k in [1..initial_module_size] do
      M0[k]:=M0[k]*Delta^(char-1);
   end for;
   H:=[S|0: i in [1..#G]];

//////////////////////////////////////////////////////////////

   basis_is_new:=true;
   while basis_is_new do
      basis_is_new:=false;
      i:=1;

      while i le #g do
"time",i,Cputime(t);
         INDi:=IND(g[i]);
         indi:=ind(g[i]);
//Is g[i] dominated by a basis element?///////////////////////
         skip:=false;
         if #next_module_basis ne 0 then
            for k in [1..#next_module_basis] do
               if INDi eq IND(next_module_basis[k]) then
//                  if qthroot(indi,ind(next_module_basis[k])) ne 0  then
                     skip:=true;
"skip",i,k,wt(g[i]);
                     break;
//                  end if;
               end if;
            end for;
         end if;
//////////////////////////////////////////////////////////////
         if skip then
            Remove(~g,i);
            Remove(~G,i);
            Remove(~H,i);
            i+:=-1;
         else         
//normal form mod old M///////////////////////////////////////
t1:=Cputime();
            G[i]:=NormalForm(G[i],I1);
            j:=#M0;
            while j ne 0 do
               redo:=false;
               if IND(G[i]) eq IND(M0[j]) then
                  if gt_equal(ind(G[i]),ind(M0[j])) then
                     h:=LeadingTerm(G[i]) div LeadingTerm(M0[j]);
                     G[i]+:=-h*M0[j];
                     H[i]+:=h*M0[j];
                     if G[i] eq 0 then 
                        break; 
                     end if;
                     redo:=true;
                     j:=#M0;
                  end if;
               end if;
               if redo eq false then
                  j+:=-1;
               end if;
            end while;
//////////////////////////////////////////////////////////////

//reduction of g,G,H//////////////////////////////////////////
            j:=i-1;
            while j ne 0 and G[i] ne 0 do
               redo:=false;
               if IND(G[i]) eq IND(G[j]) then
                  if gt_equal(ind(G[i]),ind(G[j])) then
                     qij:=qthroot(ind(G[i]),ind(G[j]));
                     if qij ne 0 then
                        lc:=LeadingCoefficient(G[i])/LeadingCoefficient(G[j]);
                        h:=g[i]-lc*g[j]*qij;
                        if h ne 0 and INDEX(g[i]) eq INDEX(h) then 
                           g[i]:=h;
                           qijq:=qij^char;
                           G[i]:=G[i]-lc*G[j]*qijq;
                           H[i]:=H[i]-lc*H[j]*qijq;
                           if G[i] eq 0 then
                              break;
                           end if;
                           G[i]:=NormalForm(G[i],I1);
                           j:=#M0;
                           while j ne 0 do
                              redo:=false;
                              if IND(G[i]) eq IND(M0[j]) then
                                 if gt_equal(ind(G[i]),ind(M0[j])) then
                                    h:=LeadingTerm(G[i]) div LeadingTerm(M0[j]);
                                    G[i]+:=-h*M0[j];
                                    H[i]+:=h*M0[j];
                                    if G[i] eq 0 then 
                                       break; 
                                    end if;
                                    redo:=true;
                                    j:=#M0;
                                 end if;
                              end if;
                              if redo eq false then
                                 j+:=-1;
                              end if;
                           end while;
                           if G[i] eq 0 then 
                              break; 
                           end if;
                           j:=i-1;
                           redo:=true;
                        end if;//INDEX
                     end if;//qij
                  end if;//gt_equal
               end if;//IND
               if redo eq false then
                  j+:=-1;
               end if;
            end while;//j and G[i]
//////////////////////////////////////////////////////////////
i,wt(g[i]),wt(G[i]);

//Is g[i] extendable?/////////////////////////////////////////
            if G[i] ne 0 then
               nw:=true;
               k:=i-1;
               while k ne 0 do
                  if IND(g[k]) eq IND(g[i]) then
                     if LeadingTerm(G[i])*LeadingTerm(g[k])^char eq 
                        LeadingTerm(G[k])*LeadingTerm(g[i])^char 
                        and qthroot(ind(g[i]),ind(g[k])) ne 0 then
                        nw:=false;
                        break;
                     end if;
                  end if;
                  k+:=-1;
               end while;           
            end if;

//////////////////////////////////////////////////////////////
            if G[i] eq 0 then
//basis and removal///////////////////////////////////////////
               r:=false;

//new basis element?//////////////////////////////////////////
               j:=#next_module_basis;
               while j ne 0 do
                  if IND(g[i]) eq IND(next_module_basis[j]) then
                     if gt_equal(ind(g[i]),ind(next_module_basis[j])) then
                        r:=true;
                        break;
                     end if;//gt_equal
                  end if;//IND
                  j+:=-1;
               end while;//j
               if r eq false then
                  Append(~next_module_basis,g[i]);
                  Append(~C,H[i]);
//end new basis element?//////////////////////////////////////

//removal of pending elements dominated by g[i]///////////////
                  j:=i+1;
                  index:=IND(g[i]);
                  indi:=ind(g[i]);
                  while j le #g do
                     if IND(g[j]) eq index then
                        if qthroot(ind(g[j]),indi) ne 0 then
                           Remove(~g,j);
                           Remove(~G,j);
                           Remove(~H,j);
                           j+:=-1;
                        end if; 
                     end if;
                     j+:=1;
                  end while;//j
//////////////////////////////////////////////////////////////
               end if;//r
"REMOVE",wt(g[i]);
               Remove(~g,i);
               Remove(~G,i);
               Remove(~H,i);
               i+:=-1;
//end removal of pending elements dominated by g[i]///////////
            elif nw then
//spoly///////////////////////////////////////////////////////
               j:=i-1;
               while j ne 0 do 
                  if IND(G[i]) eq IND(G[j]) then
                     MAXI:=maxi(ind(G[i]),ind(G[j]));
                     qi:=qthroot(MAXI,ind(G[i]));
                     qj:=qthroot(MAXI,ind(G[j]));
                     if qi ne 0 and qj ne 0 then
                        gi:=g[i]*qi;
                        gj:=g[j]*qj;
                        if gi ge gj then
                           gl:=gi;
                           qiq:=qi^char;
                           Gl:=G[i]*qiq;
                           Hl:=H[i]*qiq;
                        else                    
                           gl:=gj;
                           qjq:=qj^char;
                           Gl:=G[j]*qjq;
                           Hl:=H[j]*qjq;
                        end if;
                        jj:=1;
                        while jj le #g do                
                           if jj eq #g then
                              if LeadingTerm(gl) ne LeadingTerm(g[jj]) then
                                 Append(~g,gl);
                                 Append(~G,Gl);
                                 Append(~H,Hl);
                              end if; 
                              break;
                           elif LeadingTerm(g[jj+1]) gt LeadingTerm(gl) then
                              if LeadingTerm(g[jj]) ne LeadingTerm(gl) then
                                 Insert(~g,jj+1,gl);
                                 Insert(~G,jj+1,Gl);
                                 Insert(~H,jj+1,Hl); 
                              end if; 
                              break;
                           else
                             jj+:=1;
                           end if;//next
                        end while;//jj
                     end if;//qi and qj
                  end if;//IND
                  j+:=-1;
               end while;//j
//end spoly///////////////////////////////////////////////////

//corners/////////////////////////////////////////////////////
               for j in [1..#M0] do
                  if IND(G[i]) eq IND(M0[j]) then
                     MAXI:=maxi(ind(G[i]),ind(M0[j]));
                     Rnd:=round_up(dif(MAXI,ind(G[i])));
                     gl:=g[i]*Rnd;
                     Rndq:=Rnd^char;
                     Gl:=G[i]*Rndq;
                     Hl:=H[i]*Rndq;
                     jj:=1;
                     while jj ne 0 do                
                        if jj eq #g then
                           if LeadingTerm(gl) ne LeadingTerm(g[jj]) then
                              Append(~g,gl);
                              Append(~G,Gl);
                              Append(~H,Hl);                          
                           end if;
                           break;
                        elif LeadingTerm(g[jj+1]) gt LeadingTerm(gl) then
                           if LeadingTerm(g[jj]) ne LeadingTerm(gl) then
                              Insert(~g,jj+1,gl);
                              Insert(~G,jj+1,Gl);
                              Insert(~H,jj+1,Hl);
                           end if;
                           break;
                        else
                           jj+:=1;
                        end if;//next
                     end while;//jj
                  end if;//IND
               end for;//j
//end  corners/////////////////////////////////////////////////
            end if;
         end if;//G[i]
         i+:=1;
      end while;//i

//define new module basis/////////////////////////////////////
      g:=next_module_basis;
      G:=C;
      H:=[S|0: i in [1..#next_module_basis]];
      if #next_module_basis ne #M0 then
         basis_is_new:=true;
      end if;//basis_number
      
      for i in [1..#next_module_basis] do
         if i gt #M0 or M0[i] ne next_module_basis[i]*Delta^(char-1) then
            basis_is_new:=true;
         end if;
      end for;
      if basis_is_new then
         M0:=[next_module_basis[i]*Delta^(char-1): i in [1..#next_module_basis]];
         next_module_basis:=[S|];
         C:=[S|];
      end if; 
//////////////////////////////////////////////////////////////
   end while;//basis_is_new


//////////////////////////////////////////////////////////////

//canonical basis/////////////////////////////////////////////
   gcdb:=GCD(next_module_basis);
   for i in [1..#next_module_basis] do
      next_module_basis[i]:=next_module_basis[i] div gcdb;
i,next_module_basis[i];
   end for;

   for i in [2..#next_module_basis] do
      j:=i-1;
      repeat
         next_module_basis[i]:=module_reduction(next_module_basis[i],next_module_basis[j]);
         j+:=-1;
      until j eq 0;
   end for;

//////////////////////////////////////////////////////////////

   modular_form:=function(g,b);
      r:=[S| 0 : i in [1..#b]];
      j:=#b;
      repeat
          if IND(g) ne IND(b[j])  
          or gt_equal(ind(g),ind(b[j])) eq false then
             j+:=-1;
          else
             lc:=LeadingTerm(g) div LeadingTerm(b[j]);
             g+:=-lc*b[j];
             r[j]+:=lc;
             j:=#b;
          end if;
      until j eq 0;
      return r;
   end function;

//////////////////////////////////////////////////////////////
   delta:=next_module_basis[1];


//ring T//////////////////////////////////////

   var_number:=#next_module_basis-1+independent;
   WT_MATRIX_T:=[[0 : j in [1..var_number]]: i in [1..independent]];
   for i in [1..independent] do
      for j in [1..var_number-independent] do
         WT_MATRIX_T[i,j]:=wt(next_module_basis[#next_module_basis+1-j])[i]-wt(next_module_basis[1])[i];
      end for;
      for j in [1..independent] do
         WT_MATRIX_T[i,j+var_number-independent]:=wt(S.(dependent+1+j))[i];
      end for;
   end for;

   T:=PolynomialRing(F,var_number,"weight",matrix_to_sequence(var_number,var_number,
                     weight_over_grevlex_matrix(independent,var_number,
                     WT_MATRIX_T)));
AssignNames(~T,
["f" cat &cat
[ "_" cat IntegerToString(WT_MATRIX_T[i,j]) 
: i in [1..independent]]
: j in [1..var_number]]);


   image1:=[next_module_basis[#next_module_basis+1-i]: i in [1..#next_module_basis-1]];
   image2:=[S.i*delta: i in [dependent+2..var_no]];
   image:=image1 cat image2;
   phi:=hom<T->S|image>;

   image3:=[T|0 : i in [1..dependent+1]];
   image4:=[T.(var_number-independent+i): i in [1..independent]];
   image5:=image3 cat image4;
   rho:=hom<S->T|image5>;

//////////////////////////////////////////////////////////////


//new ideal I/////////////////////////////////////////////////

   LIN:=[T|];
   for i in [2..#next_module_basis] do
      j:=i-1;
      repeat
         if IND(next_module_basis[i]) eq IND(next_module_basis[j]) then
            gcdij:=GCD(LeadingMonomial(next_module_basis[i]),LeadingMonomial(next_module_basis[j]));
            syzij:=next_module_basis[i]*(LeadingTerm(next_module_basis[j]) div gcdij)
                  -next_module_basis[j]*(LeadingTerm(next_module_basis[i]) div gcdij);
            sij:=modular_form(syzij,next_module_basis);
            lingen:=T.(#next_module_basis+1-i)*(LeadingTerm(next_module_basis[j]) div gcdij)@rho
                 -T.(#next_module_basis+1-j)*(LeadingTerm(next_module_basis[i]) div gcdij)@rho
                 -sij[1]@rho
                 -&+[sij[k]@rho*T.(#next_module_basis+1-k): k in [2..#next_module_basis]];
            Append(~LIN,lingen);
         end if;
         j+:=-1;
      until j eq 0;
   end for;
   ID:=ideal<T|LIN>;
   GID:=GroebnerBasis(ID);

   QUAD:=[T|];
   for i in [2..#next_module_basis] do
      j:=i;
      while j ge 2 do
         sij:=modular_form(NormalForm((next_module_basis[i]*next_module_basis[j]),I1) div delta,next_module_basis);
         quadgen:=T.(#next_module_basis+1-i)*T.(#next_module_basis+1-j)
                 -(sij[1]) @rho
                 -&+[(sij[k])@rho*T.(#next_module_basis+1-k) : k in [2..#next_module_basis]];
         red_quadgen:=NormalForm(quadgen,ID);
         Append(~QUAD,red_quadgen);
         j+:=-1;
      end while;
   end for;

   I_T:=ideal<T|LIN cat QUAD>;

 /////////////////////////////////////////////////////////////

   SS:=[[S|0: i in [1..#next_module_basis]] : i in [1..var_no]];
   for i in [1..var_no] do
      SS[i]:=modular_form(S.i*delta,next_module_basis);
   end for;

   psi:=hom<S->T|[SS[i,1]@rho+&+[SS[i,k]@rho*T.(#next_module_basis+1-k)
                 : k in [2..#next_module_basis]] 
                 : i in [1..var_no]]>;

//////////////////////////////////////////////////////////////
insert:=function(basis,element)
   for i in [1..#basis] do
      if element gt basis[i] then
         return i;
      end if;
   end for;
   return 1+#basis;
end function;
//////////////////////////////////////////////////////////////////
"WT_MATRIX_T=",WT_MATRIX_T;

   repeat
      n:=#WT_MATRIX_T;
      N:=#WT_MATRIX_T[1];
      WT_MATRIX_U:=WT_MATRIX_T;
      N_min:=N; 
//find a pair of columns to switch///////////////////////////
      c1:=WT_MATRIX_T[n,N-n+1];
      row:=0;
      col:=N-n;
      repeat
         c2:=WT_MATRIX_T[1,col];
         Row:=0;
         if c2 lt c1 then
            while row lt n do
               if WT_MATRIX_T[row+1,col] eq c2 then
                  row+:=1;
                  if row eq n then 
                     Row:=n; 
                  end if;
               elif WT_MATRIX_T[row+1,col] eq 0 then 
                  Row:=row;
                  break;
               else
                  break;
               end if;
            end while;
            if Row ne 0 then
//switch is possible////////////////////////////////

               list:=[i:i in [1..N]];
               col2:=N-Row+1;
               wty:=[0: i in [1..n]];
               for i in [1..n] do
                  wty[i]:=WT_MATRIX_T[i,col2];
                  WT_MATRIX_U[i,col2]:=WT_MATRIX_T[i,col];
               end for;
               list[col]:=col2;
               j:=col-1;
               y:=T.col2;
               while j gt 0 do
                  if y gt T.j then
                     for i in [1..n] do
                        WT_MATRIX_U[i,j+1]:=WT_MATRIX_T[i,j];
                     end for;
                     list[j]:=j+1;
                  else
                     for i in [1..n] do
                        WT_MATRIX_U[i,j+1]:=wty[i];
                     end for;
                     list[col2]:=j+1;
                     col3:=j+1; 
                     break;
                  end if;
                  j+:=-1;
               end while;
"WT_MATRIX_U=",WT_MATRIX_U;
//ring U//////////////////////////////////////////////////////

               U:=PolynomialRing(F,N,"weight",matrix_to_sequence(N,N,
                                              weight_over_grevlex_matrix(n,N,
                                              WT_MATRIX_U)));
               var_list:=[U.(list[i]):i in [1..N]];
               homTU:=hom<T->U|var_list>;
               T_list:=[T|0 : j in [1..N]];
               for j in [1..N] do
                  for k in [1..N] do
                     if var_list[k] eq U.j then
                        T_list[j]:=T.k;
                     end if;
                  end for;
               end for;
               homUT:=hom<U->T|T_list>;
               I_U:=ideal<U|[Basis(I_T)[i]@homTU:i in [1..#Basis(I_T)]]>;
               B_U:=Reduce(Basis(I_U));
               list:=[0:j in [1..N]];
               U_list:=[U|];
               prev:=[];
               rev_WT_MATRIX_V:=[[]: i in [1..n]];
               for j in [1..n] do
                  list[j]:=j;
                  Insert(~U_list,j,U.(N+1-j));
                  Insert(~prev,j,j-1-N);
//                  rev_WT_MATRIX_V[j]:=[];
               end for;
               for j in [1..n] do
                  for i in [1..n] do
                     Append(~rev_WT_MATRIX_V[i],WT_MATRIX_U[i,N+1-j]);
                  end for;
               end for;

               l:=n+1;
               for j in [n+1..N] do
                  good:=false;
                  k:=1;
                  repeat
                     if LeadingMonomial(B_U[k]) eq U.(N+1-j) then
                       // Remove(~B_U,k);
                        good:=true;
                     end if;
                     k+:=1;
                  until k gt #B_U;
                  if good eq false then
                     list[j]:=l;
                     for i in [1..n] do
                        Append(~rev_WT_MATRIX_V[i],WT_MATRIX_U[i,N+1-j]);
                     end for;                  
                     Append(~U_list,U.(N+1-j));
                     Append(~prev,j-1-N);
                     l+:=1;
                  end if;
               end for;
               J_U:=ideal<U|B_U>;
               MarkGroebner(J_U);
               j:=n+1;
               repeat
                  temp:=U_list[j]*U.col3;
                  if NormalForm(temp,J_U) eq temp then
                     l:=j+1;
                     while l le #U_list+1 do
                        if l eq #U_list +1 or temp lt U_list[l] then
                           for i in [1..n] do
                              Insert(~rev_WT_MATRIX_V[i],l,rev_WT_MATRIX_V[i,j]+wty[i]);
                           end for;
                           Insert(~U_list,l,temp);
                           Insert(~prev,l,j);
                           for k in [n+1..N] do
                              if list[k] ge l then
                                 list[k]+:=1;
                              end if;
                           end for;
                           j+:=1;
                           l:=j+1;
                           break;
                        else
                           l+:=1;
                        end if;
                     end while;
                  end if;
                  j+:=1;
               until j gt #U_list;
//               k:=1;
//               repeat
//                  if LeadingTotalDegree(B_U[k]) eq 1 then
//                     Remove(~B_U,k);
//                  else
//                     k+:=1;
//                  end if;
//               until k gt #B_U;
//ring V////////////////////////////////////////////////////////////
               N_min:=#U_list;
               WT_MATRIX_V:=[Reverse(rev_WT_MATRIX_V[i]):i in [1..n]];
"WT_MATRIX_V=",WT_MATRIX_V;
               V_sequence:=matrix_to_sequence(N_min,N_min,
                           weight_over_grevlex_matrix(n,N_min,
                           WT_MATRIX_V));
               V:=PolynomialRing(F,N_min,"weight",V_sequence);
               homVU:=hom<V->U|Reverse(U_list)>;
               V_list:=[V|0:j in [1..N]];
               for j in [1..N_min] do
                  if prev[j] lt 0 then 
                     V_list[-prev[j]]:=V.(N_min+1-j);
                  end if;
               end for;
               prehomUV:=hom<U->V|V_list>;
for ii in [1..#V_list] do
   if V_list[ii] eq 0 then

      V_list[ii]:=NormalForm(U.ii,J_U)@prehomUV;
   end if;
end for;
homUV:=hom<U->V|V_list>;
               Y:=(U.col3)@homUV;
               I_V:=ideal<V|[B_U[k]@homUV:k in [1..#B_U]]>;
GroebnerBasis(I_V);
               MarkGroebner(I_V);
               B_V:=[V|Basis(I_V)[k]: k in [1..#Basis(I_V)]];
               for i in [1..N_min] do
                  if prev[i] gt 0 then
                     b_V:=Y*V.(1+N_min-prev[i])-V.(1+N_min-i);
                     ll:=insert(B_V,b_V);
                     Insert(~(B_V),ll,b_V);
                     C_V:=Reduce(B_V); 
                     I_V:=ideal<V|C_V>;
                     MarkGroebner(I_V);
                     B_V:=[Basis(I_V)[k]: k in [1..#Basis(I_V)]];
                  end if;
               end for;
               for i in [n+1..N_min] do
                  for j in [i..N_min] do
                     Vi:=V.(1+N_min-i);
                     Vj:=V.(1+N_min-j);
                     if NormalForm(Vi*Vj,I_V) eq Vi*Vj then

                        bij:=Vi*Vj-NormalForm(NormalForm((Vi@homVU)*(Vj@homVU),I_U)@homUV,I_V);
                        lll:=insert(B_V,bij); 
                        C_V:=[V|0: ll in [1..#B_V+1]];
                        for ll in [1..#B_V+1] do
                           if ll lt lll then
                              C_V[ll]:=B_V[ll];
                           elif ll eq lll then 
                              C_V[ll]:=bij;
                           else
                              C_V[ll]:=B_V[ll-1];
                           end if;
                        end for;
                        I_V:=ideal<V|C_V>;
                        MarkGroebner(I_V);
                        B_V:=[Basis(I_V)[k]: k in [1..#Basis(I_V)]];
                     end if;
                  end for;
               end for;       
//ring W///////////////////////////////////////////////////////////
               WT_MATRIX_W:=WT_MATRIX_V;
               for i in [1..n-1] do
                  for j in [1..N_min] do
                     WT_MATRIX_W[i,j]-:=WT_MATRIX_W[i+1,j];
                 end for;
               end for;
               lcm:=LCM([WT_MATRIX_W[i,N_min+1-i]:i in [1..n]]);
               for i in [1..n] do
                  for j in [1..N_min] do
                     WT_MATRIX_W[i,j] *:=(lcm div WT_MATRIX_W[i,N_min+1-i]);
                  end for;
               end for;
7777777777777777777777777777777777777777777777;
               gcd:=GCD([GCD([WT_MATRIX_W[i,j]:j in [1..N_min]]):i in [1..n]]);
               for i in [1..n] do
                  for j in [1..N_min] do
                     WT_MATRIX_W[i,j] div:=gcd;
                  end for;
               end for;
               for i in [1..n-1] do
                  for j in [1..N_min] do
                     WT_MATRIX_W[n-i,j]+:=WT_MATRIX_W[n-i+1,j];
                  end for;
               end for; 
8888888888888888888888888888888888888888;
"WT_MATRIX_W=",WT_MATRIX_W;
               W:=PolynomialRing(F,N_min,"weight",
                                 matrix_to_sequence(N_min,N_min,
                                  weight_over_grevlex_matrix(n,N_min,
                                  WT_MATRIX_W)));
               homVW:=hom<V->W|[W.j:j in [1..N_min]]>;
               homWV:=hom<W->V|[V.j:j in [1..N_min]]>; 
               I_W:=ideal<W|Basis(I_V)@homVW>;
//ring X/////////////////////////////////////////////////////////
               transpose:=[[WT_MATRIX_W[i,j]: i in [1..n]] : j in [1..N_min]];
               W_list:=[i: i in [1..N_min]];
               for i in [1..N_min-n-1] do
                  temp_row:=transpose[1];
                  temp:=W_list[1];
                  for k in [2..N_min-n] do
                     if transpose[k] gt temp_row then
                        transpose[k-1]:=transpose[k];
                        W_list[k-1]:=W_list[k];
                     else
                        transpose[k-1]:=temp_row;
                        temp_row:=transpose[k];
                       temp:=W_list[k];
                     end if;
                  end for;
                  transpose[N_min-n]:=temp_row;
                  W_list[N_min-n]:=temp;
               end for;
               WT_MATRIX_X:=[[transpose[i,j]:i in [1..N_min]]:j in [1..n]];
               X:=PolynomialRing(F,N_min,"weight",matrix_to_sequence(N_min,N_min,
                                 weight_over_grevlex_matrix(n,N_min,
                                WT_MATRIX_X)));
"WT_MATRIX_X=",WT_MATRIX_X;
               homXW:=hom<X->W|[W.(W_list[i]): i in [1..N_min]]>;
               X_list:=[X| 0 : i in [1..N_min]];
               for k in [1..N_min] do
                  X_list[W_list[k]]:=X.k;
               end for;
               homWX:=hom<W->X|X_list>;
               I_X:=ideal<X|Basis(I_W)@homWX>;
               MarkGroebner(I_X);
AssignNames(~X,
["g" cat &cat
[ "_" cat IntegerToString(WT_MATRIX_X[i,j]) 
: i in [1..n]]
: j in [1..N_min]]);
               Spsi:=[((((S.i@psi)@homTU)@homUV)@homVW)@homWX : i in [1..Rank(S)]];

               N:=N_min;
               WT_MATRIX_T:=WT_MATRIX_X;
               T:=X;
               I_T:=I_X;
               homXT:=hom<X->T|[T.i: i in [1..Rank(X)]]>;
               homTX:=hom<T->X|[X.i: i in [1..Rank(T)]]>; 
               Xphi:=[(((X.i@homXW)@homWV)@homVU)@homUT:i in [1..N_min]];
               XXphi:=[Xphi[i]@phi div delta^(TotalDegree(Xphi[i])-1): i in [1..N_min]];
               phi:=hom<T->S|XXphi>;
               psi:=hom<S->T|Spsi@homXT>;
               break;   
            else
               col+:=-1;

            end if;
         else
            col+:=-1;

         end if;
      until col eq 0;
   until col eq 0;
   Y:=PolynomialRing(F,N_min,"weight",matrix_to_sequence(N_min,N_min,
                     grevlex_over_weight_matrix(independent,N_min,WT_MATRIX_T)));

   homTY:=hom<T->Y|[Y.i: i in [1..Rank(T)]]>;
   I_Y:=ideal<Y|[Basis(I_T)[i]@homTY: i in [1..#Basis(I_T)]]>;
   homYT:=hom<Y->T|[T.i: i in [1..Rank(Y)]]>;
   return WT_MATRIX_T,Y,I_Y,delta,phi,psi;
end function;