// rocio, last modified 19.06.09
////////////////////////////////////////////////////////////////////////////
version="$Id$";
category="Commutaive algebra";
info="
LIBRARY: binresol.lib    Combinatorial algorithm of resolution of singularities
                        of binomial ideals in arbitrary characteristic.
                        Binomial resolution algorithm of Blanco

AUTHORS:  R. Blanco,            rblanco@modulor.arq.uva.es,
@*        G. Pfister,           pfister@mathematik.uni-kl.de

PROCEDURES:
 Eresol(J);                         computes a E-resolution of singularities of (J) in char 0
 BINresol(J);      computes a E-resolution of singularities of (J) (THE SECOND PART IS NOT IMPLEMENTED YET)
 determinecenter(L1,L2,c,n,Y,a,mb,flag,control3);    computes the next blowing-up center
 Blowupcenter(L1,id,m,L2,c,n,h);    makes the blowing-up
 Nonhyp(Coef,expJ,sJ,n,flag,sums);  computes the ideal generated by the non hyperbolic generators of expJ
 inidata(K,k);                      verifies input data, a binomial ideal K of k generators
 identifyvar();                     identifies status of variables
 data(K,k,n);                       transforms data on lists of lenght n
 Edatalist(Coef,Exp,k,n,flag);      gives the E-order of each term in Exp
 EOrdlist(Coef,Exp,k,n,flag);       computes the E-order of an ideal (giving in the language of lists)
 maxEord(Coef,Exp,k,n,flag);        computes de maximum E-order of an ideal given by Coef and Exp
 ECoef(Coef,expP,sP,V,auxc,n,flag); Computes a simplified version of the E-Coeff ideal.
 elimrep(L);                        removes repeated terms from a list
 Emaxcont(Coef,Exp,k,n,flag);       computes a list of hypersurfaces of E-maximal contact
 cleanunit(mon,n,flag);             clean the units in a monomial mon
 resfunction(t,auxinv,nchart,n);    composes the E-resolution function
 calculateI(Coef,J,c,n,Y,a,b,D);    computes the order of the non monomial part of an ideal J
 Maxord(L,n);                       computes the maximum exponent of an exceptional monomial ideal
 Gamma(L,c,n);                      computes the Gamma function for an exceptional monomial ideal given by L
 convertdata(C,L,n,flag);           computes the ideal corresponding to C,L
 tradblwup(blwhist,n,Y,a,num);      composes the blowing up at this chart
 lcmofall(nchart,mobile);           computes the lcm of the denominators of the E-orders for all the charts
 computemcm(Eolist);                computes the lcm of the denominators of the E-orders for one chart
 constructH(Hhist,n,flag);                construct the list of exceptional divisors accumulated at this chart
 constructblwup(blwhist,n,chy,flag);      construct the ideal defining the composition map
 constructlastblwup(blwhist,n,chy,flag);  construct the ideal defining the last blowup leading to this chart
 genoutput(chart,mobile,nchart,nsons,n,q);             generates the output for visualization
 salida(idchart,chart,mobile,numson,previousa,n,q);    generates the output for one chart
 iniD(n);                           creates a list of lists of zeros of size n
 sumlist(L1,L2);                    sums two lists component to component
 reslist(L1,L2);                    subtracts two lists component to component
 multiplylist(L,a);                 multiplies a list by a number, component to component
 dividelist(L1,L2);                 divides two lists component to component
 createlist(L1,L2);                 creates a list of lists of two elements
 list0(n);                          creates a list of zeros of size n
";

LIB "general.lib";
LIB "qhmoduli.lib";
LIB "inout.lib";
LIB "poly.lib";
LIB "resolve.lib";
LIB "reszeta.lib";
LIB "resgraph.lib";
////////////////////////////////////////////////////////////////////////////

proc inidata(ideal K,int k)
"USAGE: inidata(K,k); K any ideal, k integer (!=0)
COMPUTE: Verifies the input data
RETURN: flag indicating if the ideal is binomial or not
EXAMPLE: example inidata; shows an example
"
{
 int i;
 for (i=1;i<=k; i++)
 {
  if (size(K[i])>2){return(0);}
 }
 return(1);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..3)),dp;
  ideal J1=x(1)^4*x(2)^2, x(1)^2+x(3)^3;
  inidata(J1,2);

  ideal J2=x(1)^4*x(2)^2, x(1)^2+x(2)^3+x(3)^5;
  inidata(J2,2);
}
/////////////////////////////////////////////////////////////////////////////////

proc identifyvar()
"USAGE: identifyvar();
COMPUTE: Asign 0 to variables x and 1 to variables y, only necessary at the beginning
RETURN: list, say l, of size the dimension of the basering
        l[i] is: 0 if the i-th variable is x(i),
                 1 if the i-th variable is y(i)
EXAMPLE: example identifyvar; shows an example
"
{
int i,n;
list flaglist;

n=nvars(basering);
flaglist=list0(n);

for (i=1;i<=n; i++){if (varstr(i)[1]=="y"){flaglist[i]=1;}}

return(flaglist);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp;
  identifyvar();
}
////////////////////////////////////////////////////////////////////////////////

proc data(ideal K,int k,int n)
"USAGE: data(K,k,n); K any ideal, k integer (!=0), n integer (!=0)
COMPUTE: Construcs a list with the coefficients and exponents of one ideal
RETURN: lists of coefficients and exponents of K
EXAMPLE: example data; shows an example
"
{int i,j,lon;
 number aa;
 intvec cc;
 list bb,dd,aux,ddaux,Coef,Exp;

 for (i=1;i<=k; i++)
 { lon=size(K[i]);

// binomial
if (lon==2){aa=leadcoef(K[i][1]);
                   bb=aa;
                   Coef[i]=bb;              // coefficients
                   cc=leadexp(K[i][1]);     // exponents

// cc is an intvec, transform cc in dd, a list of lists
                   dd=cc[1..n];
                   aux[1]=dd;
// the same for the second term

                   aa=leadcoef(K[i][2]);
                   bb=aa;
                   Coef[i]=Coef[i] + bb;  // all the coefficients of i-th generator of K
                   cc=leadexp(K[i][2]);

                   dd=cc[1..n];
                   aux[2]=dd;
                   Exp[i]=aux;}

// monomial
if (lon==1){aux=list();
            aa=leadcoef(K[i][1]);
            bb=aa;
            Coef[i]=bb;
            cc=leadexp(K[i][1]);
            dd=cc[1..n];
            aux[1]=dd;
            Exp[i]=aux;}
} //end for
return(Coef,Exp);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..3)),dp;
  ideal J=x(1)^4*x(2)^2, x(1)^2-x(3)^3;
  data(J,2,3);
}
//////////////////////////////////////////////////////

proc Edatalist(list Coef,list Exp,int k,int n,list flaglist)
"USAGE: Edatalist(Coef,Exp,k,n,flaglist);
        Coef,Exp,flaglist lists, k,n, integers
        Exp is a list of lists of exponents, k=size(Exp)
COMPUTE: computes a list with the E-order of each term
RETURN: a list with the E-order of each term
EXAMPLE: example Edatalist; shows an example
"
{int i,j,lon,mm;
 list dd,ss,sums;
 number aux,aux1,aux2;

 for (i=1;i<=k;i++)
 {
   lon=size(Coef[i]);
   if (lon==1)
   {
     for (j=1;j<=n;j++)
     {
       if (flaglist[j]==0){aux=aux+Exp[i][1][j];}
     }
     ss=aux; aux=0;
   }            // monomial
   else
   {
     for (j=1;j<=n;j++)
     {
       if (flaglist[j]==0)
       {
         aux1=aux1+Exp[i][1][j];
         aux2=aux2+Exp[i][2][j];
       }
     }
     ss=aux1,aux2; aux1=0; aux2=0;
   }   // binomial
   sums[i]=ss;}
   return(sums);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp;
  list flag=identifyvar();
  ideal J=x(1)^3*x(3)-y(2)*y(4)^2,x(5)*y(2)-x(7)*y(4)^2,x(6)^2*(1-y(4)*y(8)^5);
  list L=data(J,3,8);
  list EL=Edatalist(L[1],L[2],3,8,flag);
  EL; // E-order of each term


  ring r = 2,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp;
  list flag=identifyvar();
  ideal J=x(1)^3*x(3)-y(2)*y(4)^2,x(5)*y(2)-x(7)*y(4)^2,x(6)^2*(1-y(4)*y(8)^5);
  list L=data(J,3,8);
  list EL=Edatalist(L[1],L[2],3,8,flag);
  EL; // E-order of each term IN CHAR 2, COMPUTATIONS NEED TO BE DONE IN CHAR 0


  ring r = 0,(x(1..3)),dp;
  list flag=identifyvar();
  ideal J=x(1)^4*x(2)^2, x(1)^2-x(3)^3;
  list L=data(J,2,3);
  list EL=Edatalist(L[1],L[2],2,3,flag);
  EL; // E-order of each term
}
///////////////////////////////////////////////////////////////////////////////////

proc EOrdlist(list Coef,list Exp,int k,int n,list flaglist)
"USAGE: EOrdlist(Coef,Exp,k,n,flaglist);
        Coef,Exp,flaglist lists, k,n, integers
        Exp is a list of lists of exponents, k=size(Exp)
COMPUTE: computes de E-order of an ideal given by a list (Coef,Exp) and extra information
RETURN: maximal E-order, and its position=number of generator and term
EXAMPLE: example EOrdlist; shows an example
"
{int i,can,canpost,lon;
 number canmin;
 list sums;

 sums=Edatalist(Coef,Exp,k,n,flaglist);

 canmin=sums[1][1];                            // inicializating, works also with a monomial
 for (i=1;i<=k; i++)
 {
   lon=size(sums[i]);         // this is 2 for binomial and 1 for monomial generators
   if (sums[i][1]<=canmin and Coef[i][1]!=0)
   {
     canmin=sums[i][1];
     can=i; canpost=1;
   }

// if the generator is a binomial we check the second term

    if (lon==2)
    {
      if (sums[i][2]<canmin and Coef[i][2]!=0)
      {
        canmin=sums[i][2];
        can=i; canpost=2;
      }
    }
  }
  return(canmin,can,canpost);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp;
  list flag=identifyvar();
  ideal J=x(1)^3*x(3)-y(2)*y(4)^2,x(5)*y(2)-x(7)*y(4)^2,x(6)^2*(1-y(4)*y(8)^5),x(7)^4*y(8)^2;
  list L=data(J,4,8);
  list Eo=EOrdlist(L[1],L[2],4,8,flag);
  Eo[1]; // E-order
  Eo[2]; // generator giving the E-order
  Eo[3]; // term giving the E-order
}

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

proc maxEord(list Coef,list Exp,int k,int n,list flaglist)
"USAGE: maxEord(Coef,Exp,k,n,flaglist);
        Coef,Exp,flaglist lists, k,n, integers
        Exp is a list of lists of exponents, k=size(Exp)
RETURN: computes de maximal E-order of an ideal given by Coef,Exp
EXAMPLE: example maxEord; shows an example
"
{
  int i,lon;
  number canmin;  // THE ASSIGNMENT IS NOT OK BECAUSE IT IS OF TYPE NUMBER
  list sums;

  sums=Edatalist(Coef,Exp,k,n,flaglist);

  canmin=sums[1][1];                             // inicializating, works also with a monomial
  for (i=1;i<=k; i++)
  {
    lon=size(sums[i]);         // this is 2 for binomial and 1 for monomial generators
    if (sums[i][1]<=canmin and Coef[i][1]!=0)
    {
      canmin=sums[i][1];
    }

// if the generator is a binomial we check the second term

    if (lon==2)
    {
      if (sums[i][2]<canmin and Coef[i][2]!=0)
      {
        canmin=sums[i][2];
      }
    }
  }
  return(canmin,sums);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp;
  list flag=identifyvar();
  ideal J=x(1)^3*x(3)-y(2)*y(4)^2*x(3),x(5)*y(2)-x(7)*y(4)^2,x(6)^2*(1-y(4)*y(8)^5),x(7)^4*y(8)^2;
  list L=data(J,4,8);
  list M=maxEord(L[1],L[2],4,8,flag);
  M[1];    // E-order
}
//////////////////////////////////////////////////////

proc elimrep(list maxvar)
"USAGE: elimrep(L); L is a list
COMPUTE: Eliminate repeated terms from a list
RETURN: the same list without repeated terms
EXAMPLE: example elimrep; shows an example
"
{
  int i,j;
  list aux2;

  aux2=maxvar;
  for (i=1;i<=size(aux2); i++)
  {
    for (j=i+1;j<=size(aux2); j++)
    {
      if (aux2[i]==aux2[j] and i!=j){aux2=delete(aux2,j);}
    }
  }
  maxvar=aux2;
  return(maxvar);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..3)),dp;
  list L=4,5,2,5,7,8,6,3,2;
  elimrep(L);
}

proc cleanunit(list mon,int n,list flaglist)
"USAGE: cleanunit(mon,n,flaglist);
        mon, flaglist lists, n integer
COMPUTE: We clean (or forget) the units in a monomial, given by "y" variables
RETURN: The list defining the monomial ideal already cleaned
EXAMPLE: example cleanunit; shows an example
"
{
  int i;
  for (i=1;i<=n;i++)
  {
    if (flaglist[i]==1){mon[i]=0;}
  }
  // coef[1]=coef[1]*y(i)^mon[i]; IS NOT ALLOWED because mon[i] can be a number
  // therefore, the coefficients remain constant
  return(mon);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1),y(2),x(3),y(4)),dp;
  list flag=identifyvar();
  ideal J=x(1)^3*y(2)*x(3)^5*y(4)^8;
  list L=data(J,1,4);
  L[2][1][1];  // list of exponents of the monomial J
  list M=cleanunit(L[2][1][1],4,flag);
  M;           // new list without units
}
//////////////////////////////////////////////////////
// Classification of the ideal E-Coeff_V(P):
// ccase=1,    E-Coeff_V(P)=0
//       2,3   Bold regular case
//       4     P=1 monomial case (detected before)
//       0     Otherwise

proc ECoef(list Coef,list expP,int sP,int V,number auxc,int n,list flaglist)
"USAGE: ECoef(Coef,expP,sP,V,auxc,n,flaglist);
        Coef, expP, flaglist lists, sP, V, n integers, auxc number
COMPUTE: The ideal E-Coeff_V(P), where V is a permissible hypersurface which belongs to the center
RETURN: list of exponents, list of coefficients and classification of the ideal E-Coeff_V(P)
EXAMPLE: example ECoef; shows an example
"
{
  int i,j,k,l,numg,ccase,cont2,cont3,val;
  number aa;
  list Eco,newcoef,auxexp,newL,rs,rs2,aux,aux2,aux3,aux4,L;

  auxexp=expP;

  l=1;
  for (i=1;i<=sP;i++)
  {
    rs[i]=size(Coef[i]);
    if (rs[i]==2)
    {                                   // binomials
      if (auxexp[i][1][V]!=auxexp[i][2][V])  // no common factors for the variable in V

      {
        for (j=1;j<=2;j++)
        {
          if (auxexp[i][j][V]<auxc)
          {
            aa=auxc/(auxc-auxexp[i][j][V]);
            auxexp[i][j][V]=0;
            aux4[1]=multiplylist(auxexp[i][j],aa);
            Eco[l]=aux4;
            // newcoef[l]=Coef[i][j]^aa; IT IS NO ALLOWED!!!
            newcoef[l]=Coef[i][j]; // we leave it constant
            l=l+1;
          }
        }
      }
      else                               // common factors for the variable in V, of zero in both terms

      {
        if (auxexp[i][1][V]<auxc)
        {
          aa=auxc/(auxc-auxexp[i][1][V]);
          auxexp[i][1][V]=0; auxexp[i][2][V]=0;

          // this generator is a power of a binomial
          // one possibility is Eco[l]=auxexp[i]; we leave it constant and add some extra number aa, or
          // define a binomial again. The E-order coincides!!!

          aux=multiplylist(auxexp[i][1],aa);
          aux2=multiplylist(auxexp[i][2],aa);
          aux3[1]=aux;
          aux3[2]=aux2;
          Eco[l]=aux3;                                                                    newcoef[l]=Coef[i];
          l=l+1;
        }
      }
    }
    else                                               // monomials
    {
      if (auxexp[i][1][V]<auxc)
      {
        aa=auxc/(auxc-auxexp[i][1][V]);
        auxexp[i][1][V]=0;
        aux4=list();
        aux4[1]=multiplylist(auxexp[i][1],aa);
        Eco[l]=aux4;
        newcoef[l]=Coef[i];
        l=l+1;
      }
    }
  }

// cleaning units from the monomial generators of Eco
// If there are hyperbolic equations in Eco, such that Eco=1, we detect it later, computing the E-order

  numg=size(Eco);
  for (k=1;k<=numg;k++)
  {
    if (size(newcoef[k])==1){Eco[k][1]=cleanunit(Eco[k][1],n,flaglist);}
  }

// checking Eco

  ccase=0;
  cont2=0;
  cont3=0;
  val=0;

// CASE Eco=0: If Eco=empty list then as ideal Eco=0

  if (numg==0){ccase=1;}
  else
  {
    for (i=1;i<=numg;i++)
    {
      rs2[i]=size(newcoef[i]);
      if (rs2[i]==1)
      {
        val=val+n;                                               // monomials
        for (l=1;l<=n; l++)
        {
          if (Eco[i][1][l]==0) {cont2=cont2+1;}
        }
      }
      else
      {
        val=val+(2*n);                                                     // binomials
        for (l=1;l<=n; l++)
        {
          if (Eco[i][1][l]==0) {cont2=cont2+1;}
          if (Eco[i][2][l]==0) {cont2=cont2+1;}
        }
      }
    }

// If cont2=val then all the entries of Eco are zero!! As ideal Eco=1

    for (i=1;i<=sP;i++)
    {
      if (rs[i]==2)
      {                                                           // binomials
        for (l=1;l<=n;l++)
        {
          if (expP[i][1][l]!=0) {cont3=cont3+1;}
          if (expP[i][2][l]!=0) {cont3=cont3+1;}
        }
      }
      else
      {                                                                    // monomials
        for (l=1;l<=n;l++)
        {
          if (expP[i][1][l]!=0)
          {
            cont3=cont3+1;
           }
         }
       }
     }

// If cont3=0 all the entries of expP are zero!! As ideal P=1 this is detected before
// If cont3=1 then P is bold regular


// CASE Eco=1

     if (cont2==val and cont3==1){ccase=2;}    // BOLD REGULAR CASE
     if (cont2==val and cont3>1){ccase=3;}     // CASE P=x^{\alpha},x^{\beta}, IN FACT, BOLD REGULAR
     if (cont2==val and cont3==0){ccase=4;}    // P=1, then I=1 monomial case

// Case BOLD REGULAR P=x^{\alpha}*(1-\mu y^{\delta})
// IT IS NON NECESSARY TO CHECK IT, Eco=empty list, already done!

     L=maxEord(newcoef,Eco,numg,n,flaglist);         // L[1] is the E-order of Eco
     if (L[1]==0){ccase=2; print("E-order zero!");}  // BOLD REGULAR CASE

// we leave it to check the computations

  } // close else

  return(Eco,newcoef,ccase);
}
example
{"EXAMPLE:"; echo = 2;
ring r = 0,(x(1),y(2),x(3),y(4),x(5..7)),dp;
list flag=identifyvar();
ideal P=x(1)^2*x(3)^5-x(5)^7*y(4),x(6)^3*y(2)^5-x(7)^5,x(5)^3*x(6)-y(4)^3*x(1)^5;
list L=data(P,3,7);
list L2=ECoef(L[1],L[2],3,1,3,7,flag);
L2[1];  // exponents of the E-Coefficient ideal respect to x(1)
L2[2];  // its coefficients
L2[3];  // classify the type of ideal obtained

ring r = 0,(x(1),y(2),x(3),y(4)),dp;
list flag=identifyvar();
ideal J=x(1)^3*(1-2*y(2)*y(4)^2);  // Bold regular case
list L=data(J,1,4);
list L2=ECoef(L[1],L[2],1,1,3,4,flag);
L2;

ring r = 0,(x(1),y(2),x(3),y(4),x(5..7)),dp;
list flag=identifyvar();
ideal J=x(1)^3-x(3)^2*y(4)^2,x(1)*x(7)*y(2)-x(6)^3*x(5)*y(4)^3,x(5)^3-x(5)^3*y(2)^2;
list L=data(J,3,7);
list L2=ECoef(L[1],L[2],3,1,2,7,flag);
L2;

ring r = 3,(x(1),y(2),x(3),y(4),x(5..7)),dp;
list flag=identifyvar();
ideal J=x(1)^3-x(3)^2*y(4)^2,x(1)*x(7)*y(2)-x(6)^3*x(5)*y(4)^3,x(5)^3-x(5)^3*y(2)^2;
list L=data(J,3,7);
list L2=ECoef(L[1],L[2],3,1,2,7,flag);
L2; // THE COMPUTATIONS ARE NOT CORRECT IN CHARACTERISTIC p>0
    // because numbers are treated as 0 in assignments

}
////////////////////////////////////////////////////////////////////////////
// The intvec a indicates the previous center
// Hhist = intvec of exceptional divisors of the parent chart

proc determinecenter(list Coef,list expJ,number c,int n,int Y,intvec a,list listmb,list flag,int control3,intvec Hhist)
"USAGE: determinecenter(Coef,expJ,c,n,Y,a,listmb,flag,control3,Hhist);
        Coef, expJ, listmb, flag lists, c number, n, Y, control3 integers, a, Hhist intvec
COMPUTE: next center of blowing up and related information, see example
RETURN: several lists defining the center and related information
EXAMPLE: example determinecenter; shows an example
"
{
  int i,j,rstep,l,mm,cont,cont1,cont2,cont3,a4,sI,sP,V,V2,ccase,b,Mindx,tip;
   number auxc,a1,a2,ex,maxEo,aux;

   list D,H,auxJ; // lists of D_n,D_n-1,...,D_1; H_n,H_n-1,...,H_1; J_n,J_n-1,...,J_1

   list oldOlist,oldC,oldt,oldD,oldH,allH;  // information of the previous step

   list Olist,C,t,Dstar,center,expI,expP,newJ,maxset;

   list maxvar,auxlist,aux3,auxD,auxolist,auxdiv,auxaux,L,rs,auxgamma,auxg2,aux1;   // auxiliary lists
   list auxinvlist,newcoef,EL,Ecoaux,Hplus,transH,Hsum,auxset,sumnewH;              // auxiliary lists
   list auxcoefI;

   intvec oldinfobo7,infobo7;
   int infaux,leh,leh2,leh3;

  tip=listmb[1];   // It is not used in this procedure, it is used to compute the lcm of the denominators
  oldOlist=listmb[2];
  oldC=listmb[3];
  oldt=listmb[4];  // t= resolution function
  oldD=listmb[5];

  oldH=listmb[6];
  allH=listmb[7];

  oldinfobo7=listmb[8];  // auxiliary intvec, it is used to define BO[7]

  // inicializating lists
   Olist=list();
   C=list();
   auxinvlist=list();

   auxJ[1]=expJ;
   rstep=n;             // we are in dimension rstep
   auxc=c;
   cont=1;

  if (Y==0) {D=iniD(n); H=iniD(n); infobo7=-1;} // first center, inicializate previous information

   if (Y!=0 and rstep==n)           // In dimension n, D'_n is always of this form
     { auxdiv=list0(n);
       Dstar[1]=oldD[1];

       b=size(a);
       for (i=1;i<=n;i++) {for (j=1;j<=b;j++) {if (a[j]==i) {aux=aux+oldD[1][i];}}}
       Dstar[1][Y]=aux;
       aux=0;

       auxdiv[Y]=oldOlist[1]-oldC[1];
       D[1]=sumlist(Dstar[1],auxdiv);}  // list defining D_n

  // computing strict transforms of the exceptional divisors H

  if (Y!=0){transH=iniD(n);
            for (i=1;i<=size(oldH);i++){transH[i]=oldH[i]; transH[i][Y]=0;}         // Note: size(oldH)<=n
            allH[Y]=1;}                                                             // transform of |H|=H_nU...UH_1

  // We put here size(oldH) instead of n because maybe we have not
  // calculated all the dimensions in the previous step

  // STARTING THE LOOP

   while (rstep>=1)
    {
      if (Y!=0 and rstep!=n) // transformation law of D_i for i<n
      {
        if (cont!=0)      // the resolution function did not drop in higher dimensions
        {
         if (oldt[n-rstep]==a1/a2 and c==oldC[1] and control3==0)
          {auxD=list0(n);
           auxD[Y]=oldOlist[n-rstep+1]-oldC[n-rstep+1];
            Dstar[n-rstep+1]=oldD[n-rstep+1];

             for (i=1;i<=n;i++) {for (j=1;j<=b;j++) {if (a[j]==i) {aux=aux+oldD[n-rstep+1][i];}}}
             Dstar[n-rstep+1][Y]=aux;
             aux=0;

             D[n-rstep+1]=sumlist(Dstar[n-rstep+1],auxD);

          }
         else
             {cont=0;
              for (j=n-rstep+1;j<=n; j++){D[j]=list0(n);}
             }
        }
      }

  // Factorizing J=M*I

     cont1=0;
     for (i=1;i<=n;i++) {if (D[n-rstep+1][i]==0) {cont1=cont1+1;}}  // if it fails write: listO(n)[i]

     if (cont1==n) {expI=expJ;}    // D[n-rstep+1]=0 (is a list of zeros)
     else {
           for (i=1;i<=size(expJ);i++)
           {rs[i]=size(Coef[i]);
            if (rs[i]==2){ aux1=list();
                           aux1[1]=reslist(expJ[i][1],D[n-rstep+1]);
                           aux1[2]=reslist(expJ[i][2],D[n-rstep+1]);
                           expI[i]=aux1;}                             // binomial
            else {aux1=list();
                  aux1[1]=reslist(expJ[i][1],D[n-rstep+1]);
                  expI[i]=aux1;}}                                     // monomial
          }

  // NOTE: coeficients of I = coeficients of J, because I and J differ in a monomial

  // Detecting errors, negative exponents in expI

  sI=size(expI);

  for (i=1;i<=sI;i++)
  {
    rs[i]=size(Coef[i]);
    if (rs[i]==2)
    {
      for (j=1;j<=2;j++)i
      {
        for (l=1;l<=n; l++)
        {
          if (expI[i][j][l]<0)
          {
            ERROR("the BINOMIAL ideal I has negative components");
          }
        }
      }
    }
    else
    {
      for (l=1;l<=n; l++)
      {
        if (expI[i][1][l]<0)
        {
          "the MONOMIAL ideal I has negative components";
          "M ideal";
          print(D[n-rstep+1]); print(expI);
          "dimension"; print(rstep);
        }
      }
    }
  }

  // Compute the maximal E-order of I

  L=maxEord(Coef,expI,sI,n,flag);
  maxEo=L[1]; // E-order of I

  // Inicializating information

     auxolist=maxEo;
     a1=maxEo;
     a2=auxc;
     Olist=Olist+auxolist;  // list of new maximal E-orders o_n,o_{n-1},...o_1
     aux3=auxc;
     C=C+aux3;              // list of new critical values c=c_{n+1},c_{n},...c_2

  // It is necessary to check if the first coordinate of the invariant has dropped or not
  // NOTE: By construction, the first coordinate is always 1 !!
  // It has dropped is equivalent to: CURRENT C<c of the previous step

  // Calculate new H, this is done for every dimension

  if (Y!=0){a4=size(oldt);
            if (n-rstep+1>a4){cont=0; oldt[n-rstep+1]=0; }            // VERIFICAR!!!!

            if (cont!=0 and oldt[n-rstep+1]==a1/a2 and c==oldC[1] and control3==0){H[n-rstep+1]=transH[n-rstep+1];

                // we fill now the value for BO[7]
                      if (oldinfobo7[n-rstep+1]==-1){leh=size(Hhist);
                                                     infobo7[n-rstep+1]=Hhist[leh];} // suitable index !!!
                      else{ infaux=oldinfobo7[n-rstep+1];
                            infobo7[n-rstep+1]=infaux;}      // the same as the previous step

                                                                                  }
            else {
                  if (rstep<n) {sumnewH=list0(n);
                                for (i=1;i<n-rstep+1;i++){sumnewH=sumlist(sumnewH,H[i]);}
                                H[n-rstep+1]=reslist(allH,sumnewH);}
                  else {H[n-rstep+1]=allH;}

                   // we fill the value for BO[7] too, we complete it at the end if necessary
                  infobo7[n-rstep+1]=-1;
                 }
            }

  // It is necessary to detect the monomial case AFTER inicializate the information
  // OTHERWISE WE WILL HAVE EMPTY COMPONENTS IN THE RESOLUTION FUNCTION

  // If maxEo=0 but maxo!=0 MONOMIAL CASE (because E-Sing(J,c) still !=emptyset)
  // If maxEo=0 and maxo=0 then I=1, (real) monomial case, the same case for us
  // NOTE THAT IT DOESN'T MATTER IF THERE IS A p-TH POWER OF A HYPERBOLIC EQ, THE E-ORDER IS ZERO ANYWAY

  if (maxEo==0){auxgamma=Gamma(D[n-rstep+1],auxc,n); // Gamma gives (maxlist,gamma,center)
                auxg2=auxgamma[3];
                center=center+auxg2;
                center=elimrep(center);
                auxinvlist=auxgamma[2]; break;}

  // Calculate P    // P=I+M^{o/(c-o)} with weight o

  if (maxEo>=auxc) {expP=expI; Mindx=0;}                 // The coefficients also remain constant
     else {ex=maxEo/(auxc-maxEo);
           auxlist=list();
           Mindx=1;
           auxlist[1]=multiplylist(D[n-rstep+1],ex);     // weighted monomial part: D[n-rstep+1]^ex;
           expP=insert(expI,auxlist);                    // P=I+D[n-rstep+1]^ex;
           auxcoefI=Coef;
           Coef=insert(Coef,list(1));}                   // Adding the coefficient for M

  // NOTE: IT IS NECESSARY TO ADD COEFFICIENT 1 TO THE MONOMIAL PART M
  // E-ord(P_i)=E-ord(I_i) so to compute the E-order of P_i we can compute E-ord(I_i)

  // Calculate variables of E-maximal contact, ALWAYS WITH RESPECT TO THE IDEAL I !!

  sP=size(expP);     // Can be different from size(expI)

  if (Mindx==1){ maxvar=Emaxcont(auxcoefI,expI,sI,n,flag);}
  else{ maxvar=Emaxcont(Coef,expP,sP,n,flag);}

  auxc=maxvar[1];     // E-order of P, critical value for the next step, ALSO VALID auxc=maxEo;
  if (auxc!=maxEo)
  {
    ERROR("the E-order is not well computed");
  }

  maxset=maxvar[2];
  center=center + maxset;

  // Cleaning the center: eliminating repeated variables

  center=elimrep(center);

  if (rstep==1) {break;}   // Induction finished, is not necessary to compute the rest

  // Calculate Hplus=set of non permissible hypersurfaces
  // RESET Hplus if c has dropped or we have eliminated hyperbolic generators

  // ES NECESARIO PONER CONDICION DE SI INVARIANTE BAJA O NO??? SI BAJA HPLUS NO SE USA...

  if (Y==0 or c<oldC[1] or control3==1) {Hplus=list0(n);}
  else {Hsum=list0(n);
        Hplus=allH;
        for (i=1;i<=n-rstep+1;i++){Hsum=sumlist(Hsum,H[i]);}
        Hplus=reslist(Hplus,Hsum);                             // CHEQUEAR QUE NO SALEN -1'S
       }

  // Taking into account variables of maxset outside of Hplus (so inside Hminus)

  if (Y==0 or c<oldC[1] or control3==1){V=maxset[1];                  // Hplus=0 so any variable is permissible
                                        maxset=delete(maxset,1);}     // eliminating this variable V from maxset
  else{
       // If the invariant remains constant V comes from the previous step

      if (cont!=0 and oldt[n-rstep+1]==a1/a2 and c==oldC[1]){
                                                             if (Mindx==1){
  //----------------------------USING HPLUS----------------------------------------
  // REMIND THAT IN THIS CASE maxset=HYPERSURFACES OF E-MAXIMAL CONTACT FOR I, INSTEAD OF P

                                                        cont2=1;
                                                        cont3=1;
                                                        auxset=maxset;
                                                          while (cont2!=0){mm=auxset[1];
                                                            if (Hplus[mm]!=0) {auxset=delete(auxset,1); cont3=cont3+1;}
                                                                 // eliminating non permissible variables from maxset
                                                            else {cont2=0;}}
                                                        V=maxset[cont3];        // first permissible variable
                                                        maxset=delete(maxset,cont3);
                          V2=a[n-rstep+1];        // can be different from the variable coming from the previous step
                                                                           }

  //-------------------------------------------------------------------------------
                                                             else{ V=a[n-rstep+1];}
                                                            }
      else {V=maxset[1];     // Hplus=0 so any variable is permissible
            maxset=delete(maxset,1);
           }

       }

  // Calculate Eco=E-Coeff_V(P) where V is a permissible hypersurface which belongs to the center
  // Eco can have rational exponents

    Ecoaux=ECoef(Coef,expP,sP,V,auxc,n,flag);

  // SPECIAL CASES: BOLD REGULAR CASE
  //--------------------------------------------------------------------

    if (Ecoaux[3]==1)
    {                      // Eco=EMPTY LIST, Eco=0 AS IDEAL
      aux1[1]=list0(n);
      newJ[1]=aux1;         // monomial with zero entries, newJ=1 as ideal
      newcoef[1]=list(1);   // the new coefficient is only 1
      auxaux=list();
      auxaux[1]=newJ;
      auxJ=auxJ+auxaux;     // auxJ list of ideals J_i
      auxinvlist=1;
      break;
    }
  //-----------------------------------------------------------
  // THIS CASE IS NOT GOING TO APPEAR, BUT WE LEAVE IT TO CHECK COMPUTATIONS
    if (Ecoaux[3]==2 or Ecoaux[3]==3)
    {                     // Eco=0 LIST, Eco=1 AS IDEAL
      aux1[1]=list0(n);
      newJ[1]=aux1;
      newcoef[1]=list(1);  print("Strange case happens");
      auxaux=list();
      auxaux[1]=newJ;
      auxJ=auxJ + auxaux;  // auxJ list of ideals J_i
      auxinvlist=1;
      break;
    }
  //-----------------------------------------------------------
  // THIS CASE IS NOT GOING TO APPEAR, BUT WE LEAVE IT TO CHECK COMPUTATIONS

  // P=1 THIS CANNOT HAPPEN SINCE P=1 IFF I=1 (or I is equivalent to 1)
  // and this is the monomial case, already checked

    if (Ecoaux[3]==4){ERROR("ERROR in ECoef");}
  //-----------------------------------------------------------
  // If we are here Ecoaux[3]=0, then continue
  // Filling the list of "ideals", auxJ=J_n,J_{n-1},...
    newJ=Ecoaux[1];
    newcoef=Ecoaux[2];

    auxJ=insert(auxJ,newJ,n-rstep+1); // newJ is inserted after n-rstep+1 position, so in position n-rstep+2

  // New input for the loop, if we are here newJ is different from 0

    expJ=newJ;
    Coef=newcoef;

    newJ=list();
    expI=list();
    expP=list();
    rstep=rstep-1;  // print(size(auxJ));
  }
  // EXIT LOOP "while"
  // we do NOT construct the center as an ideal because WE USE LISTS
  t=dividelist(Olist,C);    // resolution function t
  // Complete the intvec infobo7 if necessary
  if (control3==1){infobo7=-1;} // We reset the value after clean hyperbolic equations
  leh2=size(Olist);
  leh3=size(infobo7);
  if (leh3<leh2){for (j=leh3+1;j<=leh2; j++){infobo7[j]=-1;}}
  // Auxiliary list to complete the resolution function in special cases
  if (size(auxinvlist)==0) {auxinvlist[1]=0;}
  return(center,auxJ,Olist,C,t,D,H,allH,auxinvlist,infobo7);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..4)),dp;
  list flag=identifyvar();
  ideal J=x(1)^2-x(2)^2*x(3)^5, x(1)*x(3)^3+x(4)^6;
  list Lmb=1,list0(4),list0(4),list0(4),iniD(4),iniD(4),list0(4),-1;
  list L=data(J,2,4);
  list LL=determinecenter(L[1],L[2],2,4,0,0,Lmb,flag,0,-1); // Compute the first center
  LL[1];  // index of variables in the center
  LL[2];  // exponents of ideals J_4,J_3,J_2,J_1
  LL[3];  // list of orders of J_4,J_3,J_2,J_1
  LL[4];  // list of critical values
  LL[5];  // components of the resolution function t
  LL[6];  // list of D_4,D_3,D_2,D_1
  LL[7];  // list of H_4,H_3,H_2,H_1 (exceptional divisors)
  LL[8];  // list of all exceptional divisors acumulated
  LL[9];  // auxiliary invariant
  LL[10]; // intvec pointing out the last step where the function t has dropped

  ring r= 0,(x(1..4)),dp;
  list flag=identifyvar();
  ideal J=x(1)^3-x(2)^2*x(3)^5, x(1)*x(3)^3+x(4)^5;
  list Lmb=2,list0(4),list0(4),list0(4),iniD(4),iniD(4),list0(4),-1;
  list L2=data(J,2,4);
  list L3=determinecenter(L2[1],L2[2],2,4,0,0,Lmb,flag,0,-1); // Example with rational exponents in E-Coeff
  L3[1]; // index of variables in the center
  L3[2]; // exponents of ideals J_4,J_3,J_2,J_1
  L3[3]; // list of orders of J_4,J_3,J_2,J_1
  L3[4]; // list of critical values
  L3[5]; // components of the resolution function
}
////////////////////////////////////////////////////////
// idchart= identity number of the current chart
// infochart=chart[idchart] information related to the chart to blow up
// infochart= int parent,int Y,intvec a,list expJ,list Coef, list flag,                // NEEDED FOR THE RESOLUTION
//            intvec Hhist, list blwhist, module path, list hipercoef, list hiperexp   // NEEDED FOR THE OUTPUT

// NOTE: IT IS NOT NECESSARY TAKE INTO ACCOUNT "y" VARIABLES BECAUSE THE CENTER IS ALREADY GIVEN

proc Blowupcenter(list center,int idchart,int nchart,list infochart,number c,int n,int currentstep)
"USAGE: Blowupcenter(center,id,nchart,infochart,c,n,cstep);
        center, infochart lists, id, nchart, n, cstep integers, c number
COMPUTE: The blowing up at the chart IDCHART along the given center
RETURN:  new affine charts and related information, see example
EXAMPLE: example Blowupcenter; shows an example
"
{int num,i,j,k,l,parent,Y,lon,m,m2;
 intvec a,Hhist,auxHhist;
 number auxsum, auxsum2;
 list sons,aux,expJ,blexpJ,blD;
 list auxstep,Coef;
 list auxchart,auxchart1,info,flaglist;
 list auxblwhist,blwhist,hipercoef,hiperexp;
module auxpath,auxp2;

parent=idchart;
num=size(center);

// Transform to intvec the list of variables defining the center
  a=center[1];
  for (i=2;i<=num;i++){a=a,center[i];}

expJ=infochart[4];
Coef=infochart[5];
flaglist=infochart[6];
Hhist=infochart[7];
blwhist=infochart[8];
auxpath=infochart[9];
hipercoef=infochart[10];
hiperexp=infochart[11];

l=size(expJ);

// input for the loop
blexpJ=expJ;

// making the blowing up in the i-th chart
for (i=1;i<=num;i++)
{
// we assign the current number of charts +1 to the i-th chart
idchart=nchart+1;
nchart=nchart+1;
aux=idchart;
sons=sons+aux;

auxstep[i]=currentstep+1;

Y=center[i];

// The blowing up

 for (j=1;j<=l;j++){lon=size(Coef[j]);
                   if (lon==1){for (m=1;m<=n;m++){for (m2=1;m2<=num;m2++){
                                                    if (m==center[m2]){auxsum=auxsum+ expJ[j][1][m];}}}
                               blexpJ[j][1][Y]=auxsum-c;
                               auxsum=0;}                   // monomial
                   else {for (m=1;m<=n;m++){for (m2=1;m2<=num;m2++){
                                              if (m==center[m2]){auxsum=auxsum+expJ[j][1][m];
                                                                auxsum2=auxsum2+expJ[j][2][m];}}}
                         blexpJ[j][1][Y]=auxsum-c;
                         blexpJ[j][2][Y]=auxsum2-c;
                         auxsum=0; auxsum2=0;}               // binomial
                  }


auxHhist=Hhist,Y;                        // history of the exceptional divisors in this chart
auxblwhist=tradblwup(blwhist,n,Y,a,num); // history of the blow ups in this chart

auxp2=auxpath,[parent,i];

auxchart1=parent,Y,a,blexpJ,Coef,flaglist,auxHhist,auxblwhist,auxp2,hipercoef,hiperexp;

// Coef, flaglist are not modified after the blowing-up, the hyperbolic information is the same as in the parent chart

auxchart[i]=auxchart1;

// Inicializating the exponents of J for the next chart

blexpJ=expJ;
}
// end of the loop

// we add its sons to the current chart
infochart=infochart+sons;
info[1]=infochart;

return(info,auxchart,nchart,auxstep,num);
}
example
{"EXAMPLE:"; echo = 2;
ring r = 0,(x(1),y(2),x(3),y(4),x(5..7)),dp;
list flag=identifyvar();
ideal J=x(1)^3-x(3)^2*y(4)^2,x(1)*x(7)*y(2)-x(6)^3*x(5)*y(4)^3,x(5)^3-x(5)^3*y(2)^2;
list Lmb=2,list0(7),list0(7),list0(7),iniD(7),iniD(7),list0(7),-1;
list L=data(J,3,7);
list L2=determinecenter(L[1],L[2],2,7,0,0,Lmb,flag,0,-1); // Computing the center
module auxpath=[0,-1];
list infochart=0,0,0,L[2],L[1],flag,0,list0(7),auxpath,list(),list();
list L3=Blowupcenter(L2[1],1,1,infochart,2,7,0);
L3[1]; // current chart (parent,Y,center,expJ,Coef,flag,Hhist,blwhist,path,hipercoef,hiperexp) with sons: [12],...,[16]
L3[2][1]; // information of its first son, write L3[2][2],...,L3[2][5] to see the other sons
L3[3];    // current number of charts
L3[4];    // step/level associated to each son
L3[5];    // number of variables in the center
}
//////////////////////////////////////////////////////////////

proc tradblwup(list blwhist,int n,int Y,intvec a,int num)
"Internal procedure - no help and no example available
"
{
int i,j,blwnew;
intvec aux,aux2;

for (j=1;j<=n;j++){
                   for (i=1;i<=num;i++){
                                       if (j==a[i] and a[i]!=Y){blwnew=Y; break;}
                                       else {blwnew=0;}
                                       }
                   aux=blwhist[j];
                   aux2=aux,blwnew;
                   blwhist[j]=aux2;
                  }
return(blwhist);
}
//////////////////////////////////////////////////////////////
// It is called only when Eord(J)=0, and J!=1 it is checked inside
// SO IT IS CALLED AFTER: maxEord(Coef,expJ,sJ,n,flaglist); --> gives (max E-order,sums)

proc Nonhyp(list Coef,list expJ,int sJ,int n,list flaglist,list sums)
"USAGE: Nonhyp(Coef,expJ,sJ,n,flaglist,sums);
        Coef, expJ, flaglist, sums lists, sJ, n integers
COMPUTE: The "ideal" generated by the non hyperbolic generators of J
RETURN: lists with the following information
        newcoef,newJ: coefficients and exponents of the non hyperbolic generators
        totalhyp,totalgen: coefficients and exponents of the hyperbolic generators
        flaglist: new list saying status of variables
NOTE: the basering r is supposed to be a polynomial ring K[x,y],
      in fact, we work in a localization of K[x,y], of type K[x,y]_y with y invertible variables.
EXAMPLE: example Nonhyp; shows an example
"
{
int i,j,k,h,lon,lon2,cont;
number eordcontrol;
list genhyp,listgen,listid,posnumJ,newJ,newcoef,hypcoef,hyp,aux1,aux2,aux3,aux,midlist;
list totalhyp,totalgen;

eordcontrol=0;

while (eordcontrol==0 and sJ!=0)
{

// Give a positional number/flag to each generator of expJ

for (i=1;i<=sJ; i++){listgen=expJ[i]; listid=i; posnumJ[i]=listgen+listid; }

// Select the non hyperbolic and hyperbolic generators

for (j=1;j<=sJ; j++){lon=size(Coef[j]);
                     if (lon==1){

// IS NOT NECESSARY TO CHECK IF THERE EXIST A MONOMIAL WITH ONLY UNITS, ALREADY DONE!!

                                     aux1=aux1+posnumJ[j];
                                     aux3=list();
                                     aux3[1]=expJ[j];
                                     newJ=newJ+aux3;
                                     aux3[1]=Coef[j];
                                     newcoef=newcoef+aux3;
                                }

                     else{   // CHECKING BINOMIALS, ONE TERM WITH E-ORDER ZERO GIVES HYPERBOLIC EQ

                          if (sums[j][1]==0 or sums[j][2]==0){aux2=aux2+posnumJ[j];
                                                              aux3=list();
                                                              aux3[1]=expJ[j];
                                                              genhyp=genhyp+aux3;
                                                              aux3[1]=Coef[j];
                                                              hypcoef=hypcoef+aux3;
                                                              if (sums[j][1]==0){aux3[1]=expJ[j][2]; hyp=hyp+aux3;}
                                                              if (sums[j][2]==0){aux3[1]=expJ[j][1]; hyp=hyp+aux3;}
                                                             }
                          else {aux1=aux1+posnumJ[j];
                                aux3=list();
                                aux3[1]=expJ[j];
                                newJ=newJ+aux3;
                                aux3[1]=Coef[j];
                                newcoef=newcoef+aux3;}

                          }
                    }

// NOTE: aux1 and aux2 are no needed right now!

// Identify new y variables, that is, x variables in the monomials contained in hyp

h=size(hyp);

for (k=1;k<=h; k++){ for(i=1;i<=n; i++){ if (hyp[k][i]!=0 and flaglist[i]==0) {flaglist[i]=1;}}}

// To replace x by y IT IS NECESSARY TO CHANGE THE BASERING!!! We change only the list flaglist

// CHECK IF THE IDEAL IS ALREADY GENERATED BY MONOMIALS, in this case
// WE HAVE FINISHED THE E-RESOLUTION PART, J GENERATED BY MONOMIALS AND HYPERBOLIC EQS

cont=0;
lon2=size(newJ);
for (j=1;j<=lon2; j++){if (size(newJ[j])==1){cont=cont+1;}}

if (cont==lon2){newcoef=list();
                newJ=list();
                totalgen=totalgen+genhyp;
                totalhyp=totalhyp+hypcoef;
                break;}

// CHECK IF THERE ARE MORE HYPERBOLIC EQUATIONS AFTER UPDATE THE FLAG LIST
// CHECK THE MAXIMAL E-ORDER AGAIN

if (lon2==0){   // we are in the previous case, newJ=empty list, save values and exit

             totalgen=totalgen+genhyp;
             totalhyp=totalhyp+hypcoef;
             break;
             }

midlist=maxEord(newcoef,newJ,lon2,n,flaglist);

eordcontrol=midlist[1];

  if (eordcontrol==0){                  // new input for the loop
                      Coef=newcoef;
                      expJ=newJ;
                      sJ=lon2;
                      sums=midlist[2]; // flaglist is already updated

                      totalgen=totalgen+genhyp;
                      totalhyp=totalhyp+hypcoef;

                      hypcoef=list();
                      genhyp=list();

                      newJ=list();
                      newcoef=list();
                    }
else{  // If the process is already finished we save the values and exit

     totalgen=totalgen+genhyp;
     totalhyp=totalhyp+hypcoef;
    }

} // closing while

return(newcoef,newJ,totalhyp,totalgen,flaglist);
}
example
{"EXAMPLE:"; echo = 2;
ring r = 0,(x(1),y(2),x(3),y(4),x(5..7)),dp;
list flag=identifyvar();  // List giving flag=1 to invertible variables: y(2),y(4)
ideal J=x(1)^3-x(3)^2*y(4)^2,x(1)*x(7)*y(2)-x(6)^3*x(5)*y(4)^3,1-x(5)^2*y(2)^2;
list L=data(J,3,7);
list L2=maxEord(L[1],L[2],3,7,flag);
L2[1];     // Maximum E-order
list New=Nonhyp(L[1],L[2],3,7,flag,L2[2]);
New[1];    // Coefficients of the non hyperbolic part
New[2];    // Exponents of the non hyperbolic part
New[3];    // Coefficients of the hyperbolic part
New[4];    // New hyperbolic equations
New[5];    // New list giving flag=1 to invertible variables: y(2),y(4),y(5)

ring r = 0,(x(1..4)),dp;
list flag=identifyvar();
ideal J=1-x(1)^5*x(2)^2*x(3)^5, x(1)^2*x(3)^3+x(1)^4*x(4)^6;
list L=data(J,2,4);
list L2=maxEord(L[1],L[2],2,4,flag);
L2[1];     // Maximum E-order
list New=Nonhyp(L[1],L[2],2,4,flag,L2[2]);
New;

}
//////////////////////////////////////////////////////////////

proc calculateI(list Coef,list expJ,number c,int n,int Y,intvec a,number oldordI,list oldD)
"USAGE: calculateI(Coef,expJ,c,n,Y,a,b,D);
        Coef, expJ, D lists, c, b numbers, n,Y integers, a intvec
RETURN: ideal I, non monomial part of J
EXAMPLE: example calculateI; shows an example
"
{
 int i,cont1,b,j;
 number EordI,aux;
 list D,L,expI;
 list auxdiv,Dstar,aux1,rs;

// WE NEED THE MONOMIAL PART, BUT ONLY IN DIMENSION n

     auxdiv=list0(n);
     auxdiv[Y]=oldordI-c;
     Dstar[1]=oldD[1];

     b=size(a);
     for (i=1;i<=n;i++) {for (j=1;j<=b;j++) {if (a[j]==i) {aux=aux+oldD[1][i];}}}
     Dstar[1][Y]=aux;
     aux=0;

     D[1]=sumlist(Dstar[1],auxdiv);

   cont1=0;
   for (i=1;i<=n;i++) {if (D[1][i]==0) {cont1=cont1+1;}}   // if it fails write listO(n)[i]

   if (cont1==n) {expI=expJ;}
   else {
         for (i=1;i<=size(expJ);i++)
         {rs[i]=size(Coef[i]);
          if (rs[i]==2){ aux1=list();
                         aux1[1]=reslist(expJ[i][1],D[1]);
                         aux1[2]=reslist(expJ[i][2],D[1]);
                         expI[i]=aux1;}                            // binomial
          else {aux1=list();
                aux1[1]=reslist(expJ[i][1],D[1]);
                expI[i]=aux1;}}                                    // monomial
        }

return(expI);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..3)),dp;
  list flag=identifyvar();
  ideal J=x(1)^4*x(2)^2, x(3)^3;
  list Lmb=1,list0(3),list0(3),list0(3),iniD(3),iniD(3),list0(3),-1;
  list L=data(J,2,3);
  list LL=determinecenter(L[1],L[2],3,3,0,0,Lmb,flag,0,-1); // Calculate the center
  module auxpath=[0,-1];
  list infochart=0,0,0,L[2],L[1],flag,0,list0(3),auxpath,list(),list();
  list L3=Blowupcenter(LL[1],1,1,infochart,3,3,0); // blowing-up and looking to the x(3) chart
  calculateI(L3[2][1][5],L3[2][1][4],3,3,3,L3[2][1][3],3,iniD(3)); //  (I_3)
  // looking to the x(1) chart
  calculateI(L3[2][2][5],L3[2][2][4],3,3,1,L3[2][2][3],3,iniD(3)); //  (I_3)
}
//////////////////////////////////////////////////////////////////////////////////////
//                                                                                  //
//  E-RESOLUTION: Eresol(J) subroutine computing the E-resolution of J, char 0      //
//                                                                                  //
//////////////////////////////////////////////////////////////////////////////////////

proc Eresol(ideal J)
"USAGE: Eresol(J); J ideal
RETURN: The E-resolution of singularities of J in terms of the affine charts, see example
EXAMPLE: example Eresol; shows an example
"
{int i,n,k,idchart,nchart,parent,Y,oldid,tnum,s,cont,control,control2,control3,cont2,val,rs2,l,cont3,tip;
 intvec a,Hhist;
 number c,EordJ,EordI,oldordI;
 list L,LL,oldD,t,auxL,finalchart,chart,auxchart,newL,auxp,auxfchart,L2;
 list Coef,expJ,expI,sons,oldOlist,oldC,oldt,oldH,allH,auxordJ,auxordI,auxmb,mobile,invariant;
 list step,nsons,auxinv,extraL,totalinv,auxsum;
 string empstring;
 module auxpath;

// ADDED LATER

 list flag,newflag,blwhist,hipercoef,hiperexp,hipercoefson,hiperexpson;
 intvec infobo7;

export finalchart;
// export nsons;
// export tnum;
// export nchart;
// export step;
export invariant;
 export auxinv;
export mobile;

 n=nvars(basering);
 flag=identifyvar();

 k=size(J);
// Checking input data
 if (inidata(J,k)==0){return("This library only works for binomial ideals.");}

idchart=1;
nchart=1;
parent=0;
step=0;
control=0;
control2=0;
control3=0;

// Translate the input ideal to a list
 auxL=data(J,k,n);                       // data gives (Coef,Exp)

// THEREAFTER WE WORK ALL THE TIME WITH LISTS

L=maxEord(auxL[1],auxL[2],k,n,flag);   // gives (max E-ord, sums)
EordJ=L[1];

// before the first blow up I=J
EordI=EordJ;

//  main loop AT EACH CHART WE MUST INICIALIZATE ALL THE VALUES AND
// CONSTRUCT THE FIRST CHART chart[1] BEFORE THE LOOP

// at the first step, before the blow up, there are no exceptional divisors, Y=0
Y=0;
expJ=auxL[2];
Coef=auxL[1];
Hhist=0;
blwhist=list0(n);
auxpath=[0,-1];
hipercoef=list();   // this is for the first chart
hiperexp=list();
auxp=parent,Y,a,expJ,Coef,flag,Hhist,blwhist,auxpath,hipercoef,hiperexp;
chart[1]=auxp;      // information of the first chart

tip=1;
oldOlist=list0(n);
oldC=list0(n);
oldC[1]=EordJ;    // non necessary here
c=EordJ;          // the value c is given by the previous step
oldt=list0(n);
oldD=iniD(n);
oldH=iniD(n);
allH=list0(n);

for (i=1;i<=n;i++){infobo7[i]=-1;}

auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7;
mobile[1]=auxmb;           // mobile corresponding to the first chart
auxinv[1]=list(0);

// NOTE: oldC[1] is the value c to classify the chart in one of the next cases

// HERE BEGIN THE LOOP

while (idchart<=nchart)   // WE PROCEED WHILE THERE EXIST UNSOLVED CHARTS
{
if (idchart!=1)           // WE ARE NOT IN THE FIRST CHART, INICIALIZATE ALL THE VALUES
{

parent=chart[idchart][1];
Y=chart[idchart][2];
a=chart[idchart][3];
expJ=chart[idchart][4];
Coef=chart[idchart][5];
flag=chart[idchart][6];
Hhist=chart[idchart][7]; // it is not necessary for the computations
blwhist=chart[idchart][8];
auxpath=chart[idchart][9];
hipercoef=chart[idchart][10];
hiperexp=chart[idchart][11];

k=size(Coef);   // IT IS NECESSARY TO COMPUTE IT BECAUSE IT DECREASES IF THERE ARE HYPERBOLIC EQS

auxordJ=maxEord(Coef,expJ,k,n,flag);
EordJ=auxordJ[1];

if (control==0){c=mobile[parent+1][3][1];}  // we keep c from the last step
else {c=EordJ; control=0; }                  // we reset the value of c

 if (control2==1){c=EordJ; control2=0; control3=1;}    // we reset the value of c

// NOTE: oldC[1] is the value c to classify the chart in one of the next cases

}

// The E-order must be computed here

oldid=idchart;

if (EordJ<0) {print("ERROR in J in chart"); print(idchart); break;}


//-------------------------------------------------------------
// CASE J=1, if we reset c, can happen Eord=c=0

// or if there are hyperbolic equations at the beginning!!!

// if (EordJ==0){auxfchart[1]=chart[idchart];       // WE HAVE FINISHED
//              finalchart=finalchart+auxfchart;
//              empstring="#";                   print("reset c and Eord=c=0"); print(idchart);
//              invariant[idchart]=empstring;
//              auxinv[idchart]=list(0);
//              nsons[idchart]=0;
//              idchart=idchart+1;}


//----------------------------------------------------------------------
if (EordJ>=c and EordJ!=0)      // subroutine: E-RESOLUTION OF PAIRS
{
  if (parent>0)
  { LL=determinecenter(Coef,expJ,c,n,Y,a,mobile[parent+1],flag,control3,chart[parent][7]); }
  else { LL=determinecenter(Coef,expJ,c,n,Y,a,mobile[parent+1],flag,control3,Hhist); }

// determinecenter gives (center,auxJ,Olist,C,t,D,H,allH,auxinvlist,infobo7)

// save current values, before the blow up
oldOlist=LL[3];
tip=computemcm(oldOlist);
oldC=LL[4];
oldt=LL[5];
oldD=LL[6];
oldH=LL[7];
allH=LL[8];
auxinv[idchart]=LL[9];
infobo7=LL[10];

auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7;
mobile[idchart+1]=auxmb;
invariant[idchart]=oldt;

newL=Blowupcenter(LL[1],idchart,nchart,chart[idchart],c,n,step[idchart]);

// Blowupcenter gives (info,auxchart,nchart,auxstep,num)

// IMPORTANT: ADD THE NEW CHARTS AFTER EACH BLOW UP, IN ORDER TO KEEP THEM CORRECTLY

step=step+newL[4];
nsons[idchart]=newL[5];

chart=chart+newL[2];
finalchart=finalchart+newL[1];

// new input for the loop

idchart=idchart+1;
nchart=newL[3];

control3=0;

} // END OF CASE EordJ>=c
//---------------------------------------------------------------------

else{

// compute EordI=max E-order(I)

expI=calculateI(Coef,expJ,c,n,Y,a,mobile[parent+1][2][1],mobile[parent+1][5]);
k=size(expJ);                     // probably non necessary
auxordI=maxEord(Coef,expI,k,n,flag);
EordI=auxordI[1];
auxsum=auxordI[2];

// CASE EordI>0  DROP c AND CONTINUE

if (EordI>0){idchart=idchart;  // keep the chart and back to the main loop while, dropping the value of c
             control=1;}
else{                          // EordI=0, so check if I=1 or not

cont2=0;                       // If cont2=val then all the entries of expI are zero!!
val=0;

for (i=1;i<=k;i++) {rs2=size(Coef[i]);
                    if (rs2==1){if (auxsum[i][1]==0){cont2=val; break;} // THERE EXIST A MONOMIAL WITH ONLY UNITS

                                val=val+n;                                               // monomials
                                for (l=1;l<=n; l++) {if (expI[i][1][l]==0) {cont2=cont2+1;}}
                               }
                    else{val=val+(2*n);                                                     // binomials
                         for (l=1;l<=n; l++) {if (expI[i][1][l]==0) {cont2=cont2+1;}
                                              if (expI[i][2][l]==0) {cont2=cont2+1;}}
                        }
                    }


// CASE EordI==0 AND I=1 THIS CHART IS DONE, FINISH

// NOTE: THIS CASE IS NOT MONOMIAL BECAUSE E-Sing(J,c) is empty

if (cont2==val){auxfchart[1]=chart[idchart];
                finalchart=finalchart+auxfchart;
                empstring="#";
                invariant[idchart]=empstring;
                auxinv[idchart]=list(0);
                nsons[idchart]=0;

// information for the mobile
                tip=1;
                oldOlist=list(0);
                oldC=list(0);
                oldt=list(0);
                oldD=list(0);
                oldH=list(0);
                allH=list(0);  // the value of the parent + the new one
                infobo7=-1;

                auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7;
                mobile[idchart+1]=auxmb;

                idchart=idchart+1;}

else{         // CASE EordI==0 AND I!=1 --> HYPERBOLIC EQUATIONS

// COMPUTE THE IDEAL OF NON HYPERBOLIC ELEMENTS

extraL=Nonhyp(Coef,expI,k,n,flag,auxordI[2]);    // gives (newcoef,newI,hypcoef,genhyp,flaglist)

// CHECK IF ALL THE VARIABLES ARE ALREADY INVERTIBLE

newflag=extraL[5];
chart[idchart][6]=extraL[5];          // update the status of variables

cont3=0;
for (i=1;i<=n;i++){if (newflag[i]==1){cont3=cont3+1;}}

if (cont3==n){    // ALL THE VARIABLES ARE INVERTIBLE
              auxfchart[1]=chart[idchart];
              finalchart=finalchart+auxfchart;
              empstring="@";
              invariant[idchart]=empstring;
              auxinv[idchart]=list(0);
              nsons[idchart]=0;

// information for the mobile
                tip=1;
                oldOlist=list(0);
                oldC=list(0);
                oldt=list(0);
                oldD=list(0);
                oldH=list(0);
                allH=list(0);
                infobo7=-1;

                auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7;
                mobile[idchart+1]=auxmb;

              idchart=idchart+1;}
else{     // OTHERWISE, CONTINUE CHEKING IF newI=0 or not

Coef=extraL[1];
expI=extraL[2];

hipercoefson=extraL[3];  // Information about hyperbolic generators
hiperexpson=extraL[4];

k=size(expI);

if (k==0){auxfchart[1]=chart[idchart];       // WE HAVE FINISHED
          finalchart=finalchart+auxfchart;
          empstring="#";                  //  no more non-hyperbolic generators in this chart
          invariant[idchart]=empstring;
          auxinv[idchart]=list(0);
          nsons[idchart]=0;

// information for the mobile
                tip=1;
                oldOlist=list(0);
                oldC=list(0);
                oldt=list(0);
                oldD=list(0);
                oldH=list(0);
                allH=list(0);
                infobo7=-1;

                auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7;
                mobile[idchart+1]=auxmb;

          idchart=idchart+1;}

else{           // CONTINUE WITH THE IDEAL OF NON HYPERBOLIC EQS

     chart[idchart][4]=expI;   // new input ideal and coefficients
     chart[idchart][5]=Coef;
     chart[idchart][10]=hipercoef+hipercoefson;
     chart[idchart][11]=hiperexp+hiperexpson;

     idchart=idchart;
     control2=1;   // it is necessary to reset the value of c
     control3=1;   // and the previous exceptional divisors
    }

// PROBABLY IT IS NEC MORE INFORMATION !!!

} // closing else otherwise

    } // closing else case I!=1

} // closing else for EordI=0

if (EordI<0) {print("ERROR in chart"); print(idchart); break;}

//----------------------- guardar de momento--------
// if (EordI==0) {auxfchart[1]=chart[idchart];
//         finalchart=finalchart+auxfchart;
//         L2=Gamma(expJ,c,n); // HAY QUE APLICARLO AL M NO AL J
//         invariant[idchart]=L2[2];
//         auxinv[idchart]=list(0);
//         nsons[idchart]=0;
//         idchart=idchart+1;}
//------------------------------------------------


} // END ELSE
//---------------------------------------------------

} // END LOOP WHILE

tnum=step[nchart];

totalinv=resfunction(invariant,auxinv,nchart,n);

return(chart,finalchart,invariant,nchart,step,nsons,auxinv,mobile,totalinv);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..2)),dp;
  ideal J=x(1)^2-x(2)^3;
  list L=Eresol(J);
  L[1][1]; // information of the first chart, L[1] list of charts
  L[2]; // list of charts with information about sons
  L[3]; // invariant, "#" means solved chart
  L[4]; // number of charts, 7 in this example
  L[5]; // height corresponding to each chart
  L[6]; // number of sons
  L[7]; // auxiliary invariant
  L[8]; // H exceptional divisors and more information
  L[9]; // complete resolution function

"Second example, write L[i] to see the i-th component of the list";
  ring r = 0,(x(1..3)),dp;
  ideal J=x(1)^2*x(2),x(3)^3; // SOLVED!
  list L=Eresol(J);
  L[4];  // 16 charts
  L[9]; // complete resolution function

"Third example, write L[i] to see the i-th component of the list";
  ring r = 0,(x(1..2)),dp;
  ideal J=x(1)^3-x(1)*x(2)^3;
  list L=Eresol(J);
  L[4];  // 8 charts, rational exponents
  L[9]; // complete resolution function
}

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

proc resfunction(list invariant, list auxinv, int nchart,int n)
"USAGE: resfunction(invariant,auxinv,nchart,n);
        invariant, auxinv lists, nchart, n integers
COMPUTE: Patch the resolution function
RETURN: The complete resolution function
EXAMPLE: example resfunction; shows an example
"
{
int i,j,l,k;
list patchfun,aux;

for (i=1;i<=nchart;i++){patchfun[i]=invariant[i];}

for (i=1;i<=nchart;i++){if (auxinv[i][1]!=0 and size(auxinv[i])==3){l=size(invariant[i]);
                                                                    for (j=1;j<=l;j++){
                                                                         if (invariant[i][j]==0){aux=auxinv[i];
                                                                                                 patchfun[i][j]=aux;
                                                                    if (l<n){for (k=j+1;k<=n;k++){patchfun[i][k]="*";}}}}

                                                                    }
                        else{
                        if (auxinv[i][1]==1 and size(auxinv[i])==1){l=size(invariant[i]);
                                                                    if (l<n){for (k=l+1;k<=n;k++){patchfun[i][k]="*";}}
                                                                    }
                            }
                        }

return(patchfun);
}
example
{"EXAMPLE:"; echo = 2;
 ring r = 0,(x(1..2)),dp;
 ideal J=x(1)^2-x(2)^3;
 list L=Eresol(J);
 L[3]; // incomplete resolution function
 resfunction(L[3],L[7],7,2); // complete resolution function
}
//////////////////////////////////////////////////////////////////////////////////////
//                                                                                  //
//                              MAIN PROCEDURE                                      //
//                                                                                  //
//////////////////////////////////////////////////////////////////////////////////////

proc BINresol(ideal J)
"USAGE: BINresol(J); J ideal
RETURN: E-resolution of singularities of a binomial ideal J in terms of the affine charts, see example
EXAMPLE: example BINresol; shows an example
"
{

int p,n;

p=char(basering);
n=nvars(basering);  // YA SE CALCULA EN Eresol, MEJORAR?

if (p>0){list Lring=ringlist(basering);
         Lring[1]=0;
         def r=basering;
         def Rnew=ring(Lring);
         setring Rnew;
         ideal chy=maxideal(1);
         map fRnew=r,chy;
         ideal J=fRnew(J);

// E-RESOLUTION, Computations in char 0

         list L=Eresol(J);

// STEP 2: WRITE THE LOCALLY MONOMIAL IDEAL AS A MONOMIAL IDEAL

// not implemented yet, CHAR p !!!!

// STEP 3: DO THE E-RESOLUTION AGAIN (char 0 again)


// generating output

         int q=lcmofall(L[4],L[8]);    // lcm of the denominators

         list B=genoutput(L[1],L[8],L[4],L[6],n,q); // generate output needed for visualization


//         setring r;                   // Back to the basering
//         ideal chy=maxideal(1);
//         map fr=Rnew,chy;
//         list L=fr(L);
//         list B=fr(B);

         }

else{

// E-RESOLUTION

 list L=Eresol(J);

// STEP 2: WRITE THE LOCALLY MONOMIAL IDEAL AS A MONOMIAL IDEAL

// not implemented yet

// STEP 3: DO THE E-RESOLUTION AGAIN


// generating output

      int q=lcmofall(L[4],L[8]);

      list B=genoutput(L[1],L[8],L[4],L[6],n,q);

     }

return(B);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..2)),dp;
  ideal J=x(1)^2-x(2)^3;
  list B=BINresol(J);
  B[1]; // list of final charts
  B[2]; // list of all charts

  ring r2 = 2,(x(1..3)),dp;
  ideal J=x(1)^2-x(2)^2*x(3)^2;
  list B2=BINresol(J);
  B2[2]; // list of all charts
}
///////////////////////////////////////////////////////

proc Maxord(list L,int n)
"USAGE: Maxord(L,n); L list, n integer
COMPUTE: Find the maximal entry of a list, input is a list defining a monomial
RETURN: maximum entry of a list and its position
EXAMPLE: example Maxord; shows an example
"
{int i,can;
 number canmax;
 list aux;

canmax=1;
can=1;
for (i=1;i<=n;i++)
{ if (L[i]>=canmax and i>=can)
 {canmax=L[i]; can=i;}}

return(canmax,can);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..3)),dp;
  ideal J=x(1)^2*x(2)*x(3)^5;
  list L=data(J,1,3);
  L[2]; // list of exponents
  Maxord(L[2][1][1],3);
}
///////////////////////////////////////////////////////

proc Gamma(list expM,number c,int n)
"USAGE: Gamma(L,c,n); L list, c number, n integer
COMPUTE: The Gamma function, resolution function corresponding to the monomial case
RETURN: lists of maximum exponents in L, value of Gamma function, center of blow up
EXAMPLE: example Gamma; shows an example
"
{int i,j,k,l,cont,can;
 intvec upla;
 number canmax;
 list expM2,gamma,L,aux,maxlist,center,aux2;

 i=1;
 cont=0;
 expM2=expM;

 while (cont==0 and i<=n)
{

  L=Maxord(expM2,n);
  aux=L[1];
  maxlist=maxlist + aux;
  can=L[2];

if (i==1) {upla=can; center=can;}
else {upla=upla,can; aux2=can; center=center+aux2;}

  canmax=sum(maxlist);
  if (canmax>=c)
  {gamma[1]=-i; gamma[2]=canmax/c; gamma[3]=upla; cont=1;}
  else {expM2[can]=0;}
  i=i+1;
}
 return(maxlist,gamma,center);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..5)),dp;
  ideal J=x(1)^2*x(2)*x(3)^5*x(4)^2*x(5)^3;
  list L=data(J,1,5);
  list G=Gamma(L[2][1][1],9,5);   // critical value c=9
  G[1]; // maximum exponents in the ideal
  G[2]; // maximal value of Gamma function
  G[3]; // center given by Gamma
}
///////////////////////////////////////////////////////

proc convertdata(list C,list L, int n, list flaglist)
"USAGE: convertdata(C,L,n,flaglist);
        C, L, flaglist lists, n integer
COMPUTE: Compute the ideal corresponding to the given lists C,L
RETURN: an ideal whose coefficients are given by C, exponents given by L
EXAMPLE: example convertdata; shows an example
"
{int i,j,k,a,b,lon;
 poly aux,aux1,aux2,aux3,f;
 ideal J;

aux=poly(0);
aux1=poly(1);
aux2=poly(0);
aux3=poly(1);


k=size(L);
for (i=1;i<=k;i++){lon=size(C[i]);
if (lon==1){                                                // variables in the monomial
            for (j=1;j<=n;j++){a=int(poly(L[i][1][j]));
                               if (a!=0){
                                         if (flaglist[j]==0){aux=poly(x(j)^a);
                                                             aux1=aux1*aux;}
                                         else {aux=poly(y(j)^a);
                                               aux1=aux1*aux;}
                                        }
                              }
            if (C[i][1]!=0){aux1=C[i][1]*aux1;}            // we add the coefficient
            else {aux1=0;}

            J[i]=aux1;
            aux1=poly(1);
           }

else{                                                    // variables in the binomial

     for (j=1;j<=n;j++){a=int(poly(L[i][1][j])); b=int(poly(L[i][2][j]));

                        if (a!=0){
                                  if (flaglist[j]==0){aux=poly(x(j)^a);
                                                      aux1=aux1*aux;}
                                  else {aux=poly(y(j)^a);
                                        aux1=aux1*aux;}
                                 }

                        if (b!=0){
                                  if (flaglist[j]==0){aux2=poly(x(j)^b);
                                                      aux3=aux3*aux2;}
                                  else {aux2=poly(y(j)^b);
                                        aux3=aux3*aux2;}
                                 }
                          }
                                                       // we add the coefficients
   if (C[i][1]!=0){aux1=C[i][1]*aux1;}
    else {aux1=0;}
   if (C[i][2]!=0){aux3=C[i][2]*aux3;}
    else {aux3=0;}

   f=aux1+aux3;
   J[i]=f;
   aux1=poly(1);
   aux3=poly(1);

     }
}
return(J);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..4),y(5)),dp;
  list M=identifyvar();
  ideal J=x(1)^2*y(5)^2-x(2)^2*x(3)^2,6*x(4)^2;
  list L=data(J,2,5);
  L[1]; // Coefficients
  L[2]; // Exponents
  ideal J2=convertdata(L[1],L[2],5,M);
  J2;
}

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

proc lcmofall(int nchart,list mobile)
"USAGE: lcmofall(nchart,mobile);
        nchart integer, mobile list of lists
COMPUTE: Compute the lcm of the denominators of the E-orders of all the charts
RETURN: an integer given the lcm
NOTE: CALL BEFORE salida
EXAMPLE: example lcmofall; shows an example
"

{
int i,m,tip,mcmall;
intvec numall;

for (i=2;i<=nchart+1;i++){
                          tip=mobile[i][1];
                          if (tip!=1){numall=numall,tip;}
                         }
m=size(numall);

if (m==1){mcmall=1;}
else{
     if (numall[1]==0){numall=numall[2..m];}
     mcmall=lcm(numall);}

return(mcmall);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..2)),dp;
  ideal J=x(1)^3-x(1)*x(2)^3;
  list L=Eresol(J);
  L[4];  // 8 charts, rational exponents
  L[8][2][2]; // E-orders at the first chart
  lcmofall(8,L[8]);
}
/////////////////////////////////////////////////////////////////////////////

proc salida(int idchart,list chart,list mobile,int numson,intvec previousa,int n,int q)
"USAGE: salida(idchart,chart,mobile,numson,previousa,n,q);
        idchart, numson, n, q integers, chart, mobile, lists, previousa intvec
COMPUTE: CONVERT THE OUTPUT OF A CHART IN A RING, WHERE DEFINE A BASIC OBJECT (BO)
RETURN: the ring corresponding to the chart
EXAMPLE: example salida; shows an example
"
{
int l,i,m,aux,parent,m4,j;
intvec Hhist,EOhist,aux7,aux9;
list expJ,Coef,BO,blwhist,Eolist,hipercoef,hiperexp;
list flag;

// chart gives: parent,Y,a,expJ,Coef,flag,Hhist,blwhist,path,hipercoef,hiperexp
// mobile gives: tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7;  NOTE: Eolist=mobile[2];

// we need to define the suitable ring at this chart

list Lring=ringlist(basering);
def RR2=basering;

flag=chart[6];
string newl;

for (l=1;l<=n; l++){if (flag[l]==1){newl=string(l);
                                    Lring[2][l]="y("+newl+")";} }


         def RRnew=ring(Lring);
         setring RRnew;
         ideal chy=maxideal(1);
         map fRnew=RR2,chy;

         list chart=fRnew(chart);

         list mobile2=fRnew(mobile);


flag=chart[6];

// we need to convert expJ and Coef to an ideal

expJ=chart[4];
Coef=chart[5];
Hhist=chart[7];
blwhist=chart[8];

// now the ideal will be correctly defined in the ring Rnew

ideal J2=convertdata(Coef,expJ,n,flag);        // Computations in RRnew

//------------------------------------------------------------------------------
// START TO CREATE THE BO corresponding to this chart

BO=createBO(J2);

// MODIFY BO WITH THE INFORMATION OF THE CHART

// BO[1] an ideal, say W_i, defining the ambient space of the i-th chart of the blowing up
// If there are hyperbolic equations, we put them here

hipercoef=chart[10];
hiperexp=chart[11];

if (size(hipercoef)!=0){
                        ideal ambJ=convertdata(hipercoef,hiperexp,n,flag);
                        BO[1]=ambJ;
                       }

// BO[2] an ideal defining the controlled transform

BO[2]=J2;

// BO[3] intvec, tupla containing the maximal E-order of BO[2]

if (numson==0){BO[3]=1;}    // we write 1 if the chart is a final chart
else{
     Eolist=mobile2[2];     // otherwise, convert the list of E-orders in an intvec
     m=size(Eolist);
     aux=int(Eolist[1]*q);
     EOhist=aux;

     if (m>1){for (i=2;i<=m;i++){aux=int(Eolist[i]*q); EOhist=EOhist,aux;}}

     BO[3]=EOhist;
     }

// BO[4] the list of exceptional divisors given by Hhist

 BO[4]=constructH(Hhist,n,flag);

// BO[5] an ideal defining the map K[W] ----> K[Wi] given by blwhist

BO[5]=constructblwup(blwhist,n,chy,flag);

// BO[6] an intvec, BO[6][j]=1 indicates that <BO[4][j],BO[2]>=1, i.e. the
//                  strict transform does not meet the j-th exceptional divisor

m4=size(BO[4]);
ideal auxydeal;
ideal Jint;

for (j=1;j<=m4;j++){

                    auxydeal=BO[4][j]+J2;
                    Jint=std(auxydeal);

                    if (size(Jint)==1 and Jint[1]==1){BO[6][j]=1;}
                    else{BO[6][j]=0;}
                   }

// BO[7] intvec, the index of the first blown-up object in the resolution process
//                leading to this object for which the value of b was BO[3]
//                the subsequent ones are the indices for the Coeff-Objects
//                of BO[2] used when determining the center
//   index of last element of H^- in H


if (numson!=0){BO[7]=mobile2[8];}  // it is always -1 at the final charts

// BO[8] a matrix indicating that BO[4][i] meets BO[4][j] by BO[8][i,j]=1 for i < j

if (m4>0){
matrix aux8[m4][m4];

BO[8]=aux8;

ideal auxydeal2;
ideal Jint2;

for (i=1;i<=m4;i++){
                    for (j=i+1;j<=m4;j++){
                                        auxydeal2=BO[4][i]+BO[4][j];
                                        Jint2=std(auxydeal2);

                                        if (size(Jint2)==1 and Jint2[1]==1){BO[8][i,j]=0;}
                                        else{ for (l=1;l<j;l++){BO[8][l,j]=1;} }
                                        }

                    }
}
else{ matrix aux8[1][1];
      BO[8]=aux8;}


// BO[9] INTERNAL DATA, second component of Villamayor resolution function,
// only needed to use the visualization procedures

 int m3=size(BO[3]);

if (m3==1){aux9=intvec(0);}
else{ aux9[1]=0;
      for (i=2;i<=m3;i++){aux9=aux9,0;}
    }

BO[9]=aux9;

//------------------------------------------------------------------------------

// START TO CREATE THE extra information corresponding to this chart

/////////////// Short description of data in a chart ///////////////////
// All chart data is stored in an object of type ring, the following
// variables are always present in such a ring:

// BO:      already created

// cent:  ideal, describing the upcoming center determined by the algorithm

   ideal cent=tradtoideal(previousa,J2,flag);
   export cent;

// path= module (autoconverted to matrix)
// path[1][idchart]=parent[idchart] index of the parent-chart in resolution history of this chart
// path[2][idchart]=index of this chart in relation with its brother-charts

   module path=chart[9];
   export path;

// lastMap: ideal, describing the preceding blow up leading to this chart

   ideal lastMap=constructlastblwup(blwhist,n,chy,flag);
   export lastMap;

//------------------------------------------------------------------------------

// EXTRA INFORMATION NEEDED

  list invSat=ideal(0),aux9;
  export(invSat);

export BO;

return(RRnew);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..2)),dp;
  ideal J=x(1)^2-x(2)^3;
  list L=Eresol(J);
  list B=salida(5,L[1][5],L[8][6],2,L[1][3][3],2,1); // chart 5
  def RR=B[1];
  setring RR;
  BO;

  ring r = 0,(x(1..2)),dp;
  ideal J=x(1)^2-x(2)^3;
  list L=Eresol(J);
  list B=salida(7,L[1][7],L[8][8],0,L[1][5][3],2,1); // chart 7
  def RR=B[1];
  setring RR;
  BO;
  showBO(BO);

  ring r = 0,(x(1..2)),dp;
  ideal J=x(1)^3-x(1)*x(2)^3;
  list L=Eresol(J); // 8 charts, rational exponents
  list B=salida(1,L[1][1],L[8][2],2,0,2,2); // CHART 1
  def RR=B[1];
  setring RR;
  BO;

}

/////////////////////////////////////////////////////////////////////////////
// CONVERT THE OUTPUT OF Eresol IN A LIST OF RINGS, WHERE A BASIC OBJECT (BO) IS DEFINED
// IN ORDER TO INTEGRATE THIS LIBRARY INSIDE THE LIBRARY resolve.lib

proc genoutput(list chart,list mobile,int nchart,list nsons,int n,int q)
"USAGE: genoutput(chart,mobile,nchart,nsons,n,q);
        chart, mobile, nsons lists, nchart, n, q integers
RETURN: two lists, the first one gives the rings corresponding to the final charts,
        the second one is the list of all rings corresponding to the affine charts of the resolution process
EXAMPLE: example genoutput; shows an example
"
{
  int idchart,parent;
  list auxlist,solvedrings,totalringlist,previousa;

  // chart gives: parent,Y,a,expJ,Coef,flag,Hhist,blwhist,path,hipercoef,hiperexp
  // mobile gives: tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7;  NOTE: Eolist=mobile[2];

  idchart=1;

  // first loop, construct list previousa

  while (idchart<=nchart)
  {
    if (idchart==1){previousa[1]=chart[2][3];}
    else
    {
      // if there are no sons, the next center is nothing
      if (nsons[idchart]==0){previousa[idchart]=0;}
      // always fill the parent
      parent=chart[idchart][1];
      previousa[parent]=chart[idchart][3];
    }
    idchart=idchart+1;
  }

  // HERE BEGIN THE LOOP

  idchart=1;

  while (idchart<=nchart)
  {
    def auxexit=salida(idchart,chart[idchart],mobile[idchart+1],nsons[idchart],previousa[idchart],n,q);
    // we add the ring to the list of all rings
    auxlist[1]=auxexit;
    totalringlist=totalringlist+auxlist;
    // if the chart has no sons, add it to the list of final charts
    if (nsons[idchart]==0){solvedrings=solvedrings+auxlist;}
    auxlist=list();
    kill auxexit;
    idchart=idchart+1;
  } // EXIT WHILE
  return(solvedrings,totalringlist);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..2)),dp;
  ideal J=x(1)^3-x(1)*x(2)^3;
  list L=Eresol(J);         // 8 charts, rational exponents
  list B=genoutput(L[1],L[8],L[4],L[6],2,2); // generates the output
  presentTree(B);
  list iden0=collectDiv(B);
  ResTree(B,iden0[1]);        // generates the resolution tree

  // Use presentTree(B); to see the final charts
  // To see the tree type in another shell
  //    dot -Tjpg ResTree.dot -o ResTree.jpg
  //   /usr/bin/X11/xv ResTree.jpg
}
/////////////////////////////////////////////////////////////////////

proc computemcm(list Eolist)
"USAGE: computemcm(Eolist); Eolist list
RETURN: an integer, the least common multiple of the denominators of the E-orders
NOTE: Make the same as lcmofall but for one chart. NECESSARY BECAUSE THE E-ORDERS ARE OF TYPE NUMBER!!
EXAMPLE: example computemcm; shows an example
"
{
  int m,i,aux,mcmchart;
  intvec num;

  m=size(Eolist);

  if (m==1)
  {
    mcmchart=int(denominator(Eolist[1]));
    return(mcmchart);
  }

  if (m>1)
  {
    num=int(denominator(Eolist[1]));
    for (i=2;i<=m;i++)
    {
      aux=int(denominator(Eolist[i]));
      num=num,aux;
    }
  }
  mcmchart=lcm(num);
  return(mcmchart);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..2)),dp;
  ideal J=x(1)^3-x(1)*x(2)^3;
  list L=Eresol(J);   // 8 charts, rational exponents
  L[8][2][2];         // maximal E-order at the first chart
  computemcm(L[8][2][2]);
}
/////////////////////////////////////////////////////////////////////

proc constructH(intvec Hhist,int n,list flag)
"USAGE: constructH(Hhist,n,flag);
        Hhist intvec, n integer, flag list
RETURN: the list of exceptional divisors accumulated at this chart
EXAMPLE: example constructH; shows an example
"
{
int i,j,m,l;
list exceplist;
ideal aux;

m=size(Hhist);
if (Hhist[1]==0 and m>1){Hhist=Hhist[2..m]; m=m-1;

                         for (i=1;i<=m;i++){
                                            l=Hhist[i];
                                            if (flag[l]==0){aux=ideal(poly(x(l))); }
                                            else {aux=ideal(poly(y(l))); }

                                            exceplist[i]=aux;
                                            }
// eliminate repeated variables
for (i=1;i<=m;i++){for (j=1;j<=m;j++){
                                      if (Hhist[i]==Hhist[j] and i!=j){
                                                                       if (i<j){exceplist[i]=ideal(1);}
                                                                       if (i>j){exceplist[j]=ideal(1);}
                                                                      }
                                     }
                  }

                         }
else  {exceplist=list();}

// else  {exceplist=list(ideal(0));} // IF IT FAILS USE THIS

return(exceplist);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..3)),dp;
  list flag=identifyvar();
  ideal J=x(1)^4*x(2)^2, x(1)^2+x(3)^3;
  list L=Eresol(J);   // 7 charts
  // history of the exceptional divisors at the 7-th chart
  L[1][7][7]; // blow ups at x(3)-th, x(1)-th and x(1)-th charts
  constructH(L[1][7][7],3,flag);

}
/////////////////////////////////////////////////////////////////////

proc constructblwup(list blwhist,int n,ideal chy,list flag)
"USAGE: constructblwup(blwhist,n,chy,flag);
        blwhist, flag lists, n integer, chy ideal
RETURN: the ideal defining the map K[W] --> K[Wi],
        which gives the composition map of all the blowing up leading to this chart
NOTE: NECESSARY START WITH COLUMNS
EXAMPLE: example constructblwup; shows an example
"
{
int i,j,m,m2;
poly aux2;

m=size(blwhist[1]);

for (j=1;j<=m;j++){
                   for (i=1;i<=n;i++){ m2=blwhist[i][j];

// If m2!=0 this variable changes. First decide if the variable to multiply is invertible or not

                                       if  (m2!=0){
                                                  if (flag[m2]==0){aux2=poly(x(m2));}
                                                  else {aux2=poly(y(m2));}

// And then substitute this variable for the corresponding product in the whole ideal

                                                  if (flag[i]==0){chy=subst(chy,x(i),x(i)*aux2);}
                                                  else {chy=subst(chy,y(i),y(i)*aux2);}

                                                  }
                                      }

                   }

return(chy);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..3)),dp;
  list flag=identifyvar();
  ideal chy=maxideal(1);
  ideal J=x(1)^4*x(2)^2, x(1)^2+x(3)^3;
  list L=Eresol(J);   // 7 charts
  // history of the blow ups at the 7-th chart, center {x(1)=x(3)=0} every time
  L[1][7][8]; // blow ups at x(3)-th, x(1)-th and x(1)-th charts
  constructblwup(L[1][7][8],3,chy,flag);
}
/////////////////////////////////////////////////////////////////////

proc constructlastblwup(list blwhist,int n,ideal chy,list flag)
"USAGE: constructlastblwup(blwhist,n,chy,flag);
        blwhist, flag lists, n integer, chy ideal
RETURN: the ideal defining the last blow up
NOTE: NECESSARY START WITH COLUMNS
EXAMPLE: example constructlastblwup; shows an example
"
{
int i,j,m,m2;
poly aux2;

m=size(blwhist[1]);

if (m>0){
for (i=1;i<=n;i++){ m2=blwhist[i][m];

// If m2!=0 this variable changes. First decide if the variable to multiply is invertible or not

                                       if  (m2!=0){
                                                  if (flag[m2]==0){aux2=poly(x(m2));}
                                                  else {aux2=poly(y(m2));}

// And then substitute this variable for the corresponding product in the whole ideal

                                                  if (flag[i]==0){chy=subst(chy,x(i),x(i)*aux2);}
                                                  else {chy=subst(chy,y(i),y(i)*aux2);}

                                                  }
                  }
}

return(chy);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..3)),dp;
  list flag=identifyvar();
  ideal chy=maxideal(1);
  ideal J=x(1)^4*x(2)^2, x(1)^2+x(3)^3;
  list L=Eresol(J);   // 7 charts
  // history of the blow ups at the 7-th chart, center {x(1)=x(3)=0} every time
  L[1][7][8]; // blow ups at x(3)-th, x(1)-th and x(1)-th charts
  constructlastblwup(L[1][7][8],3,chy,flag);
}
/////////////////////////////////////////////////////////////////////

proc tradtoideal(intvec a,ideal J2,list flag)
"USAGE: tradtoideal(a,J2,flag);
        a intvec, J2 ideal, flag list
COMPUTE: traslate to an ideal the intvec defining the center
RETURN: the ideal of the center, given by the intvec a, or J2 if a=0
EXAMPLE: example tradtoideal; shows an example
"
{
int i,m;
ideal acenter,aux2;

if (a==0){acenter=J2;}
else{
     m=size(a);
     for (i=1;i<=m;i++){
                        if (flag[a[i]]==0){aux2=poly(x(a[i]));}
                        else {aux2=poly(y(a[i]));}

                        acenter=acenter+aux2;
                        }
    }
return(acenter);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..3)),dp;
  list flag=identifyvar();
  ideal J=x(1)^4*x(2)^2, x(1)^2+x(3)^3;
  intvec a=1,3; // first center of blowing up
  tradtoideal(a,J,flag);
}
//////////////////////////////////////////////////////////////////////////////////////
//              OPERATIONS WITH LISTS
//////////////////////////////////////////////////////////////////////////////////////

 proc iniD(int n)
"USAGE: iniD(n);   n integer
RETURN: list of lists of zeros of size n
EXAMPLE: example iniD; shows an example
"
 {int i,j;
  list D,auxD;
 for (j=1;j<=n; j++) {auxD[j]=0;}
 for (i=1;i<=n; i++) {D[i]=auxD;}
 return(D);
 }
example
{"EXAMPLE:"; echo = 2;
  iniD(3);
}
/////////////////////////////////////////////////////////

proc sumlist(list L1,list L2)
"USAGE: sumlist(L1,L2);   L1,L2 lists, (size(L1)==size(L2))
RETURN: a list, sum of L1 and L2
EXAMPLE: example sumlist; shows an example
"
{
 int i,k;
 list sumL;
k=size(L1);
if (size(L2)!=k) {return("ERROR en sumlist, lists must have the same size");}
for (i=1;i<=k;i++) {sumL[i]=L1[i]+L2[i];}
return(sumL);
}
example
{"EXAMPLE:"; echo = 2;
  list L1=1,2,3;
  list L2=5,9,7;
  sumlist(L1,L2);
}
///////////////////////////////////////////////////////

proc reslist(list L1,list L2)
"USAGE: reslist(L1,L2);  L1,L2 lists, (size(L1)==size(L2))
RETURN: a list, subtraction of L1 and L2
EXAMPLE: example reslist; shows an example
"
{
 int i,k;
 list resL;
k=size(L1);
if (size(L2)!=k) {return("ERROR en reslist, lists must have the same size");}
for (i=1;i<=k;i++) {resL[i]=L1[i]-L2[i];}
return(resL);
}
example
{"EXAMPLE:"; echo = 2;
  list L1=1,2,3;
  list L2=5,9,7;
  reslist(L1,L2);
}
//////////////////////////////////////////////////////

proc multiplylist(list L,number a)
"USAGE: multiplylist(L,a); L list, a number
RETURN: list of elements of type number, multiplication of L times a
EXAMPLE: example multiplylist; shows an example
"
{int i,k;
 list newL,bb;
 number b;
 k=size(L);
 for (i=1;i<=k;i++) {b=L[i]*a; bb=b; newL=newL+bb;}
return(newL);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..3)),dp;
  list L=1,2,3;
  multiplylist(L,1/5);
}
///////////////////////////////////////////////////////

proc dividelist(list L1,list L2)
"USAGE: dividelist(L1,L2);  L1,L2 lists
RETURN: list of elements of type number, division of L1 by L2
EXAMPLE: example dividelist; shows an example
"
{int i,k,k1,k2;
list LL,bb;
number a1,a2,b;
k1=size(L1);
k2=size(L2);
if (k2!=k1) {print("ERROR en dividelist, lists must have the same size");}
if (k1<=k2) {k=k1;}
else {k=k2;}
 for (i=1;i<=k;i++)
  {a1=L1[i]; a2=L2[i]; b=a1/a2; bb=b; LL=LL+bb;}
return(LL);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1..3)),dp;
  list L1=1,2,3;
  list L2=5,9,7;
  dividelist(L1,L2);
}
///////////////////////////////////////////////////////

proc createlist(list L1,list L2)
"USAGE: createlist(L1,L2);  L1,L2 lists, (size(L1)==size(L2))
RETURN: list of lists of two elements, the first one of L1 and the second of L2
EXAMPLE: example createlist; shows an example
"
{
  int i,k;
  list L,aux;
  k=size(L1);
  if (size(L2)!=k)
  {ERROR ("createlist: lists must have the same size");}
  L=list0(k);
  for (i=1;i<=k;i++)
  {
    if (L1[i]!=0)
    {
      aux=L1[i],L2[i]; L[i]=aux;
    }
    else {L=delete(L,i);}
  }
  return(L);
}
example
{"EXAMPLE:"; echo = 2;
  list L1=1,2,3;
  list L2=5,9,7;
  createlist(L1,L2);
}
///////////////////////////////////////////////////////
proc list0(int n)
"USAGE: list0(n); n integer
RETURN: list of n zeros
EXAMPLE: example list0; shows an example
"
{
  int i;
  list L0;
  for (i=1;i<=n;i++) {L0[i]=0;}
  return(L0);
}
example
{"EXAMPLE:"; echo = 2;
  list0(4);
}
////////////////////////////////////////////////////////////

proc Emaxcont(list Coef,list Exp,int k,int n,list flag)
"USAGE: Emaxcont(Coef,Exp,k,n,flag);
        Coef,Exp,flag lists, k,n, integers
        Exp is a list of lists of exponents, k=size(Exp)
COMPUTE: Identify ALL the variables of E-maximal contact
RETURN: a list with the indexes of the variables of E-maximal contact
EXAMPLE: example Emaxcont; shows an example
"
{
  int i,j,lon;
  number maxEo;
  list L,sums,bx,maxvar;

  L=maxEord(Coef,Exp,k,n,flag);

  maxEo=L[1];
  sums=L[2];

  if (maxEo>0)
  {
    for (i=1;i<=k; i++)
    {
      lon=size(sums[i]);
      if (lon==2)
      {
        if (sums[i][1]==maxEo)     // variables of the first term
        {
          for (j=1;j<=n; j++)
          {
            if(Exp[i][1][j]!=0 and flag[j]==0)
            {
              bx=j; maxvar=maxvar + bx;
            }
          }
        }
        if (sums[i][2]==maxEo)     // variables of the second term
        {
          for (j=1;j<=n; j++)
          {
            if(Exp[i][2][j]!=0 and flag[j]==0){bx=j; maxvar=maxvar + bx;}
          }
        }
      }
      else
      {
        if (sums[i][1]==maxEo)
        {
          for (j=1;j<=n; j++)
          {
            if(Exp[i][1][j]!=0 and flag[j]==0)
            {
              bx=j; maxvar=maxvar + bx;
            }
          }
        }
      }
    }
  }
  else {maxvar=list();}

   // eliminating repeated terms
  maxvar=elimrep(maxvar);
   // It is necessary to check if flag[j]==0 in order to avoid the selection of y variables
  return(maxEo,maxvar);
}
example
{"EXAMPLE:"; echo = 2;
  ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp;
  list flag=identifyvar();
  ideal J=x(1)^3*x(3)-y(2)*y(4)^2,x(5)*y(2)-x(7)*y(4)^2,x(6)^2*(1-y(4)*y(8)^5),x(7)^4*y(8)^2;
  list L=data(J,4,8);
  list hyp=Emaxcont(L[1],L[2],4,8,flag);
  hyp[1]; // max E-order=0
  hyp[2]; // There are no hypersurfaces of E-maximal contact

  ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp;
  list flag=identifyvar();
  ideal J=x(1)^3*x(3)-y(2)*y(4)^2*x(3),x(5)*y(2)-x(7)*y(4)^2,x(6)^2*(1-y(4)*y(8)^5),x(7)^4*y(8)^2;
  list L=data(J,4,8);
  list hyp=Emaxcont(L[1],L[2],4,8,flag);
  hyp[1]; // the E-order is 1
  hyp[2]; // {x(3)=0},{x(5)=0},{x(7)=0} are hypersurfaces of E-maximal contact
}
///////////////////////////////////////////////////////