//////////////////////////////////////////////////////////////////////////////
version="$Id$";
category="General purpose";
info="
LIBRARY:  grobcov.lib   Groebner Cover for parametric ideals.
PURPOSE:  Comprehensive Groebner Systems, Groebner Cover, Canonical Forms.
          The library contains Montes's algorithms to compute the
          canonical Groebner cover of a parametric ideal as described in
          the paper:

          Montes A., Wibmer M.,
          Groebner Bases for Polynomial Systems with parameters.
          Journal of Symbolic Computation 45 (2010) 1391-1425.

          The central routine is grobcov. Given a parametric
          ideal, grobcov outputs its canonical Groebner cover, consisting
          of a set of pairs of (basis, segment). The basis (after
          normalization) is the reduced Groebner basis for each point
          of the segment. The segments are disjoint, locally closed
          and correspond to constant lpp (leading power product)
          of the basis, and are represented in canonical prime
          representation. The segments are disjoint and cover the
          whole parameter space. The output is canonical, it only
          depends on the given parametric ideal and the monomial order.
          This is much more than a simple comprehensive Groebner system.
          The algorithm grobcov allows options to solve partially the
          problem when the whole automatic algorithm does not finish
          in reasonable time.

          grobcov uses a first algorithm cgsdr that outputs a disjoint
          reduced comprehensive Groebner system with constant lpp.
          cgsdr can be called directly if only a disjoint reduced
          comprehensive Groebner system is required.

          Two other routines: gencase1 and multigrobcov can be used
          in problems with basis of the generic case equal to 1
          (for example in automatic geometric theorem discovering)
          that allow to obtain partial results even when grobcov does
          not finish in reasonable time.

          For completeness, the library also contains the algorithms
          with similar purposes contained in the old library redcgs.lib.
          These algorithms are, in general, less efficient and do not
          ensure a canonical results, even if they are similar to the
          results obtained with grobcov.
          The old routines are no more recommended and remain in
          this library for didactic purposes. These are
          cgsdrold, grobcovold, buildtreetoMaple, cantreetoMaple.

AUTHORS:  Antonio Montes , Hans Schoenemann.
OVERVIEW: see \"Groebner Bases for Polynomial Systems with parameters\"
          Montes A., Wibmer M.,
          Journal of Symbolic Computation 45 (2010) 1391-1425.
          (http://www-ma2.upc.edu/~montes/).

NOTATIONS: All given and determined polynomials and ideals are in the
@*         basering Q[a][x]; (a=parameters, x=variables)
@*         After defining the ring, the main routines
@*         grobcov, cgsdr, gencase1, multigrobcov
@*         generate the global rings
@*         @R   (Q[a][x]),
@*         @P   (Q[a]),
@*         @RP  (Q[x,a])
@*         that are used inside and killed before the output.
@*         If you want to use some internal routine you must
@*         create before the above rings by calling setglobalrings();
@*         because most of the internal routines use these rings.
@*         The call to the basic routines grobcov, cgsdr, gencase1, multigrobcov
@*         or even the older grobcovold, cgsdrold will kill these rings.

PROCEDURES:

grobcov(F);          Is the basic routine giving the canonical
                     Groebner cover of the parametric ideal F.
                     This routine accepts many options, that
                     allow to obtain results even when the canonical
                     computation does not finish in reasonable time.

cgsdr(F);            Is the procedure for obtaining a first disjoint,
                     reduced comprehensive Groebner system that
                     is used in grobcov, but that can be used
                     independently if only the CGS is required.
                     It is a more efficient version of buildtree
                     that does not output the complete discussion tree
                     but only the terminal vertices giving the
                     disjoint reduced comprehensive Groebner system.

gencase1(F);         Returns the segment of the generic case when his
                     basis is 1. This is useful for automatic discovering
                     of geometrical theorems, as it gives the components
                     where a solution exists and is much more efficient
                     than the complete computation of grobcov.

multigrobcov(F);     In problems like automatic discovery of theorems,
                     when grobcov does not give the answer in reasonable
                     time, and the generic case is expected to
                     have basis 1, one can try with multigrobcov procedure
                     to obtain an answer over the different irreducible
                     components: the generic case with basis 1, and the
                     components not corresponding to the generic case. To
                     deduce from its result the true Groebner cover one
                     must discuss theoretically in which segment
                     must be located the intersecting parts in the
                     different irreducible components.

setglobalrings();    Generates the global rings @R, @P and @PR that are
                     respectively the rings Q[a][x], Q[a], Q[x,a].
                     It is called inside each of the fundamental routines of the
                     library: grobcov, cgsdr, gencase1, multigrobcov, as well as
                     by the old routines cgsdrold, grobcovold and killed
                     before the output.
                     If the user want to use some other internal routine,
                     then setglobalrings() is to be called before, as
                     the rings @R, @P and @RP are needed in most of them.
                     globally, and more internal routines can be used, but
                     These rings are destroyed by the call to any of the basic
                     routines.

pdivi(f,F);          Performs a pseudodivision of a parametric polynomial
                     by a parametric ideal.

pnormalform(f,N,W);  Reduces a parametric polynomial f by a reduced-representation
                     (N,W) of null and non-null conditions over the parameters.
                     Before using it setglobalrings() must be called.

Also included from the old library redcgs.lib the following routines

cgsdrold(F);         Similar to cgsdr using the algorithm buildtree
                     of the old library.
grobcovold(F);       Similar to grobcov with the algorithms of the old
                     library.
buildtreetoMaple(T); Writes into a file the output of cgsdrold called
                     with option ('old',0) into a text file that is Maple
                     readable and can be plotted in Maple using
                     the tplot routine of the library dpgb.
cantreetoMaple(M);   Writes into a text file the output of grobcovold called
                     with  option ('out',1), that is readable
                     in Maple and can be plotted using the routine
                     plotcantree of the Maple library dpgb.

SEE ALSO: compregb_lib
";

LIB "primdec.lib";
LIB "qhmoduli.lib";

// ************ Begin of the grobcov library *********************

// Library grobcov.lib
// (Groebner cover):
// Initial data: 6-9-2009
// Release 1:
// Final data: 30-12-2010
// Contains also the old redcgs.lib library that was created
// Initial data: 21-1-2008
// Release 1:
// Final data: 3-7-2008
// Given and determined polynomials and ideals are in the
// basering Q[a][x];

// ************ Begin of buildtree ******************************

proc setglobalrings()
"USAGE:   setglobalrings();
          No arguments
RETURN:   After its call the rings @R=Q[a][x], @P=Q[a], @RP=Q[x,a] are
          defined as global variables.
NOTE:     It is called internally by the fundamental routines of the
          library grobcov, cgsdr, gencase1, muligrobcov as well as by the
          old ones grobcovold,cgsdrold, and killed before the output.
          The user does not need to call it, except when it is interested
          in using some internal routine of the library that
          uses these rings.
          The basering R, must be of the form Q[a][x], a=parameters,
          x=variables, and should be defined previously.
KEYWORDS: ring, rings
EXAMPLE:  setglobalrings; shows an example"
{
  if (defined(@P)==1)
  {
    kill @P; kill @R; kill @RP;
  }
  def RR=basering;
  def @R=basering;  // must be of the form K[a][x], a=parameters, x=variables
  def Rx=ringlist(RR);
  def @P=ring(Rx[1]);
  list Lx;
  Lx[1]=0;
  Lx[2]=Rx[2]+Rx[1][2];
  Lx[3]=Rx[1][3];
  Lx[4]=Rx[1][4];
  Rx[1]=0;
  def D=ring(Rx);
  def @RP=D+@P;
  exportto(Top,@R);      // global ring K[a][x]
  exportto(Top,@P);      // global ring K[a]
  exportto(Top,@RP);     // global ring K[x,a] with product order
  setring(RR);
}
example
{ "EXAMPLE:"; echo = 2;
  ring R=(0,a,b),(x,y,z),dp;
  setglobalrings();
  @R;
  @P;
  @RP;
}

//*************Auxilliary routines**************

// cld : clears denominators of an ideal and normalizes to content 1
//       can be used in @R or @P or @RP
// input:
//    ideal J (J can be also poly), but the output is an ideal;
// output:
//    ideal Jc (the new form of ideal J without denominators and
//       normalized to content 1)
static proc cld(ideal J)
{
  if (size(J)==0){return(ideal(0));}
  def RR=basering;
  setring(@RP);
  def Ja=imap(RR,J);
  ideal Jb;
  if (size(Ja)==0){return(ideal(0));}
  int i;
  def j=0;
  for (i=1;i<=ncols(Ja);i++){if (size(Ja[i])!=0){j++; Jb[j]=cleardenom(Ja[i]);}}
  setring(RR);
  def Jc=imap(@RP,Jb);
  return(Jc);
}

static proc memberpos(f,J)
//"USAGE:  memberpos(f,J);
//         (f,J) expected (polynomial,ideal)
//               or       (int,list(int))
//               or       (int,intvec)
//               or       (intvec,list(intvec))
//               or       (list(int),list(list(int)))
//               or       (ideal,list(ideal))
//               or       (list(intvec),  list(list(intvec))).
//         The ring can be @R or @P or @RP or any other.
//RETURN:  The list (t,pos) t int; pos int;
//         t is 1 if f belongs to J and 0 if not.
//         pos gives the position in J (or 0 if f does not belong).
//EXAMPLE: memberpos; shows an example"
{
  int pos=0;
  int i=1;
  int j;
  int t=0;
  int nt;
  if (typeof(J)=="ideal"){nt=ncols(J);}
  else{nt=size(J);}
  if ((typeof(f)=="poly") or (typeof(f)=="int"))
  { // (poly,ideal)  or
    // (poly,list(poly))
    // (int,list(int)) or
    // (int,intvec)
    i=1;
    while(i<=nt)
    {
      if (f==J[i]){return(list(1,i));}
      i++;
    }
    return(list(0,0));
  }
  else
  {
    if ((typeof(f)=="intvec") or ((typeof(f)=="list") and (typeof(f[1])=="int")))
    { // (intvec,list(intvec)) or
      // (list(int),list(list(int)))
      i=1;
      t=0;
      pos=0;
      while((i<=nt) and (t==0))
      {
        t=1;
        j=1;
        if (size(f)!=size(J[i])){t=0;}
        else
        {
          while ((j<=size(f)) and t)
          {
            if (f[j]!=J[i][j]){t=0;}
            j++;
          }
        }
        if (t){pos=i;}
        i++;
      }
      if (t){return(list(1,pos));}
      else{return(list(0,0));}
    }
    else
    {
      if (typeof(f)=="ideal")
      { // (ideal,list(ideal))
        i=1;
        t=0;
        pos=0;
        while((i<=nt) and (t==0))
        {
          t=1;
          j=1;
          if (ncols(f)!=ncols(J[i])){t=0;}
          else
          {
            while ((j<=ncols(f)) and t)
            {
              if (f[j]!=J[i][j]){t=0;}
              j++;
            }
          }
          if (t){pos=i;}
          i++;
        }
        if (t){return(list(1,pos));}
        else{return(list(0,0));}
      }
      else
      {
        if ((typeof(f)=="list") and (typeof(f[1])=="intvec"))
        { // (list(intvec),list(list(intvec)))
          i=1;
          t=0;
          pos=0;
          while((i<=nt) and (t==0))
          {
            t=1;
            j=1;
            if (size(f)!=size(J[i])){t=0;}
            else
            {
              while ((j<=size(f)) and t)
              {
                if (f[j]!=J[i][j]){t=0;}
                j++;
              }
            }
            if (t){pos=i;}
            i++;
          }
          if (t){return(list(1,pos));}
          else{return(list(0,0));}
        }
      }
    }
  }
}
//example
//{ "EXAMPLE:"; echo = 2;
//  list L=(7,4,5,1,1,4,9);
//  memberpos(1,L);
//  >
//}


static proc subset(J,K)
//"USAGE:   subset(J,K);
//          (J,K)  expected (ideal,ideal)
//                   or     (list, list)
//RETURN:   1 if all the elements of J are in K, 0 if not.
//EXAMPLE:  subset; shows an example;"
{
  int i=1;
  int nt;
  if (typeof(J)=="ideal"){nt=ncols(J);}
  else{nt=size(J);}
  if (size(J)==0){return(1);}
  while(i<=nt)
  {
    if (memberpos(J[i],K)[1]){i++;}
    else {return(0);}
  }
  return(1);
}
//example
//{ "EXAMPLE:"; echo = 2;
//  list J=list(7,3,2);
//  list K=list(1,2,3,5,7,8);
//  subset(J,K);
//}

// elimintfromideal: elimine the constant numbers from an ideal
//     (designed for W, nonnull conditions)
// input: ideal J
// output:ideal K with the elements of J that are non constants, in the ring @P
static proc elimintfromideal(ideal J)
{
  int i;
  int j=0;
  ideal K;
  if (size(J)==0){return(ideal(0));}
  for (i=1;i<=ncols(J);i++){if (size(variables(J[i])) !=0){j++; K[j]=J[i];}}
  return(K);
}

// simpqcoeffs : simplifies a quotient of two polynomials
// input: two coeficients (or terms), that are considered as a quotient
// output: the two coeficients reduced without common factors
static proc simpqcoeffs(poly n,poly m)
{
  def nc=content(n);
  def mc=content(m);
  def gc=gcd(nc,mc);
  ideal s=n/gc,m/gc;
  return (s);
}

// pdivi : pseudodivision of a poly f by an ideal F in a parametric ideal
//         Q[a][x]
// input:
//   poly f0  (in the parametric ring @R)
//   ideal F0 (in the parametric ring @R)
// output:
//   list (poly r, ideal q, poly mu)
proc pdivi(poly f,ideal F)
"USAGE:   pdivi(f,F);
          poly f: the polynomial to be divided
          ideal F: the divisor ideal
RETURN:   A list (poly r, ideal q, poly m). r is the remainder of the
          pseudodivision, q is the set of quotients, and m is the
          factor by which f is to be multiplied.
NOTE:     pseudodivision of a poly f by an ideal F in @R. Returns a
          list (r,q,m) such that m*f=r+sum(q.G), and no lpp of a divisor
          divides a pp of r.
KEYWORDS: division, reduce
EXAMPLE:  pdivi; shows an example"
{
  int i;
  int j;
  poly r=0;
  poly mu=1;
  def p=f;
  ideal q;
  for (i=1; i<=size(F); i++){q[i]=0;}
  ideal lpf;
  ideal lcf;
  for (i=1;i<=size(F);i++){lpf[i]=leadmonom(F[i]);}
  for (i=1;i<=size(F);i++){lcf[i]=leadcoef(F[i]);}
  poly lpp;
  poly lcp;
  poly qlm;
  poly nu;
  poly rho;
  int divoc=0;
  ideal qlc;
  while (p!=0)
  {
    i=1;
    divoc=0;
    lpp=leadmonom(p);
    lcp=leadcoef(p);
    while (divoc==0 and i<=size(F))
    {
      qlm=lpp/lpf[i];
      if (qlm!=0)
      {
        qlc=simpqcoeffs(lcp,lcf[i]);
        nu=qlc[2];
        mu=mu*nu;
        rho=qlc[1]*qlm;
        p=nu*p-rho*F[i];
        r=nu*r;
        for (j=1;j<=size(F);j++){q[j]=nu*q[j];}
        q[i]=q[i]+rho;
        divoc=1;
      }
      else {i++;}
    }
    if (divoc==0)
    {
      r=r+lcp*lpp;
      p=p-lcp*lpp;
    }
  }
  list res=r,q,mu;
  return(res);
}
example
{ "EXAMPLE:"; echo = 2;
  ring R=(0,a,b,c),(x,y),dp;
  "Divisor=";
  poly f=(ab-ac)*xy+(ab)*x+(5c);
  "Dividends=";
  ideal F=ax+b,cy+a;
  "(Remainder, quotients, factor)=";
  def r=pdivi(f,F);
  r;
  "Verifying the division: r[3]*f-(r[2][1]*F[1]+r[2][2]*F[2])-r[1] =";
  r[3]*f-(r[2][1]*F[1]+r[2][2]*F[2])-r[1];
}

// pspol : S-poly of two polynomials in @R
// @R
// input:
//   poly f  (given in the ring @R)
//   poly g (given in the ring @R)
// output:
//   list (S, red):  S is the S-poly(f,g) and red is a Boolean variable
//                if red==1 then S reduces by Buchberger 1st criterion (not used)
static proc pspol(poly f,poly g)
{
  def lcf=leadcoef(f);
  def lcg=leadcoef(g);
  def lpf=leadmonom(f);
  def lpg=leadmonom(g);
  def v=gcd(lpf,lpg);
  def s=simpqcoeffs(lcf,lcg);
  def vf=lpf/v;
  def vg=lpg/v;
  poly S=s[2]*vg*f-s[1]*vf*g;
  return(S);
}

// facvar: Returns all the free-square factors of the elements
//         of ideal J (non repeated). Integer factors are ignored,
//         even 0 is ignored. It can be called from ideal @R, but
//         the given ideal J must only contain poynomials in the
//         parameters.
//         Operates in the ring @P, but can be called from ring @R,
//         and the ideal @P must be defined calling first setglobalrings();
// input:   ideal J
// output:  ideal Jc: Returns all the free-square factors of the elements
//         of ideal J (non repeated). Integer factors are ignored,
//         even 0 is ignored. It can be called from ideal @R, but
//         the given ideal J must only contain poynomials in the
//         parameters.
static proc facvar(ideal J)
//"USAGE:   facvar(J);
//          J: an ideal in the parameters
//RETURN:   all the free-square factors of the elements
//          of ideal J (non repeated). Integer factors are ignored,
//          even 0 is ignored. It can be called from ideal @R, but
//          the given ideal J must only contain poynomials in the
//          parameters.
//NOTE:     Operates in the ring @P, and the ideal J must contain only
//          polynomials in the parameters, but can be called from ring @R.
//KEYWORDS: factor
//EXAMPLE:  facvar; shows an example"
{
  int i;
  def RR=basering;
  setring(@P);
  def Ja=imap(RR,J);
  if(size(Ja)==0){return(ideal(0));}
  Ja=elimintfromideal(Ja); // also in ideal @P
  ideal Jb;
  if (size(Ja)==0){Jb=ideal(0);}
  else
  {
    for (i=1;i<=ncols(Ja);i++){if(size(Ja[i])!=0){Jb=Jb,factorize(Ja[i],1);}}
    Jb=simplify(Jb,2+4+8);
    Jb=cld(Jb);
    Jb=elimintfromideal(Jb); // also in ideal @P
  }
  setring(RR);
  def Jc=imap(@P,Jb);
  return(Jc);
}
//example
//{ "EXAMPLE:"; echo = 2;
//  ring R=(0,a,b,c),(x,y,z),dp;
//  setglobalrings();
//  ideal J=a2-b2,a2-2ab+b2,abc-bc;
//  facvar(J);
//}

// Wred: eliminate the factors in the polynom f that are in W
//       in ring @RP
// input:
//   poly f:
//   ideal W  of non-null conditions (already supposed that it is facvar)
// output:
//   poly f2  where the non-null conditions in W have been dropped from f
static proc Wred(poly f, ideal W)
{
  if (f==0){return(f);}
  def RR=basering;
  setring(@RP);
  def ff=imap(RR,f);
  def RPW=imap(RR,W);
  def l=factorize(ff,2);
  int i;
  poly f1=1;
  for(i=1;i<=size(l[1]);i++)
  {
    if ((memberpos(l[1][i],RPW)[1]) or (memberpos(-l[1][i],RPW)[1])){;}
    else{f1=f1*((l[1][i])^(l[2][i]));}
  }
  setring(RR);
  def f2=imap(@RP,f1);
  return(f2);
}

// pnormalform: reduces a polynomial wrt a red-spec dividing by N and eliminating factors in W.
//              called in the ring @R
//              operates in the ring @RP
//              both ideals must be defined calling first setglobalrings();
// input:
//         poly  f
//         ideal N  (depends only on the parameters)
//         ideal W  (depends only on the parameters)
//                   (N,W) must be a red-spec (depends only on the parameters)
// output: poly f2 reduced wrt to the red-spec (N,W)
// note:   for security a lot of work is done. If (N,W) is already a red-spec
//         it should be simplified
proc pnormalform(poly f, ideal N, ideal W)
"USAGE:   pnormalform(f,N,W);
          f: the polynomial to be reduced modulo (N,W) a reduced representation
          of a segment in the parameters.
          N: the null conditions ideal
          W: the non-null conditions (set of irreducible polynomials)
RETURN:   a reduced polynomial g of f, whose coefficients are reduced
          modulo N and having no factor in W.
NOTE:     Should be called from ring Q[a][x], and the global rings @R, @P
          and @RP must be defined. These rings can be created by calling
          previously setglobalrings();
          Ideals N and W must be given by polynomials
          in the parameters forming a reduced-representation (see
          definition in the paper).
KEYWORDS: division, pdivi, reduce
EXAMPLE:  pnormalform; shows an example"
{
    def RR=basering;
    setglobalrings();
    setring(@RP);
    def fa=imap(RR,f);
    def Na=imap(RR,N);
    def Wa=imap(RR,W);
    option(redSB);
    Na=std(Na);
    def r=cld(reduce(fa,Na));
    def f1=Wred(r[1],Wa);
    setring(RR);
    def f2=imap(@RP,f1);
    return(f2)
}
example
{ "EXAMPLE:"; echo = 2;
  ring R=(0,a,b,c),(x,y),dp;
  setglobalrings();
  poly f=(b^2-1)*x^3*y+(c^2-1)*x*y^2+(c^2*b-b)*x+(a-bc)*y;
  ideal N=(ab-c)*(a-b),(a-bc)*(a-b);
  ideal W=a^2-b^2,bc;
  def r=redspec(N,W);
  pnormalform(f,r[1],r[2]);
}

// idint: ideal intersection
//        in the ring @P.
//        it works in an extended ring
// input: two ideals in the ring @P
// output the intersection of both (is not a GB)
static proc idint(ideal I, ideal J)
{
  def RR=basering;
  ring T=0,t,lp;
  def K=T+RR;
  setring(K);
  def Ia=imap(RR,I);
  def Ja=imap(RR,J);
  ideal IJ;
  int i;
  for(i=1;i<=size(Ia);i++){IJ[i]=t*Ia[i];}
  for(i=1;i<=size(Ja);i++){IJ[size(Ia)+i]=(1-t)*Ja[i];}
  ideal eIJ=eliminate(IJ,t);
  setring(RR);
  return(imap(K,eIJ));
}

// redspec: generates a red-representation
//          called in any ring
//          it changes to the ring @P
//          So the globalrings @P, @RP, @R, must be created before
//          using it by calling setglobalrings();
// input:
//   ideal N : the ideal of null-conditions
//   ideal W : set of non-null polynomials:
//             if W corresponds to no non null conditions then W=ideal(0)
//             otherwise it should be given as an ideal.
// returns: list (Na,Wa,DGN)
// the completely reduced representation:
//   Na = ideal reduced and radical of the red-spec
//   facvar(Wa) = ideal the reduced non-null set of polynomials of the red-spec.
//             if it corresponds to no non null conditions then it is ideal(0)
//             otherwise the ideal is returned.
//   DGN = the list of prime ideals associated to Na (uses primASSGTZ in "primdec.lib")
//   none of the polynomials in facvar(Wa) are contained in none of the ideals in DGN
//   If the given conditions are not compatible, then N=ideal(1) and DGN=list(ideal(1))
proc redspec(ideal Ni, ideal Wi)
//"USAGE:   redspec(N,W);
//          N: null conditions ideal
//          W: set of non-null polynomials (ideal)
//RETURN:   a list (N1,W1,L1) containing a red-representation of the segment (N,W).
//          N1 is the radical reduced ideal characterizing the segment.
//          V(N1) is the Zariski closure of the segment (N,W).
//          The segment S=V(N1) \ V(h), where h=prod(w in W1)
//          N1 is uniquely determined and no prime component of N1 contains none of
//          the polynomials in W1. The polynomials in W1 are prime and reduced
//          wrt N1, and are considered non-null on the segment.
//          L1 contains the list of prime components of N1.
//NOTE:     Called from ring @R it works in ring @P, that must be defined
//          by the call to setglobalrings();
//          Used in the old library redcgs.lib.
//KEYWORDS: representation
//EXAMPLE:  redspec; shows an example"
{
  ideal Nc;
  ideal Wc;
  def RR=basering;
  setring(@P);
  def N=imap(RR,Ni);
  def W=imap(RR,Wi);
  ideal Wa;
  ideal Wb;
  if(size(W)==0){Wa=ideal(0);}
     //when there are no non-null conditions then W=ideal(0)
  else
  {
    Wa=facvar(W);
  }
  if (size(N)==0)
  {
    setring(RR);
    Wc=imap(@P,Wa);
    return(list(ideal(0), Wc, list(ideal(0))));
  }
  int i;
  list LNb;
  list LNa;
  def LN=minGTZ(N);
  for (i=1;i<=size(LN);i++)
  {
    option(redSB);
    LNa[i]=std(LN[i]);
  }
  poly h=1;
  if (size(Wa)!=0)
  {
    for(i=1;i<=size(Wa);i++){h=h*Wa[i];}
  }
  ideal Na;
  intvec save_opt=option(get);
  if (size(N)!=0 and (size(LNa)>0))
  {
    option(returnSB);
    Na=intersect(LNa[1..size(LNa)]);
    option(redSB);
    Na=std(Na);
    option(set,save_opt);
  }
  attrib(Na,"isSB",1);
  if (reduce(h,Na,1)==0)
  {
    setring(RR);
    Wc=imap(@P,Wa);
    return(list (ideal(1),Wc,list(ideal(1))));
  }
  i=1;
  while(i<=size(LNa))
  {
    if (reduce(h,LNa[i],1)==0){LNa=delete(LNa,i);}
    else{ i++;}
  }
  if (size(LNa)==0)
  {
    setring(RR);
    return(list(ideal(1),ideal(0),list(ideal(1))));
  }
  option(returnSB);
  ideal Nb=intersect(LNa[1..size(LNa)]);
  option(redSB);
  Nb=std(Nb);
  option(set,save_opt);
  if (size(Wa)==0)
  {
    setring(RR);
    Nc=imap(@P,Nb);
    Wc=imap(@P,Wa);
    LNb=imap(@P,LNa);
    return(list(Nc,Wc,LNb));
  }
  Wb=ideal(0);
  attrib(Nb,"isSB",1);
  for (i=1;i<=size(Wa);i++){Wb[i]=reduce(Wa[i],Nb);}
  Wb=facvar(Wb);
  if (size(LNa)!=0)
  {
    setring(RR);
    Nc=imap(@P,Nb);
    Wc=imap(@P,Wb);
    LNb=imap(@P,LNa);
    return(list(Nc,Wc,LNb))
  }
  else
  {
    setring(RR);
    Nd=imap(@P,Nb);
    Wc=imap(@P,Wb);
    kill LNb;
    list LNb;
    return(list(Nd,Wc,LNb))
  }
}
//example
//{ "EXAMPLE:"; echo = 2;
//  ring r=(0,a,b,c),(x,y),dp;
//  setglobalrings();
//  ideal N=(ab-c)*(a-b),(a-bc)*(a-b);
//  ideal W=a^2-b^2,bc;
//  redspec(N,W);
//}

// lesspol: compare two polynomials by its leading power products
// input:  two polynomials f,g in the ring @R
// output: 0 if f<g,  1 if f>=g
static proc lesspol(poly f, poly g)
{
  if (leadmonom(f)<leadmonom(g)){return(1);}
  else
  {
    if (leadmonom(g)<leadmonom(f)){return(0);}
    else
    {
      if (leadcoef(f)<leadcoef(g)){return(1);}
      else {return(0);}
    }
  }
}

// delfromideal: deletes the i-th polynomial from the ideal F
static proc delfromideal(ideal F, int i)
{
  int j;
  ideal G;
  if (size(F)<i){ERROR("delfromideal was called incorrect arguments");}
  if (size(F)<=1){return(ideal(0));}
  if (i==0){return(F);}
  for (j=1;j<=size(F);j++)
  {
    if (j!=i){G[size(G)+1]=F[j];}
  }
  return(G);
}

// delidfromid: deletes the polynomials in J that are in I
// input: ideals I,J
// output: the ideal J without the polynomials in I
static proc delidfromid(ideal I, ideal J)
{
  int i; list r;
  ideal JJ=J;
  for (i=1;i<=size(I);i++)
  {
    r=memberpos(I[i],JJ);
    if (r[1])
    {
      JJ=delfromideal(JJ,r[2]);
    }
  }
  return(JJ);
}

// sortideal: sorts the polynomials in an ideal by lm in ascending order
static proc sortideal(ideal Fi)
{
  def RR=basering;
  setring(@RP);
  def F=imap(RR,Fi);
  def H=F;
  ideal G;
  int i;
  int j;
  poly p;
  while (size(H)!=0)
  {
    j=1;
    p=H[1];
    for (i=1;i<=size(H);i++)
    {
      if(lesspol(H[i],p)){j=i;p=H[j];}
    }
    G[size(G)+1]=p;
    H=delfromideal(H,j);
  }
  setring(RR);
  def GG=imap(@RP,G);
  return(GG);
}

// mingb: given a basis (gb reducing) it
// order the polynomials is ascending order and
// eliminate the polynomials whose lpp is divisible by some
// smaller one
static proc mingb(ideal F)
{
  int t; int i; int j;
  def H=sortideal(F);
  ideal G;
  if (ncols(H)<=1){return(H);}
  G=H[1];
  for (i=2; i<=ncols(H); i++)
  {
    t=1;
    j=1;
    while (t and (j<i))
    {
      if((leadmonom(H[i])/leadmonom(H[j]))!=0) {t=0;}
      j++;
    }
    if (t) {G[size(G)+1]=H[i];}
  }
  return(G);
}

// redgb: given a minimal basis (gb reducing) it
// reduces each polynomial wrt to the others
static proc redgb(ideal F, ideal N, ideal W)
{
  ideal G;
  ideal H;
  int i;
  if (size(F)==0){return(ideal(0));}
  for (i=1;i<=size(F);i++)
  {
    H=delfromideal(F,i);
    G[i]=pnormalform(pdivi(F[i],H)[1],N,W);
  }
  return(G);
}

//********************Main routines for buildtree******************

// splitspec: a new leading coefficient f is given to a red-spec
//            then splitspec computes the two new red-spec by
//            considering it null, and non null.
// in ring @P
// given f, and the red-spec (N,W)
//     it outputs the null and the non-null red-spec adding f.
//     if some of the output representations has N=1 then
//     there must be no split and buildtree must continue on
//     the compatible red-spec
// input:  poly f coefficient to split if needed
//         list r=(N,W,LN) redspec
// output: list L = list(ideal N0, ideal W0), list(ideal N1, ideal W1), cond
static proc splitspec(poly fi, list ri)
{
  def RR=basering;
  def Ni=ri[1];
  def Wi=ri[2];
  setring(@P);
  def f=imap(RR,fi);
  def N=imap(RR,Ni);
  def W=imap(RR,Wi);
  f=Wred(f,W);
  def N0=N;
  def W1=W;
  N0[size(N0)+1]=f;
  def r0=redspec(N0,W);
  W1[size(W1)+1]=f;
  def r1=redspec(N,W1);
  setring(RR);
  def ra0=imap(@P,r0);
  def ra1=imap(@P,r1);
  def cond=imap(@P,f);
  return (list(ra0,ra1,cond));
}

// redcoefs
// 15/09/2010
static proc redcoefs(poly f, ideal N)
{
  def f1=f; int test0=1; poly lc; poly lm;
  poly lc1;
  def RR=basering;
  setring(@P);
  poly lcp;
  def Np=imap(RR,N);
  attrib(Np,"isSB",1);
  setring(RR);
  while((test0==1) and (f1<>0))
  {
    lc=leadcoef(f1);
    lm=leadmonom(f1);
    setring(@P);
    lcp=imap(RR,lc);
    lcp=reduce(lcp,Np);
    setring(RR);
    lc1=imap(@P,lcp);
    if(lc1<>0){test0=0;}
    f1=f1+(lc1-lc)*lm;
  }
  return(f1);
}

// discusspolys: given a basis B and a red-spec (N,W), it analyzes the
//               leadcoef of the polynomials in B until it finds
//               that one of them can be either null or non null.
//               If at the end only the non null option is compatible
//               then the reduced B has all the leadcoef non null.
//               Else recbuildtree must split.
// ring @R
// input:  ideal B
//         ideal N
//         ideal W (a reduced-representation)
// output: list of ((N0,W0,LN0),(N1,W1,LN1),Br,cond)
//         cond is the condition to branch
static proc discusspolys(ideal B, list r)
{
  poly f;     poly f1;    poly f2;
  poly cond;
  def N=r[1]; def W=r[2]; def LN=r[3];
  def Ba=B;   def F=B;
  ideal N0=1; def W0=W;   list LN0=ideal(1);
  def N1=N;   def W1=W;   def LN1=LN;
  list L;
  list M;     list M0;    list M1;
  list rr;
  if (size(B)==0)
  {
    M0=N0,W0,LN0; // incompatible
    M1=N1,W1,LN1;
    M=M0,M1,B,poly(1);
    return(M);
  }
  while ((size(F)!=0) and ((N0[1]==1) or (N1[1]==1)))
  {
    f=F[1];
    F=delfromideal(F,1);
    f1=pnormalform(f,N,W);
    rr=memberpos(f,Ba);
    if (f1!=0)
    {
      Ba[rr[2]]=f1;
      if (pardeg(leadcoef(f1))!=0)
      {
        f2=Wred(leadcoef(f1),W);
        L=splitspec(f2,list(N,W,LN));
        N0=L[1][1]; W0=L[1][2]; LN0=L[1][3]; N1=L[2][1]; W1=L[2][2]; LN1=L[2][3];
        cond=L[3];
      }
    }
    else
    {
      Ba=delfromideal(Ba,rr[2]);
      N0=ideal(1); //F=ideal(0);
    }
  }
  M0=N0,W0,LN0;
  M1=N1,W1,LN1;
  M=M0,M1,Ba,cond;
  return(M);
}

// discussSpolys: given a basis B and a red-spec (N,W), it analyzes the
//                leadcoef of the polynomials in B until it finds
//                that one of them can be either null or non null.
//                If at the end only the non null option is compatible
//                then the reduced B has all the leadcoef non null.
//                Else recbuildtree must split.
// ring @R
// input:  ideal B
//         ideal N
//         ideal W (a reduced-representation)
//         list  P current set of pairs of polynomials from B to be tested.
// output: list of (N0,W0,LN0),(N1,W1,LN1),Br,Pr,cond]
//         list Pr the not checked list of pairs.
static proc discussSpolys(ideal B, list r, list P)
{
  int i; int j; int k;
  int npols; int nSpols; int tt;
  poly cond=1;
  poly lm; poly lpf; poly lpg;
  def F=B; def Pa=P; list Pa0;
  def N=r[1]; def W=r[2]; def LN=r[3];
  ideal N0=1; def W0=W; list LN0=ideal(1);
  def N1=N; def W1=W; def LN1=LN;
  ideal Bw;
  poly S;
  list L; list L0; list L1;
  list M; list M0; list M1;
  list pair;
  list KK; int loc;
  int crit;
  poly h;
  if (size(B)==0)
  {
    M0=N0,W0,LN0;
    M1=N1,W1,LN1;
    M=M0,M1,ideal(0),Pa,cond;
    return(M);
  }
  tt=1;
  i=1;
  while ((tt) and (i<=size(B)))
  {
    h=B[i];
    for (j=1;j<=npars(@R);j++)
    {
      h=subst(h,par(j),0);
    }
    if (h!=B[i]){tt=0;}
    i++;
  }
  if (tt)
  {
    //"T_ a non parametric system occurred";
    def RR=basering;
    def RL=ringlist(RR);
    RL[1]=0;
    def LRR=ring(RL);
    setring(LRR);
    def BP=imap(RR,B);
    option(redSB);
    BP=std(BP);
    setring(RR);
    B=imap(LRR,BP);
    M0=ideal(1),W0,LN0;
    M1=N1,W1,LN1;
    M=M0,M1,B,list(),cond;
    return(M);
  }
  if (size(Pa)==0){npols=size(B); Pa=orderingpairs(F); nSpols=size(Pa);}
  while ((size(Pa)!=0) and (N0[1]==1) or (N1[1]==1))
  {
    pair=Pa[1];
    i=pair[1];
    j=pair[2];
    Pa=delete(Pa,1);
    // Buchberger 1st criterion (not needed here, it is already eliminated
    // when creating the list of pairs
    for (k=1;k<=size(Pa);k++){Pa0[k]=delete(Pa[k],3);}
    crit=0;
    if (not(crit))
    {
      S=pspol(F[i],F[j]);
      KK=pdivi(S,F);
      S=KK[1];
      if (S!=0)
      {
        S=pnormalform(S,N,W);
        if (S!=0)
        {
          L=discusspolys(ideal(S),list(N,W,LN));
          N0=L[1][1];
          W0=L[1][2];
          LN0=L[1][3];
          N1=L[2][1];
          W1=L[2][2];
          LN1=L[2][3];
          S=L[3][1];
          cond=L[4];
          if (S==1)
          {
            M0=ideal(1),W0,list(ideal(1));
            M1=N1,W1,LN1;
            M=M0,M1,ideal(1),list(),cond;
            return(M);
          }
          if (S!=0)
          {
            F[size(F)+1]=S;
            npols=size(F);
            for (k=1;k<size(F);k++)
            {
              lm=lcmlmonoms(F[k],S);
              // Buchberger 1st criterion
              lpf=leadmonom(F[k]);
              lpg=leadmonom(S);
              if (lpf*lpg!=lm)
              {
                pair=k,size(F),lm;
                Pa=placepairinlist(pair,Pa);
                nSpols=size(Pa);
              }
            }
            if (N0[1]==1){N=N1; W=W1; LN=LN1;}
          }
        }
      }
    }
  }
  M0=N0,W0,LN0;
  M1=N1,W1,LN1;
  M=M0,M1,F,Pa,cond;
  return(M);
}

// lcmlmonoms: computes the lcm of the leading monomials
//             of the polynomils f and g
// ring @R
static proc lcmlmonoms(poly f,poly g)
{
  def lf=leadmonom(f);
  def lg=leadmonom(g);
  def gls=gcd(lf,lg);
  return((lf*lg)/gls);
}

// placepairinlist
// 15/09/2010
// input:  given a new pair of the form (i,j,lmij)
//         and a list of pairs of the same form
// ring @R
// output: it inserts the new pair in ascending order of lmij
static proc placepairinlist(list pair,list P)
{
  list Pr;
  if (size(P)==0){Pr=insert(P,pair); return(Pr);}
  if (pair[3]<P[1][3]){Pr=insert(P,pair); return(Pr);}
  if (pair[3]>=P[size(P)][3]){Pr=insert(P,pair,size(P)); return(Pr);}
  kill Pr;
  list Pr;
  int j;
  int i=1;
  int loc=0;
  while((i<=size(P)) and (loc==0))
  {
    if (pair[3]>=P[i][3]){j=i; i++;}
    else{loc=1; j=i-1;}
  }
  Pr=insert(P,pair,j);
  return(Pr);
}

// orderingpairs:
// input:  ideal F
// output: list of ordered pairs (i,j,lcmij) of F in ascending order of lcmij
//         if a pair verifies Buchberger 1st criterion it is not stored
// ring @R
static proc orderingpairs(ideal F)
{
  int i;
  int j;
  poly lm;
  poly lpf;
  poly lpg;
  list P;
  list pair;
  if (size(F)<=1){return(P);}
  for (i=1;i<=size(F)-1;i++)
  {
    for (j=i+1;j<=size(F);j++)
    {
      lm=lcmlmonoms(F[i],F[j]);
      // Buchberger 1st criterion
      lpf=leadmonom(F[i]);
      lpg=leadmonom(F[j]);
      if (lpf*lpg!=lm)
      {
        pair=(i,j,lm);
        P=placepairinlist(pair,P);
      }
    }
  }
  return(P);
}

// Buchberger 2nd criterion
// input:  integers i,j
//         list P of pairs of the form (i,j) not yet verified
// ring @R
//         not used (it increases time)
static proc criterion(int i, int j, list P, ideal B)
{
  def lcmij=lcmlmonoms(B[i],B[j]);
  int crit=0;
  int k=1;
  list ik; list jk;
  while ((k<=size(B)) and (crit==0))
  {
    if ((k!=i) and (k!=j))
    {
      if (i<k){ik=i,k;} else{ik=k,i;}
      if (j<k){jk=i,k;} else{jk=k,j;}
      if (not((memberpos(ik,P)[1]) or (memberpos(jk,P)[1])))
      {
        if ((lcmij)/leadmonom(B[k])!=0){crit=1;}
      }
    }
    k++;
  }
  return(crit);
}

// buildtree: Basic routine of the old redcgs.lib generating a
//     first reduced CGS
//     it will define the rings @R, @P and @RP as global rings
//     and the list @T a global list that will be killed at the output
// input:  ideal F in ring K[a][x];
// output: list T of lists whose list elements are of the form
//         T[i]=list(list lab, boolean terminal, ideal B, ideal N, ideal W, list of ideals decomp of N,
//              ideal of monomials lpp);
// all the ideals are in the ring K[a][x];
static proc buildtree(ideal F, list #)
//"USAGE:   buildtree(F);
//          F: ideal in Q[a][x] (parameters and variables) to be discussed.
//          It outputs the whole discussion tree to construct the
//          first disjoint reduced CGS. It is the old version of the new
//          cgsdr routine. It remains on the library for didactic purposes
//          and is, in general, less efficient.
//          Also, for some problems where cgsdr does stack, sometimes
//          buildtree is able to obtain the result.
//          The output of buildtree contains the whole information about the discussion
//          process (the whole tree discussion) and can be reduced to
//          somewhat similar to the output of cgsdr after calling
//          setglobalrings(); then applying finalcases and then groupsegments to the
//          output of buidtree. This is automatically done by the routine
//          cgsdrold also contained in the library, that outputs only the
//          CGS like the new cgsdr.
//
//RETURN:   Returns a list T describing the complete discussion tree
//          for obtaining a reduced and disjoint comprehensive
//          Groebner system (CGS) of the ideal F of Q[a][x] with
//          constant leading power products (lpp) of the reduced Groebner
//          basis.
//          The first element of the list is the root, and contains
//            [1] label: intvec(-1)
//            [2] number of children : int
//            [3] the ideal F
//            [4], [5], [6] the red-representation of the segment
//                (null, non-null conditions, prime components of the null
//                conditions) given (as option).
//                ideal (0), ideal (1), list(ideal(0)) is assumed if
//                no optional conditions are given.
//            [7] the set of lpp of ideal F
//            [8] condition that was taken to reach the vertex
//                (poly 1, for the root).
//          The remaining elements of the list represent vertices of the tree:
//          with the same structure:
//            [1] label: intvec (1,0,0,1,...) gives its position in the tree:
//                first branch condition is taken non-null, second null,...
//            [2] number of children (0 if it is a terminal vertex)
//            [3] the specialized ideal with the previous assumed conditions
//                to reach the vertex
//            [4],[5],[6] the red-representation of the segment corresponding
//                to the previous assumed conditions to reach the vertex
//            [7] the set of lpp of the specialized ideal at this stage
//            [8] condition that was taken to reach the vertex from the
//                father's vertex (that was taken non-null if the last
//                integer in the label is 1, and null if it is 0)
//          The terminal vertices form a disjoint partition of the parameter space
//          whose bases specialize to the reduced Groebner basis of the
//          specialized ideal on each point of the segment and preserve
//          the lpp. So they form a disjoint reduced CGS.
//NOTE:     The basering R, must be of the form Q[a][x], a=parameters,
//          x=variables, and should be defined previously. The ideal must
//          be defined on R.
//          The call of finalcases applied to the output of buildtree
//          selects the terminal vertices forming the disjoint and reduced
//          CGS. To obtain the output similar
//          to that of the new cgsdr procedure one can call instead
//          cgsdrold.
//
//          The content of buildtree can be written in a file that is readable
//          by Maple in order to plot its content using buildtreetoMaple;
//          The file written by buildtreetoMaple when is read in a Maple
//          worksheet can be plotted using the dbgb routine tplot;
//
//KEYWORDS: CGS, disjoint, reduced, comprehensive Groebner system
//EXAMPLE:  buildtree; shows an example"
{
  list @T;
  exportto(Top,@T);
  setglobalrings();
  int i;
  int j;
  poly f;
  poly cond=1;
  list LN;
  LN[1]=ideal(0);
  def N=ideal(0);
  def W=ideal(1);
  int comment=0;
  list L=#;
  for(i=1;i<=size(L)/2;i++)
  {
    if(L[2*i-1]=="null"){N=L[2*i];}
    else
    {
      if(L[2*i-1]=="nonnull"){W=L[2*i];}
      else
      {
        if(L[2*i-1]=="comment"){comment=L[2*i];}
      }
    }
  }
  ideal B;
  if(equalideals(N,ideal(0))==0)
  {
    def LL=redspec(N,W);
    N=LL[1];
    W=LL[2];
    LN=LL[3];
    for (i=1;i<=size(F);i++)
    {
      f=pnormalform(F[i],N,W);
      if (f!=0){B[size(B)+1]=f;}
    }
  }
  else {B=F;}
  def lpp=ideal(0);
  if (size(B)==0){lpp=ideal(0);}
  else
  {
     for (i=1;i<=size(B);i++){lpp[i]=leadmonom(B[i]);}
    // lpp=ideal of lead power product of the polys in B
  }
  intvec lab=-1;
  int term=0;
  list root;
  root[1]=lab;
  root[2]=term;
  root[3]=B;
  root[4]=N;
  root[5]=W;
  root[6]=LN;
  root[7]=lpp;
  root[8]=cond;
  @T[1]=root;
  list P;
  recbuildtree(root,P);
  def T=@T;
  kill @T;
  kill @P; kill @RP; kill @R;
  return(T)
}
//example
//{ "EXAMPLE:"; echo = 2;
//  ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp;
//  "Casas conjecture for degree 4";
//  ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0),
//          x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1),
//          x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0),
//          x2^2+(2*a3)*x2+(a2),
//          x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0),
//          x3+(a3);
//  def T=buildtree(F); "buildtree(F)="; T;
//  setglobalrings();
//  def FC=finalcases(T);
//  "finalcases(buildtree(F))="; FC;
//  "groupsegments(finalcases(buildtree(F)))=";
//  groupsegments(FC);
//  buildtreetoMaple(T,"Tb","Tb.txt"); " ";
//  "Compare with cgsdrold"; " ";
//  def CDR=cgsdrold(F);
//  "cgsdrold(F)="; CDR;
//}

// recbuildtree: auxilliary recursive routine called by buildtree
static proc recbuildtree(list v, list P)
{
  def vertex=v;
  int i;
  int j;
  int pos;
  list P0;
  list P1;
  poly f;
  def lab=vertex[1];
  if ((size(lab)>1) and (lab[1]==-1))
  {lab=lab[2..size(lab)];}
  def term=vertex[2];
  def B=vertex[3];
  def N=vertex[4];
  def W=vertex[5];
  def LN=vertex[6];
  def lpp=vertex[7];
  def cond=vertex[8];
  def lab0=lab;
  def lab1=lab;
  if ((size(lab)==1) and (lab[1]==-1))
  {
    lab0=0;
    lab1=1;
  }
  else
  {
    lab0[size(lab)+1]=0;
    lab1[size(lab)+1]=1;
  }
  list vertex0;
  list vertex1;
  ideal B0;
  ideal lpp0;
  ideal lpp1;
  ideal N0=1;
  def W0=ideal(0);
  list LN0=ideal(1);
  def B1=B;
  def N1=N;
  def W1=W;
  list LN1=LN;
  list L;
  if (size(P)==0)
  {
    L=discusspolys(B,list(N,W,LN));
    N0=L[1][1];
    W0=L[1][2];
    LN0=L[1][3];
    N1=L[2][1];
    W1=L[2][2];
    LN1=L[2][3];
    B1=L[3];
    cond=L[4];
  }
  if ((size(B1)!=0) and (N0[1]==1))
  {
    L=discussSpolys(B1,list(N1,W1,LN1),P);
    N0=L[1][1];
    W0=L[1][2];
    LN0=L[1][3];
    N1=L[2][1];
    W1=L[2][2];
    LN1=L[2][3];
    B1=L[3];
    P1=L[4];
    cond=L[5];
    lpp=ideal(0);
    for (i=1;i<=size(B1);i++){lpp[i]=leadmonom(B1[i]);}
  }
  vertex[3]=B1;
  vertex[4]=N1; // unnecessary
  vertex[5]=W1; // unnecessary
  vertex[6]=LN1;// unnecessary
  vertex[7]=lpp;
  vertex[8]=cond;
  if (size(@T)>0)
  {
    pos=size(@T)+1;
    @T[pos]=vertex;
  }
  if ((N0[1]!=1) and (N1[1]!=1))
  {
    vertex1[1]=lab1;
    vertex1[2]=0;
    vertex1[3]=B1;
    vertex1[4]=N1;
    vertex1[5]=W1;
    vertex1[6]=LN1;
    vertex1[7]=lpp1;
    vertex1[8]=cond;
    if (size(B1)==0){B0=ideal(0); lpp0=ideal(0);}
    else
    {
      j=1;
      lpp0=ideal(0);
      for (i=1;i<=size(B1);i++)
      {
        f=pnormalform(B1[i],N0,W0);
        if (f!=0){B0[j]=f; lpp0[j]=leadmonom(f);j++;}
      }
    }
    vertex0[1]=lab0;
    vertex0[2]=0;
    vertex0[3]=B0;
    vertex0[4]=N0;
    vertex0[5]=W0;
    vertex0[6]=LN0;
    vertex0[7]=lpp0;
    vertex0[8]=cond;
    recbuildtree(vertex0,P0);
    recbuildtree(vertex1,P1);
  }
  else
  {
    if (equalideals(N1,ideal(1))==0)
    {
      vertex[2]=1;
      B1=mingb(B1);
      vertex[3]=redgb(B1,N1,W1);
      vertex[4]=N1;
      vertex[5]=W1;
      vertex[6]=LN1;
      lpp=ideal(0);
      for (i=1;i<=size(vertex[3]);i++){lpp[i]=leadmonom(vertex[3][i]);}
      vertex[7]=lpp;
      vertex[8]=cond;
      @T[pos]=vertex;
      //print(vertex);
    }
  }
}

// RtoPrep
// Computes the P-representaion of a R-representaion (N,W,L) of a set
// input:
//    ideal N (null conditions, must be radical)
//    ideal W (non-null conditions ideal)
//    list L  must contain the radical decomposition of N.
// output:
//    the ((p_1,(p_11,..,p_1k_1)),..,(p_r,(p_r1,..,p_rk_r)));
//    the Prep of V(N) \ V(h), where h=prod(w in W).
static proc RtoPrep(ideal N, ideal W, list L)
{
  int i; int j; list L0;
  if (N[1]==1)
  {
    L0[1]=list(ideal(1),list(ideal(1)));
    return(L0);
  }
  def RR=basering;
  setring(@P);
  ideal Np=imap(RR,N);
  ideal Wp=imap(RR,W);
  list Lp=imap(RR,L);
  poly h=1;
  for (i=1;i<=size(Wp);i++){h=h*Wp[i];}
  list r; list Ti; list LL;
  for (i=1;i<=size(Lp);i++)
  {
    Ti=minGTZ(Lp[i]+h);
    for(j=1;j<=size(Ti);j++)
    {
      option(redSB);
      Ti[j]=std(Ti[j]);
    }
    //list LL[i];
    LL[i]=list(Lp[i],Ti);
  }
  setring(RR);
  return(imap(@P,LL));
}

// groupRtoPrep
// input:  L (list) is the output of groupsegments
// output: LL (list) the same list but the segments are expressed
//                   in canonical representations:
//  ( (lpp, (lab BuildTree, basis,
//             ((P_1),(P_{11},...,P_{1t1}))
//             ...
//             ((P_j),(P_{j1},...,P_{jtj}))
//          )
//          ...
//          (lab BuildTree, basis,
//             ((P_1),(P_{11},...,P_{1t1}))
//             ...
//             ((P_j),(P_{j1},...,P_{jtj}))
//          )
//    )
//    ...
//    (lpp, (lab BuildTree, basis,
//             ((P_1),(P_{11},...,P_{1t1}))
//             ...
//             ((P_j),(P_{j1},...,P_{jtj}))
//          )
//          ...
//          (lab BuildTree, basis,
//             ((P_1),(P_{11},...,P_{1t1}))
//             ...
//             ((P_j),(P_{j1},...,P_{jtj}))
//          )
//    )
//  )
static proc groupRtoPrep(list L)
{
  int i; int j;
  list LL; list ct;
  // size(L)=number of lpp-segments
  for (i=1;i<=size(L);i++)
  {
    LL[i]=list();
    LL[i][1]=L[i][1];
    // L[i][1]=lpp
    LL[i][2]=list();
    for (j=1;j<=size(L[i][2]);j++)
    {
      ct=RtoPrep(L[i][2][j][3],L[i][2][j][4],L[i][2][j][5]);
      LL[i][2][j]=list();
      LL[i][2][j][1]=L[i][2][j][1];
      // L[i][2][j][1]=label
      LL[i][2][j][2]=L[i][2][j][2];
      // L[i][2][j][2]=basis
      LL[i][2][j][3]=ct;
    }
  }
  return(LL);
}

// NEW
// input:  L (list) is the output of groupsegments
// output: LL (list) the same list but the segments are expressed
//                   in canonical representations:
//  ( (lpp, (lab BuildTree, basis,
//             ((1,u1),(lab,child,P_1)),
//             ((1,1,1),(lab,child,P_{11})),
//             ...
//             ((1,1,t1),(lab,child,P_{1t1})),
//             ...
//             ((1,u1),(lab,child,P_u1)),
//             ((1,u1,1),(lab,child,P_{u1,1})),
//             ...
//             ((1,u1,tu),(lab,child,P_{u1,tu})),
//          (lab BuildTree, basis,
//             ((1,u2),(lab,child,P_2)),
//             ((1,u1+1,1),(lab,child,P_{21})),
//             ...
//             ((1,u1+1,t2),(lab,child,P_{2,t2})),
//             ...
//             ((1,u1+..+ut),(lab,child,P_ut)),
//             ((1,u1+..+ut,1),(lab,child,P_{ut,1})),
//             ...
//             ((1,u1+..+ut,tu),(lab,child,P_{ut,tu})),
// ...
static proc groupredtocan(list L)
{
  int i; int j;
  list LL; list ct;
  for (i=1;i<=size(L);i++)
  {
    LL[i]=list();
    LL[i][1]=L[i][1];
    LL[i][2]=list();
    for (j=1;j<=size(L[i][2]);j++)
    {
      ct=redtocanspec(intvec(i),j-1,list(L[i][2][j][3],L[i][2][j][4],L[i][2][j][5]));
      LL[i][2][j]=list();
      LL[i][2][j][1]=L[i][2][j][1];
      LL[i][2][j][2]=L[i][2][j][2];
      LL[i][2][j][3]=ct;
    }
  }
  return(LL);
}

//****************End of BuildTree*************************************

//****************Begin BuildTree To Maple*****************************

// buildtreetoMaple: writes the list provided by buildtree to a file
//    containing the table representing it in Maple

// writes the list L=buildtree(F) to a file "writefile" that
// is readable by Maple whith name T
// input:
//   L: the list output by buildtree
//   T: the name (string) of the output table in Maple
//   writefile: the name of the datafile where the output is to be stored
// output:
//   the result is written on the datafile "writefile" containig
//   the assignement to the table with name "T"
proc buildtreetoMaple(list L, string T, string writefile)
"USAGE:   buildtreetoMaple(T, TM, writefile);
          T: is the list provided by grobcovold called with option "old",0;
          TM: is the name (string) of the table variable in Maple that will represent
          the output of cgsdrold;
          writefile: is the name (string) of the file whereas to write the
          content.
RETURN:   writes the list provided by grobcovold called with option "old",0,
          (old buildtree) to a file containing the table representing it in
          Maple.
KEYWORDS: cgsdrold, buildtree, Maple
EXAMPLE:  buildtreetoMaple; shows an example"
{
  def R=basering;
  if(size(T[1])!=8)
  {
    "  'Warning!' cgsdrold must be called with option 'old' set to 0 to be operative";
    return();
  }
  short=0;
  poly cond;
  int i;
  link LLw=":w "+writefile;
  string La=string("table(",T,");");
  write(LLw, La);
  close(LLw);
  link LLa=":a "+writefile;
  def RL=ringlist(R);
  list p=RL[1][2];
  string param=string(p[1]);
  if (size(p)>1)
  {
    for(i=2;i<=size(p);i++){param=string(param,",",p[i]);}
  }
  list v=RL[2];
  string vars=string(v[1]);
  if (size(v)>1)
  {
    for(i=2;i<=size(v);i++){vars=string(vars,",",v[i]);}
  }
  list xord;
  list pord;
  if (RL[1][3][1][1]=="dp"){pord=string("tdeg(",param);}
  if (RL[1][3][1][1]=="lp"){pord=string("plex(",param);}
  if (RL[3][1][1]=="dp"){xord=string("tdeg(",vars);}
  if (RL[3][1][1]=="lp"){xord=string("plex(",vars);}
  write(LLa,string(T,"[[9]]:=",xord,");"));
  write(LLa,string(T,"[[10]]:=",pord,");"));
  write(LLa,string(T,"[[11]]:=true; "));
  list S;
  for (i=1;i<=size(L);i++)
  {
    if (L[i][2]==0)
    {
      cond=L[i][8];
      S=btcond(T,L[i],cond);
      write(LLa,S[1]);
      write(LLa,S[2]);
    }
    S=btbasis(T,L[i]);
    write(LLa,S);
    S=btN(T,L[i]);
    write(LLa,S);
    S=btW(T,L[i]);
    write(LLa,S);
    if (L[i][2]==1) {S=btterminal(T,L[i]); write(LLa,S);}
    S=btlpp(T,L[i]);
    write(LLa,S);
  }
  close(LLa);
}
example
{ "EXAMPLE:"; echo = 2;
  ring R=(0,a1,a2,a3,a4),(x1,x2,x3,x4),dp;
  ideal F=x4-a4+a2,
   x1+x2+x3+x4-a1-a3-a4,
   x1*x3*x4-a1*a3*a4,
   x1*x3+x1*x4+x2*x3+x3*x4-a1*a4-a1*a3-a3*a4;
  def T=cgsdrold(F,"old",0); "T="; T;
  buildtreetoMaple(T,"Tb","Tb.txt");
}

// auxiliary routine called by buildtreetoMaple
// input:
//   list L: element i of the list of buildtree(F)
// output:
//   the string of T[[lab,1]]:=label; in Maple
static proc btterminal(string T, list L)
{
  int i;
  string Li;
  string term;
  string coma=",";
  if (L[2]==0){term="false";} else {term="true";}
  def lab=L[1];
  string slab;
  if ((size(lab)==1) and lab[1]==-1)
  {slab="";coma="";} //if (size(lab)==0)
  else
  {
    slab=string(lab[1]);
    if (size(lab)>=1)
    {
      for (i=2;i<=size(lab);i++){slab=string(slab,",",lab[i]);}
    }
  }
  Li=string(T,"[[",slab,coma,"6]]:=",term,"; ");
  return(Li);
}

// auxiliary routine called by buildtreetoMaple
// input:
//   list L: element i of the list of buildtree(F)
// output:
//   the string of T[[lab,3]] (basis); in Maple
static proc btbasis(string T, list L)
{
  int i;
  string Li;
  string coma=",";
  def lab=L[1];
  string slab;
  if ((size(lab)==1) and lab[1]==-1)
  {slab="";coma="";} //if (size(lab)==0)
  else
  {
    slab=string(lab[1]);
    if (size(lab)>=1)
    {
      for (i=2;i<=size(lab);i++){slab=string(slab,",",lab[i]);}
    }
  }
  Li=string(T,"[[",slab,coma,"3]]:=[",L[3],"]; ");
  return(Li);
}

// auxiliary routine called by buildtreetoMaple
// input:
//   list L: element i of the list of buildtree(F)
// output:
//   the string of T[[lab,4]] (null conditions ideal); in Maple
static proc btN(string T, list L)
{
  int i;
  string Li;
  string coma=",";
  def lab=L[1];
  string slab;
  if ((size(lab)==1) and lab[1]==-1)
  {slab=""; coma="";}
  else
  {
    slab=string(lab[1]);
    if (size(lab)>=1)
    {
      for (i=2;i<=size(lab);i++){slab=string(slab,",",lab[i]);}
    }
  }
  if ((size(lab)==1) and lab[1]==-1)
    {Li=string(T,"[[",slab,coma,"4]]:=[ ]; ");}
  else
    {Li=string(T,"[[",slab,coma,"4]]:=[",L[4],"]; ");}
  return(Li);
}

// auxiliary routine called by buildtreetoMaple
// input:
//   list L: element i of the list of buildtree(F)
// output:
//   the string of T[[lab,5]] (null conditions ideal); in Maple
static proc btW(string T, list L)
{
  int i;
  string Li;
  string coma=",";
  def lab=L[1];
  string slab;
  if ((size(lab)==1) and lab[1]==-1)
  {slab=""; coma="";}
  else
  {
    slab=string(lab[1]);
    if (size(lab)>=1)
    {
      for (i=2;i<=size(lab);i++){slab=string(slab,",",lab[i]);}
    }
  }
  if (size(L[5])==0)
    {Li=string(T,"[[",slab,coma,"5]]:={ }; ");}
  else
    {Li=string(T,"[[",slab,coma,"5]]:={",L[5],"}; ");}
  return(Li);
}

// auxiliary routine called by buildtreetoMaple
// input:
//   list L: element i of the list of buildtree(F)
// output:
//   the string of T[[lab,12]] (lpp); in Maple
static proc btlpp(string T, list L)
{
  int i;
  string Li;
  string coma=",";;
  def lab=L[1];
  string slab;
  if ((size(lab)==1) and lab[1]==-1)
  {slab=""; coma="";}
  else
  {
    slab=string(lab[1]);
    if (size(lab)>=1)
    {
      for (i=2;i<=size(lab);i++){slab=string(slab,",",lab[i]);}
    }
  }
  if (size(L[7])==0)
  {
    Li=string(T,"[[",slab,coma,"12]]:=[ ]; ");
  }
  else
  {
    Li=string(T,"[[",slab,coma,"12]]:=[",L[7],"]; ");
  }
  return(Li);
}

// auxiliary routine called by buildtreetoMaple
// input:
//   list L: element i of the list of buildtree(F)
// output:
//   the list of strings of (T[[lab,0]]=0,T[[lab,1]]<>0); in Maple
static proc btcond(string T, list L, poly cond)
{
  int i;
  string Li1;
  string Li2;
  def lab=L[1];
  string slab;
  string coma=",";;
    if ((size(lab)==1) and lab[1]==-1)
    {slab=""; coma="";}
  else
  {
    slab=string(lab[1]);
    if (size(lab)>=1)
    {
      for (i=2;i<=size(lab);i++){slab=string(slab,",",lab[i]);}
    }
  }
  Li1=string(T,"[[",slab+coma,"0]]:=",L[8],"=0; ");
  Li2=string(T,"[[",slab+coma,"1]]:=",L[8],"<>0; ");
  return(list(Li1,Li2));
}

//*****************End of BuildtreetoMaple*********************

//*****************Begin of Selectcases************************

// given an intvec with sum=n
// it returns the list of intvect with the sum=n+1
static proc comp1(intvec l)
{
  list L;
  int p=size(l);
  int i;
  if (p==0){return(l);}
  if (p==1){return(list(intvec(l[1]+1)));}
  L[1]=intvec((l[1]+1),l[2..p]);
  L[p]=intvec(l[1..p-1],(l[p]+1));
  for (i=2;i<p;i++)
  {
    L[i]=intvec(l[1..(i-1)],(l[i]+1),l[(i+1)..p]);
  }
  return(L);
}

// comp: p-compositions of n
// input
//   int n;
//   int p;
// return
//   the list of all intvec (p-composition of n)
static proc comp(int n,int p)
{
  if (n<0){ERROR("comp was called with negative argument");}
  if (n==0){return(list(0:p));}
  int i;
  int k;
  list L1=comp(n-1,p);
  list L=comp1(L1[1]);
  list l;
  list la;
  for (i=2; i<=size(L1);i++)
  {
    l=comp1(L1[i]);
    for (k=1;k<=size(l);k++)
    {
      if(not(memberpos(l[k],L)[1]))
      {L[size(L)+1]=l[k];}
    }
  }
  return(L);
}

// given the matrices of coefficients and monomials m amd m1 of
// two polynomials (the first one contains all the terms of f
// and the second only those of f
// it returns the list with the comon monomials and the list of coefficients
// of the polynomial f with zeroes if necessary.
static proc adaptcoef(matrix m, matrix m1)
{
  int i;
  int j;
  int ncm=ncols(m);
  int ncm1=ncols(m1);
  ideal T;
  for (i=1;i<=ncm;i++){T[i]=m[1,i];}
  ideal C;
  for (i=1;i<=ncm;i++){C[i]=0;}
  for (i=1;i<=ncm;i++)
  {
    j=1;
    while((j<ncm1) and (m1[1,j]>m[1,i])){j++;}
    if (m1[1,j]==m[1,i]){C[i]=m1[2,j];}
  }
  return(list(T,C));
}

// given teh ideal of non-null conditions and an intvec lambda
// with the exponents of each w in W
// it returns the polynomial prod (w_i)^(lambda_i).
static proc WW(ideal W, intvec lambda)
{
  if (size(W)==0){return(poly(1));}
  poly w=1;
  int i;
  for (i=1;i<=ncols(W);i++)
  {
    w=w*(W[i])^(lambda[i]);
  }
  return(w);
}

// given a polynomial f and the non-null conditions W
// WPred eliminates the factors in f that are in W
// ring @PAB
// input:
//   poly f:
//   ideal W  of non-null conditions (already supposed that it is facvar)
// output:
//   poly f2  where the non-null conditions in W have been dropped from f
static proc WPred(poly f, ideal W)
{
  if (f==0){return(f);}
  def l=factorize(f,2);
  int i;
  poly f1=1;
  for(i=1;i<=size(l[1]);i++)
  {
    if (memberpos(l[1][i],W)[1]){;}
    else{f1=f1*((l[1][i])^(l[2][i]));}
  }
  return(f1);
}

//genimage
// ring @R
//input:
//   poly f1, idel N1,ideal W1,poly f2, ideal N2, ideal W2
//   corresponding to two polynomials having the same lpp
//   f1 in the redspec given by N1,W1,  f2 in the redspec given by N2,W2
//output:
//   the list of (ideal GG, list(list r1, list r2))
//   where g an ideal whose elements have the same lpp as f1 and f2
//   that specialize well to f1 in N1,W1 and to f2 in N2,W2.
//   If it doesn't exist a genimage, then g=ideal(0).
static proc genimage(poly f1, ideal N1, ideal W1, poly f2, ideal N2, ideal W2)
{
  int i; ideal W12;  poly ff1; poly g1=0; ideal GG;
  int tt=1;
  // detect weather f1 reduces to 0 on segment 2
  ff1=pnormalform(f1,N2,W2);
  if (ff1==0)
  {
    // detect weather N1 is included in N2
    def RR=basering;
    setring @P;
    def NP1=imap(RR,N1);
    def NP2=imap(RR,N2);
    attrib(NP2,"isSB",1);
    poly nr;
    i=1;
    while ((tt) and (i<=size(NP1)))
    {
      nr=reduce(NP1[i],NP2);
      if (nr!=0){tt=0;}
      i++;
    }
    setring(RR);
  }
  else{tt=0;}
  if (tt==1)
  {
    // detect weather W1 intersect W2 is non-empty
    for (i=1;i<=size(W1);i++)
    {
      if (memberpos(W1[i],W2)[1])
      {
        W12[size(W12)+1]=W1[i];
      }
      else
      {
        if (nonnull(W1[i],N2,W2))
        {
          W12[size(W12)+1]=W1[i];
        }
      }
    }
    for (i=1;i<=size(W2);i++)
    {
      if (not(memberpos(W2[i],W12)[1]))
      {
        W12[size(W12)+1]=W2[i];
      }
    }
  }
  if (tt==1){g1=extendpoly(f1,N1,W12);}
  if (g1!=0)
  {
    if (pnormalform(g1,N1,W1)==0)
    {
      GG=f1,g1;
    }
    else
    {
      GG=g1;
    }
    return(GG);
  }

  // begins the second step;
  int bound=6;
  // in ring @R
  int j; int g=0; int alpha; int r1; int s1=1; int s2=1;
  poly G;
  matrix qT;
  matrix T;
  ideal N10;
  poly GT;
  ideal N12=N1,N2;
  def varx=maxideal(1);
  int nx=size(varx);
  poly pvarx=1;
  for (i=1;i<=nx;i++){pvarx=pvarx*varx[i];}
  def m=coef(43*f1+157*f2,pvarx);
  def m1=coef(f1,pvarx);
  def m2=coef(f2,pvarx);
  list L1=adaptcoef(m,m1);
  list L2=adaptcoef(m,m2);
  ideal Tm=L1[1];
  ideal c1=L1[2];
  ideal c2=L2[2];
  poly ww1;
  poly ww2;
  poly cA1;
  poly cB1;
  matrix TT;
  poly H;
  list r;
  ideal q;
  poly mu;
  ideal N;

  // in ring @PAB
  list Px=ringlist(@P);
  list v="@A","@B";
  Px[2]=Px[2]+v;
  def npx=size(Px[3][1][2]);
  Px[3][1][2]=1:(npx+size(v));
  def @PAB=ring(Px);
  setring(@PAB);

  poly PH;
  ideal NP;
  list rP;
  def PN1=imap(@R,N1);
  def PW1=imap(@R,W1);
  def PN2=imap(@R,N2);
  def PW2=imap(@R,W2);
  def a1=imap(@R,c1);
  def a2=imap(@R,c2);
  matrix PT;
  ideal PN;
  ideal PN12=PN1,PN2;
  PN=liftstd(PN12,PT);
  list compos1;
  list compos2;
  list compos0;
  intvec comp0;
  poly w1=0;
  poly w2=0;
  poly h;
  poly cA=0;
  poly cB=0;
  int t=0;
  list l;
  poly h1;
  g=0;
  while ((g<=bound) and not(t))
  {
    compos0=comp(g,2);
    r1=1;
    while ((r1<=size(compos0)) and not(t))
    {
      comp0=compos0[r1];
      if (comp0[1]<=bound/2)
      {
        compos1=comp(comp0[1],ncols(PW1));
        s1=1;
        while ((s1<=size(compos1)) and not(t))
        {
          if (comp0[2]<=bound/2)
          {
            compos2=comp(comp0[2],ncols(PW2));
            s2=1;
            while ((s2<=size(compos2)) and not(t))
            {
              w1=WW(PW1,compos1[s1]);
              w2=WW(PW2,compos2[s2]);
              h=@A*w1*a1[1]-@B*w2*a2[1];
              h=reduce(h,PN);
              if (h==0){cA=1;cB=-1;}
              else
              {
                l=factorize(h,2);
                h1=1;
                for(i=1;i<=size(l[1]);i++)
                {
                  if ((memberpos(@A,variables(l[1][i]))[1]) or  (memberpos(@B,variables(l[1][i]))[1]))
                  {h1=h1*l[1][i];}
                }
                cA=diff(h1,@B);
                cB=diff(h1,@A);
              }
              if ((cA!=0) and (cB!=0) and (jet(cA,0)==cA) and (jet(cB,0)==cB))
              {
                t=1;
                alpha=1;
                while((t) and (alpha<=ncols(a1)))
                {
                  h=cA*w1*a1[alpha]+cB*w2*a2[alpha];
                  if (not(reduce(h,PN,1)==0)){t=0;}
                  alpha++;
                }
              }
              else{t=0;}
              s2++;
            }
          }
          s1++;
        }
      }
      r1++;
    }
    g++;
  }
  setring(@R);
  ww1=imap(@PAB,w1);
  ww2=imap(@PAB,w2);
  T=imap(@PAB,PT);
  N=imap(@PAB,PN);
  cA1=imap(@PAB,cA);
  cB1=imap(@PAB,cB);
  if (t)
  {
    G=0;
    for (alpha=1;alpha<=ncols(Tm);alpha++)
    {
      H=cA1*ww1*c1[alpha]+cB1*ww2*c2[alpha];
      setring(@PAB);
      PH=imap(@R,H);
      PN=imap(@R,N);
      rP=division(PH,PN);
      setring(@R);
      r=imap(@PAB,rP);
      if (r[2][1]!=0){ERROR("the division is not null and it should be");}
      q=r[1];
      qT=transpose(matrix(q));
      N10=N12;
      for (i=size(N1)+1;i<=size(N1)+size(N2);i++){N10[i]=0;}
      G=G+(cA1*ww1*c1[alpha]-(matrix(N10)*T*qT)[1,1])*Tm[alpha];
    }
    GG=ideal(G);
  }
  else{GG=ideal(0);}
  return(GG);
}

// purpose: given a polynomial f (in the reduced basis)
//          the null-conditions ideal N in the segment
//          end the set of non-null polynomials common to the segment and
//          a new segment,
//          to obtain an equivalent polynomial with a leading coefficient
//          that is non-null in the second segment.
// input:
// poly f:    a polynomials of the reduced basis in the segment (N,W)
// ideal N:   the null-conditions ideal in the segment
// ideal W12: the set of non-null polynomials common to the segment and
//            a second segment
static proc extendpoly(poly f, ideal N, ideal W12)
{
  int bound=4;
  ideal cfs;
  ideal cfsn;
  ideal ppfs;
  poly p=f;
  poly fn;
  poly lm; poly lc;
  int tt=0;
  int i;
  while (p!=0)
  {
    lm=leadmonom(p);
    lc=leadcoef(p);
    cfs[size(cfs)+1]=lc;
    ppfs[size(ppfs)+1]=lm;
    p=p-lc*lm;
  }
  def lcf=cfs[1];
  int r1=0; int s1;
  def RR=basering;
  setring @P;
  list compos1;
  poly w1;
  ideal q;
  def lcfp=imap(RR,lcf);
  def W=imap(RR,W12);
  def Np=imap(RR,N);
  def cfsp=imap(RR,cfs);
  ideal cfspn;
  matrix T;
  ideal H=lcfp,Np;
  def G=liftstd(H,T);
  list r;
  while ((r1<=bound) and not(tt))
  {
    compos1=comp(r1,ncols(W));
    s1=1;
    while ((s1<=size(compos1)) and not(tt))
    {
      w1=WW(W,compos1[s1]);
      cfspn=ideal(0);
      cfspn[1]=w1;
      tt=1;
      i=2;
      while ((i<=size(cfsp)) and (tt))
      {
        r=division(w1*cfsp[i],G);
        if (r[2][1]!=0){tt=0;}
        else
        {
          q=r[1];
          cfspn[i]=(T*transpose(matrix(q)))[1,1];
        }
        i++;
      }
      s1++;
    }
    r1++;
  }
  setring RR;
  if (tt)
  {
    cfsn=imap(@P,cfspn);
    fn=0;
    for (i=1;i<=size(ppfs);i++)
    {
      fn=fn+cfsn[i]*ppfs[i];
    }
  }
  else{fn=0;}
  return(fn);
}

// nonnull
// ring @P (or @R)
// input:
//   poly f
//   ideal N
//   ideal W
// output:
//   1 if f is nonnull in the segment (N,W)
//   0 if it can be zero
static proc nonnull(poly f, ideal N, ideal W)
{
  int tt;
  ideal N0=N;
  N0[size(N0)+1]=f;
  poly h=1;
  int i;
  for (i=1;i<=size(W);i++){h=h*W[i];}
  def RR=basering;
  setring(@P);
  list Px=ringlist(@P);
  list v="@C";
  Px[2]=Px[2]+v;
  def npx=size(Px[3][1][2]);
  Px[3][1][1]="dp";
  Px[3][1][2]=1:(npx+size(v));
  def @PC=ring(Px);
  setring(@PC);
  def N1=imap(RR,N0);
  def h1=imap(RR,h);
  ideal G=1-@C*h1;
  G=G+N1;
  option(redSB);
  ideal G1=std(G);
  if (G1[1]==1){tt=1;} else{tt=0;}
  setring(RR);
  return(tt);
}

// decide
// input:
//   given two corresponding polynomials g1 and g2 with the same lpp
//   g1 belonging to the basis in the segment N1,W1
//   g2 belonging to the basis in the segment N2,W2
// output:
//   an ideal (with a single polynomial or more if a sheaf is needed)
//   that specializes well on both segments to g1 and g2 respectivelly.
//   If ideal(0) is output, then no such polynomial nor sheaf exists.
static proc decide(poly g1, ideal N1, ideal W1, poly g2, ideal N2, ideal W2)
{
  poly S;
  poly S1;
  poly S2;
  S=leadcoef(g2)*g1-leadcoef(g1)*g2;
  def RR=basering;
  setring(@RP);
  def SR=imap(RR,S);
  def N1R=imap(RR,N1);
  def N2R=imap(RR,N2);
  attrib(N1R,"isSB",1);
  attrib(N2R,"isSB",1);
  poly S1R=reduce(SR,N1R);
  poly S2R=reduce(SR,N2R);
  setring(RR);
  S1=imap(@RP,S1R);
  S2=imap(@RP,S2R);
  if ((S2==0) and (nonnull(leadcoef(g1),N2,W2))){return(ideal(g1));}
  if ((S1==0) and (nonnull(leadcoef(g2),N1,W1))){return(ideal(g2));}
  if ((S1==0) and (S2==0))
  {
    return(ideal(g1,g2));
  }
  return(ideal(genimage(g1,N1,W1,g2,N2,W2)));
}

// input:  the tree (list) from buildtree output
// output: the list of terminal vertices.
static proc finalcases(list T)
//"USAGE:   finalcases(T);
//          T is the list provided by buildtree
//RETURN:   A list with the CGS determined by buildtree.
//          Each element of the list represents one segment
//          of the terminal vertices of buildtree givieng the CGS.
//          The list elements have the following structure:
//           [1]: label (an intvec(1,0,..)) that indicates the position
//                in the buildtree but that is irrelevant for the CGS
//           [2]: 1 (integer) it is also irrelevant and indicates
//                that this was a terminal vertex in buildtree.
//           [3]: the reduced basis of the segment.
//           [4], [5], [6]: the red-representation of the segment
//                [4] are the null-conditions radical ideal N,
//                [5] are the non-null polynomials set (ideal) W,
//                [6] is the set of prime components (ideals) of N.
//           [7]: is the set of lpp
//           [8]: poly 1 (irrelevant) is the condition to branch (but no
//                more branch is necessary in the discussion, so 1 is the result.
//NOTE:     It can be called having as argument the list output by buildtree
//KEYWORDS: buildtree, buildtreetoMaple, CGS
//EXAMPLE:  finalcases; shows an example"
{
  int i;
  list L;
  for (i=1;i<=size(T);i++)
  {
    if (T[i][2])
    {L[size(L)+1]=T[i];}
  }
  return(L);
}
//example
//{ "EXAMPLE:"; echo = 2;
//  ring R=(0,a1,a2,a3,a4),(x1,x2,x3,x4),dp;
//  ideal F=x4-a4+a2, x1+x2+x3+x4-a1-a3-a4, x1*x3*x4-a1*a3*a4, x1*x3+x1*x4+x2*x3+x3*x4-a1*a4-a1*a3-a3*a4;
//  def T=buildtree(F);
//  setglobalrings();
//  finalcases(T);
//}

// input:  the list of terminal vertices of buildtree (output of finalcases)
// output: the same terminal vertices grouped by lpp
static proc groupsegments(list T)
{
  int i;
  list L;
  list lpp;
  list lp;
  list ls;
  int n=size(T);
  lpp[1]=T[n][7];
  L[1]=list(lpp[1],list(list(T[n][1],T[n][3],T[n][4],T[n][5],T[n][6])));
  if (n>1)
  {
    for (i=1;i<=size(T)-1;i++)
    {
      lp=memberpos(T[n-i][7],lpp);
      if(lp[1]==1)
      {
        ls=L[lp[2]][2];
        ls[size(ls)+1]=list(T[n-i][1],T[n-i][3],T[n-i][4],T[n-i][5],T[n-i][6]);
        L[lp[2]][2]=ls;
      }
      else
      {
        lpp[size(lpp)+1]=T[n-i][7];
        L[size(L)+1]=list(T[n-i][7],list(list(T[n-i][1],T[n-i][3],T[n-i][4],T[n-i][5],T[n-i][6])));
      }
    }
  }
  //"L in groupsegments="; L;
  return(L);
}

// eliminates repeated elements form an ideal
static proc elimrepeated(ideal F)
{
  int i;
  int j;
  ideal FF;
  FF[1]=F[1];
  for (i=2;i<=ncols(F);i++;)
  {
    if (not(memberpos(F[i],FF)[1]))
    {
      FF[size(FF)+1]=F[i];
    }
  }
  return(FF);
}

// decide F is the same as decide but allows as first element a sheaf F
static proc decideF(ideal F,ideal N,ideal W, poly f2, ideal N2, ideal W2)
{
  int i;
  ideal G=F;
  ideal g;
  if (ncols(F)==1) {return(decide(F[1],N,W,f2,N2,W2));}
  for (i=1;i<=ncols(F);i++)
  {
    G=G+decide(F[i],N,W,f2,N2,W2);
  }
  return(elimrepeated(G));
}

// newredspec
// input:  two redspec in the form of N,W and Nj,Wj
// output: a redspec representing the minimal redspec segment that contains
//         both input segments.
static proc newredspec(ideal N,ideal W, ideal Nj, ideal Wj)
{
  ideal nN;
  ideal nW;
  int u;
  def RR=basering;
  setring(@P);
  list r;
  def Np=imap(RR,N);
  def Wp=imap(RR,W);
  def Njp=imap(RR,Nj);
  def Wjp=imap(RR,Wj);
  Np=intersect(Np,Njp);
  ideal WR;
  for(u=1;u<=size(Wjp);u++)
  {
    if(nonnull(Wjp[u],Np,Wp)){WR[size(WR)+1]=Wjp[u];}
  }
  for(u=1;u<=size(Wp);u++)
  {
    if((not(memberpos(Wp[u],WR)[1])) and (nonnull(Wp[u],Njp,Wjp)))
    {
      WR[size(WR)+1]=Wp[u];
    }
  }
  r=redspec(Np,WR);
  option(redSB);
  Np=std(r[1]);
  Wp=r[2];
  setring(RR);
  nN=imap(@P,Np);
  nW=imap(@P,Wp);
  return(list(nN,nW));
}

// selectcases
// input:
//   list bT: the list output by buildtree.
// output:
//   list L   it contins the list of segments allowing a common
//            reduced basis. The elements of L are of the form
//            list (lpp,B,list(list(N,W,L),..list(N,W,L)) )
static proc selectcases(list bT)
{
  list T=groupsegments(finalcases(bT));
  //NEW
  //groupredtocan(T);
  list T0=bT[1];
             // first element of the list of buildtree
  list TT0;
  TT0[1]=list(T0[7],T0[3],list(list(T0[4],T0[5],T0[6])));
             // first element of the output of selectcases
  list T1=T; // the initial list; it is only actualized (split)
             // when a segment is completly revised (all split are
             // already be considered);
             // ( (lpp, ((lab,B,N,W,L),.. ()) ), .. (..) )
  list TT;   // the output list ( (lpp,B,((N,W,L),..()) ),.. (..) )
  // case i
  list S1;   // the segments in case i T1[i][2]; ( (lab,B,N,W,L),..() )
  list S2;   // the segments in case i that are being summarized in
             // actual segment ( (N,W,L),..() )
  list S3;   // the segments in case i that cannot be summarized in
             // the actual case. When the case is finished a new case
             // is created with them ( (lab,B,N,W,L),..() )
  list s3;   // list of integers s whose segment cannot be summarized
             // in the actual case
  ideal lpp; // the summarized lpp (can contain repetitions)
  ideal lppi;// in process of sumarizing lpp (can contain repetitions)
  ideal B;   // the summarized B (can contain polynomials with
             // the same lpp (sheaves))
  ideal Bi;  // in process of summarizing B (can contain polynomials with
             // the same lpp (sheaves))
  ideal N;   // the summarized N
  ideal W;   // the summarized W
  ideal F;   // the summarized poly j (can contain a sheaf instead of
             // a single poly)
  ideal FF;  // the same as F but it can be ideal(0)
  poly lpj;
  poly fj;
  ideal Nj;
  ideal Wj;
  ideal G;
  int i;     // the index of the case i in T1;
  int j;     // the index of the polynomial j of the basis
  int s;     // the index of the segment s in S1;
  int u;
  int tests; // true if al the polynomial in segment s have been generalized;
  list r;
  // initializing the new list
  i=1;
  while(i<=size(T1))
  {
    S1=T1[i][2]; // ((lab,B,N,W,L)..) of the segments in case i
    if (size(S1)==1)
    {
      TT[i]=list(T1[i][1],S1[1][2],list(list(S1[1][3],S1[1][4],S1[1][5])));
    }
    else
    {
      S2=list();
      S3=list(); // ((lab,B,N,W,L)..) of the segments in case i to
                 // create another segment i+1
      s3=list();
      B=S1[1][2];
      Bi=ideal(0);
      lpp=T1[i][1];
      j=1;
      tests=1;
      while (j<=size(S1[1][2]))
      { // j desings the new j-th polynomial
        N=S1[1][3];
        W=S1[1][4];
        F=ideal(S1[1][2][j]);
        s=2;
        while (s<=size(S1) and not(memberpos(s,s3)[1]))
        { // s desings the new segment s
          fj=S1[s][2][j];
          Nj=S1[s][3];
          Wj=S1[s][4];
          FF=decideF(F,N,W,fj,Nj,Wj);
          if (FF[1]==0)
          {
            if (@ish)
            {
              "Warning: Dealing with an homogeneous ideal";
              "mrcgs was not able to summarize all lpp cases into a single segment";
              "Please send a mail with your Problem to antonio.montes@upc.edu";
              "You found a counterexample of the complete success of the actual mrcgs algorithm";
              //NEW
              "f1:"; F; "N1:"; N; "W1:"; W; "f2:"; fj; "N2:"; Nj; "W2:"; Wj;
            }
            S3[size(S3)+1]=S1[s];
            s3[size(s3)+1]=s;
            tests=0;
          }
          else
          {
            F=FF;
            lpj=leadmonom(fj);
            r=newredspec(N,W,Nj,Wj);
            N=r[1];
            W=r[2];
          }
          s++;
        }
        if (Bi[1]==0){Bi=FF;}
        else
        {
          Bi=Bi+FF;
        }
        j++;
      }
      if (tests)
      {
        B=Bi;
        lpp=ideal(0);
        for (u=1;u<=size(B);u++){lpp[u]=leadmonom(B[u]);}
      }
      for (s=1;s<=size(T1[i][2]);s++)
      {
        if (not(memberpos(s,s3)[1]))
        {
          S2[size(S2)+1]=list(S1[s][3],S1[s][4],S1[s][5]);
        }
      }
      TT[i]=list(lpp,B,S2);
      // for (s=1;s<=size(s3);s++){S1=delete(S1,s);}
      T1[i][2]=S2;
      if (size(S3)>0){T1=insert(T1,list(T1[i][1],S3),i);}
    }
    i++;
  }
  for (i=1;i<=size(TT);i++){TT0[i+1]=TT[i];}
  return(TT0);
}

//*****************End of Selectcases**************************

//*****************Begin of CanTree****************************

// equalideals
// input: 2 ideals F and G;
// output: 1 if they are identical (the same polynomials in the same order)
//         0 else
static proc equalideals(ideal F, ideal G)
{
  int i=1; int t=1;
  if (size(F)!=size(G)){return(0);}
  while ((i<=size(F)) and (t))
  {
    if (F[i]!=G[i]){t=0;}
    i++;
  }
  return(t);
}

// delintvec
// input: intvec V
//        int i
// output:
//        intvec W (equal to V but the coordinate i is deleted
static proc delintvec(intvec V, int i)
{
  int j;
  intvec W;
  for (j=1;j<i;j++){W[j]=V[j];}
  for (j=i+1;j<=size(V);j++){W[j-1]=V[j];}
  return(W);
}

// redtocanspec
// Computes the canonical representation of a redspec (N,W,L).
// input:
//    ideal N (null conditions, must be radical)
//    ideal W (non-null conditions ideal)
//    list L  must contain the radical decomposition of N.
// output:
//    the list of elements of the (ideal N1,list(ideal M11,..,ideal M1k))
//    determining the canonical representation of the difference of
//    V(N) \ V(h), where h=prod(w in W).
static proc redtocanspec(intvec lab, int child, list rs)
{
  ideal N=rs[1]; ideal W=rs[2]; list L=rs[3];
  intvec labi; intvec labij;
  int childi;
  int i; int j; list L0;
  L0[1]=list(lab,size(L));
  if (W[1]==0)
  {
    for (i=1;i<=size(L);i++)
    {
      labi=lab,child+i;
      L0[size(L0)+1]=list(labi,1,L[i]);
      labij=labi,1;
      L0[size(L0)+1]=list(labij,0,ideal(1));
    }
    return(L0);
  }
  if (N[1]==1)
  {
    L0[1]=list(lab,1);
    labi=lab,child+1;
    L0[size(L0)+1]=list(labi,1,ideal(1));
    labij=labi,1;
    L0[size(L0)+1]=list(labij,0,ideal(1));
  }
  def RR=basering;
  setring(@P);
  ideal Np=imap(RR,N);
  ideal Wp=imap(RR,W);
  poly h=1;
  for (i=1;i<=size(Wp);i++){h=h*Wp[i];}
  list Lp=imap(RR,L);
  list r; list Ti; list LL;
  LL[1]=list(lab,size(Lp));
  for (i=1;i<=size(Lp);i++)
  {
    Ti=minGTZ(Lp[i]+h);
    for(j=1;j<=size(Ti);j++)
    {
      option(redSB);
      Ti[j]=std(Ti[j]);
    }
    labi=lab,child+i;
    childi=size(Ti);
    LL[size(LL)+1]=list(labi,childi,Lp[i]);
    for (j=1;j<=childi;j++)
    {
      labij=labi,j;
      LL[size(LL)+1]=list(labij,0,Ti[j]);
    }
  }
  LL[1]=list(lab,size(Lp));
  setring(RR);
  return(imap(@P,LL));
}

// difftocanspec
// Computes the canonical representation of a diffspec V(N) \ V(M)
// input:
//    intvec lab: label where to hang the canspec
//    list  N ideal of null conditions.
//    ideal M ideal of the variety to be substacted
// output:
//    the list of elements determining the canonical representation of
//    the difference  V(N) \ V(M):
//      ( (intvec(i),children), ...(lab, children, prime ideal),...)
static proc difftocanspec(intvec lab, int child, ideal N, ideal M)
{
  int i; int j; list LLL;
  def RR=basering;
  setring(@P);
  ideal Np=imap(RR,N);
  ideal Mp=imap(RR,M);
  def L=minGTZ(Np);
  for(j=1;j<=size(L);j++)
  {
    option(redSB);
    L[j]=std(L[j]);
  }
  intvec labi; intvec labij;
  int childi;
  list LL;
  if ((Mp[1]==0) or ((size(L)==1) and (L[1][1]==1)))
  {
    //LL[1]=list(lab,1);
    //labi=lab,1;
    //LL[2]=list(labi,1,ideal(1));
    //labij=labi,1;
    //LL[3]=list(labij,0,ideal(1));
    setring(RR);
    return(LLL);
  }
  list r; list Ti;
  def k=0;
  LL[1]=list(lab,0);
  for (i=1;i<=size(L);i++)
  {
    Ti=minGTZ(L[i]+Mp);
    for(j=1;j<=size(Ti);j++)
    {
      option(redSB);
      Ti[j]=std(Ti[j]);
    }
    if (not((size(Ti)==1) and (equalideals(L[i],Ti[1]))))
    {
      k++;
      labi=lab,child+k;
      childi=size(Ti);
      LL[size(LL)+1]=list(labi,childi,L[i]);
      for (j=1;j<=childi;j++)
      {
        labij=labi,j;
        LL[size(LL)+1]=list(labij,0,Ti[j]);
      }
    }
    else{setring(RR); return(LLL);}
  }
  if (size(LL)>0)
  {
    LL[1]=list(lab,k);
    setring(RR);
    return(imap(@P,LL));
  }
  else {setring(RR); return(LLL);}
}

// tree
// purpose: given a label and the list L of vertices of the tree,
//          whose content
//          are of the form list(intvec lab, int children, ideal P)
//          to obtain the vertex and its position
// input:
//  intvec lab: label of the vertex
//  list:  L    the list containing the vertices
// output:
//  list   V    the vertex list(lab, children, P)
static proc tree(intvec lab,list L)
{
  int i=0; int tt=1; list V; intvec labi;
  while ((i<size(L)) and (tt))
  {
    i++;
    labi=L[i][1];
    if (labi==lab)
    {
      V=list(L[i],i);
      tt=0;
    }
  }
  if (tt==0){return(V);}
  else{return(list(list(intvec(0)),0));}
}

// GCR (generalized canonical representation)
// new structure of a GCR

// L is a list of vertices V of the GCR.
// first vertex=list(intvec lab, int children, ideal lpp, ideal B)
// other vertices=list(intvec lab, int children, ideal P)
// the individual vertices can be accessed with the function tree
// by the call  V=tree(lab,L), that outputs the vertex if it exists
// and its position in L, or nothing if it does not exist.
// The first element of the list must be the root of the tree and has
// label lab=i, and other information.

// example:
// the canonical representation
// V(a^2-ac-ba+c-abc) \ (union( V(b,a), V(c,a), V(b,a-c), V(c,a-b)))
// is represented by  the list
// L=((intvec(i),children=1,lpp,B),(intvec(i,1),4,ideal(a^2-ac-ba+c-abc)),
//    (intvec(i,1,1),0,ideal(b,a)),     (intvec(i,1,2),0,ideal(c,a)),
//    (intvec(i,1,3),0,ideal(b,a-c)),   (intvec(i,1,4),0,ideal(c,a-b))
//   )
// example:
// the canonical representation
// (V(a)\(union(V(c,a),V(b+c,a),V(b,a)))) union
// (V(b)\(union(V(b,a),V(b,a-c))))        union
// (V(c)\(union(V(c,a),V(c,a-b))))
// is represented by  the list
// L=((i,children=3,lpp,B),
//    (intvec(i,1),3,ideal(a)),
//    (intvec(i,1,1),0,(c,a)),(intvec(i,1,2),0,(b+c,a)),(intvec(i,1,3),0,(b,a)),
//    (intvec(i,2),2,ideal(b)),
//    (intvec(i,2,1),0,(b,a)),(intvec(i,2,2),0,(b,a-c)),
//    (intvec(i,3),2,ideal(c)),
//    (intvec(i,3,1),0,(c,a)),(intvec(i,3,2),0,(c,a-b))
//   )
// If L is the list in the last example, the call
// tree(intvec(i,2,1),L) will output   ((intvec(i,2,1),0,(b,a)),7)

// GCR
// input: list T is supposed to be an element L[i] of selectcases:
//        T= list( ideal lpp, ideal B, list(N,W,L),.., list(N,W,L))
// output: the list L of vertices being the GCR of the addition of
//         all the segments in T.
//         list(list(intvec lab, int children, ideal lpp, ideal B),
//              list(intvec lab, int children, ideal P),..
//         )
static proc GCR(intvec lab, list case)
{
  int i; int ii; int t;
  list @L;
  @L[1]=list(lab,0,case[1],case[2]);
  exportto(Top,@L);
  int j;
  list u; intvec labu; int childu;
  list v; intvec labv; int childv;
  list T=case[3];
  for (j=1;j<=size(T);j++)
  {
    t=addcase(lab,T[j]);
    deletebrotherscontaining(lab);
  }
  relabelingindices(lab,lab);
  list L=@L;
  kill @L;
  return(L);
}

// sorbylab:
// pupose: given the list of mrcgs to order is by increasing label
static proc sortbylab(list L)
{
  int n=L[1][2];
  int i; int j;
  list H=L;
  list LL;
  list L1;
  //LL[1]=L[1];
  //H=delete(H,1);
  while (size(H)!=0)
  {
    j=1;
    L1=H[1];
    for (i=1;i<=size(H);i++)
    {
      if(lesslab(H[i],L1)){j=i;L1=H[j];}
    }
    LL[size(LL)+1]=L1;
    H=delete(H,j);
  }
  return(LL);
}

// lesslab
// purpose: given two elements of the list of mrcgs it
// returns 1 if the label of the first is less than that of the second
static proc lesslab(list l1, list l2)
{
  intvec lab1=l1[1];
  intvec lab2=l2[1];
  int n1=size(lab1);
  int n2=size(lab2);
  int n=n1;
  if (n2<n1){n=n2;}
  int tt=0;
  int j=1;
  while ((lab1[j]==lab2[j]) and (j<n)){j++;}
  if (lab1[j]<lab2[j]){tt=1;}
  if ((j==n) and (lab1[j]==lab2[j]) and (n2>n1)){tt=1;}
  return(tt);
}

// cantree
// input:  the list provided by selectcases
// output: the list providing the canonicaltree
static proc cantree(list S)
{
  string method=" ";
  list T0=S[1];
    // first element of the list of selectcases
  int i; int j;
  list L;
  list T;
  L[1]=list(intvec(0),size(S)-1,T0[1],T0[2],T0[3][1],method);
  for (i=2;i<=size(S);i++)
  {
    T=GCR(intvec(i-1),S[i]);
    T=sortbylab(T);
    for (j=1;j<=size(T);j++)
    {L[size(L)+1]=T[j];}
  }
  return(L);
}

// addcase
// recursive routine that adds to the list @L, (an alredy GCR)
// a new redspec rs=(N,W,L);
// and returns the test t whose value is
// 0 if the new canspec is not to be hung to the fathers vertex,
// 1 if yes.
static proc addcase(intvec labu, list rs)
{
  int i; int j; int childu; ideal Pu;
  list T; int nchildu;
  def N=rs[1]; def W=rs[2]; def PN=rs[3];
  ideal NN; ideal MM;
  int tt=1;
  poly h=1; for (i=1;i<=size(W);i++){h=h*W[i];}
  list u=tree(labu,@L); childu=u[1][2];
  list v; intvec labv; int childv; list w; intvec labw;
  if (childu>0)
  {
    v=firstchild(u[1][1]);
    while(v[2][1]!=0)
    {
      labv=v[1][1];
      w=firstchild(labv);
      while(w[2][1]!=0)
      {
        labw=w[1][1];
        if(addcase(labw,rs)==0)
        {tt=0;}
        w=nextbrother(labw);
      }
      u=tree(labu,@L);
      childu=u[1][2];
      v=nextbrother(v[1][1]);
    }
    deletebrotherscontaining(labu);
    relabelingindices(labu,labu);
  }
  if (tt==1)
  {
    u=tree(labu,@L);
    nchildu=lastchildrenindex(labu);
    if (size(labu)==1)
    {
      T=redtocanspec(labu,nchildu,rs);
      tt=0;
    }
    else
    {
      NN=N;
      if (containedP(u[1][3],N)){tt=0;}
      for (i=1;i<=size(u[1][3]);i++)
      {
        NN[size(NN)+1]=u[1][3][i];
      }
      MM=NN;
      MM[size(MM)+1]=h;
      T=difftocanspec(labu,nchildu,NN,MM);
    }
    if (size(T)>0)
    {
      @L[u[2]][2]=@L[u[2]][2]+T[1][2];
      for (i=2;i<=size(T);i++){@L[size(@L)+1]=T[i];}
      if (size(labu)>1)
      {
        simplifynewadded(labu);
      }
    }
    else{tt=1;}
  }
  return(tt);
}

// reduceR
// reduces the polynomial f wrt N, in the ring @P
static proc reduceR(poly f, ideal N)
{
  def RR=basering;
  setring(@P);
  poly fP=imap(RR,f);
  ideal NP=imap(RR,N);
  attrib(NP,"isSB",1);
  poly rp=reduce(fP,NP);
  setring(RR);
  return(imap(@P,rp));
}

// containedP
// returns 1 if ideal Pu is contained in ideal Pv
// returns 0 if not
// in ring @P
static proc containedP(ideal Pu,ideal Pv)
{
  int t=1;
  int n=ncols(Pu);
  int i=0;
  poly r=0;
  while ((t) and (i<n))
  {
    i++;
    r=reduceR(Pu[i],Pv);
    if (r!=0){t=0;}
  }
  return(t);
}

// simplifynewadded
// auxiliary routine of addcase
// when a new redspec is added to a non terminal vertex,
// it is applied to simplify the addition.
// When Pu==Pv, the children of w are hung from u fathers
// and deleted the whole new addition.
// Finally, deletebrotherscontaining is applied to u fathers
// in order to eliminate branches contained.
static proc simplifynewadded(intvec labu)
{
  int t; int ii; int k; int kk; int j;
  intvec labfu=delintvec(labu,size(labu)); list fu; int childfu;
  list u=tree(labu,@L); int childu=u[1][2]; ideal Pu=u[1][3];
  list v; intvec labv; int childv; ideal Pv;
  list w; intvec labw; intvec nlab; list ww;
  if (childu>0)
  {
    v=firstchild(u[1][1]); labv=v[1][1]; childv=v[1][2]; Pv=v[1][3];
    ii=0;
    t=0;
    while ((not(t)) and (ii<childu))
    {
      ii++;
      if (equalideals(Pu,Pv))
      {
        fu=tree(labfu,@L);
        childfu=fu[1][2];
        j=lastchildrenindex(fu[1][1])+1;
        k=0;
        w=firstchild(v[1][1]);
        childv=v[1][2];
        for (kk=1;kk<=childv;kk++)
        {
          if (kk<childv){ww=nextbrother(w[1][1]);}
          nlab=labfu,j;
          @L[w[2]][1]=nlab;
          j++;
          if (kk<childv){w=ww;}
        }
        childfu=fu[1][2]+childv-1;
        @L[fu[2]][2]=childfu;
        @L[v[2]][2]=0;
        t=1;
        deleteverts(labu);
      }
    }
  }
  deletebrotherscontaining(labfu);
}

// given the the label labfu of the vertex fu it returns the last
// int of the label of the last existing children.
// if no child exists, then it ouputs 0.
static proc lastchildrenindex(intvec labfu)
{
  int i;
  int lastlabi; intvec labi; intvec labfi;
  int lastlab=0;
  for (i=1;i<=size(@L);i++)
  {
    labi=@L[i][1];
    if (size(labi)>1)
    {
      labfi=delintvec(labi,size(labi));
      if (labfu==labfi)
      {
        lastlabi=labi[size(labi)];
        if (lastlab<lastlabi)
        {
          lastlab=lastlabi;
        }
      }
    }
  }
  return(lastlab);
}

// given the the vertex u it provides the next brother of u.
// if it does not exist, then it ouputs v=list(list(intvec(0)),0)
static proc nextbrother(intvec labu)
{
  list L; int i; int j; list next;
  int lastlabu=labu[size(labu)];
  intvec labfu=delintvec(labu,size(labu));
  int lastlabi; intvec labi; intvec labfi;
  for (i=1;i<=size(@L);i++)
  {
    labi=@L[i][1];
    if (size(labi)>1)
    {
      labfi=delintvec(labi,size(labi));
      if (labfu==labfi)
      {
        lastlabi=labi[size(labi)];
        if (lastlabu<lastlabi)
        {L[size(L)+1]=list(lastlabi,list(@L[i],i));}
      }
    }
  }
  if (size(L)==0){return(list(intvec(0),0));}
  next=L[1];
  for (i=2;i<=size(L);i++)
  {
    if (L[i][1]<next[1]){next=L[i];}
  }
  return(next[2]);
}

// gives the first child of vertex fu
static proc firstchild(labfu)
{
  intvec labfu0=labfu;
  labfu0[size(labfu0)+1]=0;
  return(nextbrother(labfu0));
}

// purpose: eliminate the children vertices of fu and all its descendents
// whose prime ideal Pu contains a prime ideal Pv of some brother vertex w.
static proc deletebrotherscontaining(intvec labfu)
{
  int i; int t;
  list fu=tree(labfu,@L);
  int childfu=fu[1][2];
  list u; intvec labu; ideal Pu;
  list v; intvec labv; ideal Pv;
  u=firstchild(labfu);
  for (i=1;i<=childfu;i++)
  {
    labu=u[1][1];
    Pu=u[1][3];
    v=firstchild(fu[1][1]);
    t=1;
    while ((t) and (v[2]!=0))
    {
      labv=v[1][1];
      Pv=v[1][3];
      if (labu!=labv)
      {
        if (containedP(Pv,Pu))
        {
          deleteverts(labu);
          fu=tree(labfu,@L);
          @L[fu[2]][2]=fu[1][2]-1;
          t=0;
        }
      }
      if (t!=0)
      {
        v=nextbrother(v[1][1]);
      }
    }
    if (i<childfu)
    {
      u=nextbrother(u[1][1]);
    }
  }
}

// purpose: delete all descendent vertices from u included u
// from the list @L.
// It must be noted that after the operation, the number of children
// in fathers vertex must be decreased in 1 unitity. This operation is not
// performed inside this recursive routine.
static proc deleteverts(intvec labu)
{
  int i; int ii; list v; intvec labv;
  list u=tree(labu,@L);
  int childu=u[1][2];
  @L=delete(@L,u[2]);
  if (childu>0)
  {
    v=firstchild(labu);
    labv=v[1][1];
    for (ii=1;ii<=childu;ii++)
    {
      deleteverts(labv);
      if (ii<childu)
      {
        v=nextbrother(v[1][1]);
        labv=v[1][1];
      }
    }
  }
}

// purpose: starting from vertex olab (initially nlab=olab)
// relabels the vertices of @L to be consecutive
static proc relabelingindices(intvec olab, intvec nlab)
{
  int i;
  intvec nlabi; intvec labv;
  list u=tree(olab,@L);
  int childu=u[1][2];
  list v;
  if (childu==0){@L[u[2]][1]=nlab;}
  else
  {
    v=firstchild(u[1][1]);
    @L[u[2]][1]=nlab;
    i=1;
    while(v[2]!=0)
    {
      labv=v[1][1];
      nlabi=nlab,i;
      relabelingindices(labv,nlabi);
      v=nextbrother(labv);
      i++;
    }
  }
}

// mrcgs
// input: F = ideal in ring R=Q[a][x]
// output: a list L representing the tree of the mrcgs.
static proc mrcgs(ideal F, list #)
//"USAGE:   mrcgs(F);
//          F is the ideal from which to obtain the Minimal Reduced CGS.
//          From the old library redcgs.lib.
//          Alternatively, as option:
//          mrcgs(F,L);
//          Options: We can give a list of options in the list L
//          of the form
//          ("null",ideal N,"nonnull",ideal W,"comment",0-1).
//          One can give none till 3 of these options by giving the
//          name of the option and the content.
//          When options "null" and/or "nonnull" are given, then the
//          parameter space is restricted to V(N)\V(h), where h is the product of
//          the non null polynomials in W. If the option ("comment",1) is set,
//          then information about the total number of segments of the
//          output is printed.
//          By default N=ideal(0), W=ideal(1), ("comment",0).
//          mrcgs is the fundamental routine of the old library redcgs.lib,
//          computing the minimal reduced comprehensive Groebner system.
//RETURN:   The list T representing the Minimal Reduced CGS.
//          The description given here is identical for rcgs and crcgs.
//          The elements of the list T computed by mrcgs are lists representing
//          a rooted tree.
//          Each element of the list T has the two first entries with the following content:
//           [1]: The label (intvec) representing the position in the rooted
//                tree:  0 for the root (and this is a special element)
//                       i for the root of the segment i
//                       (i,...) for the children of the segment i
//           [2]: the number of children (int) of the vertex.
//          There thus three kind of vertices:
//           (1) the root (first element labelled 0),
//           (2) the vertices labelled with a single integer i,
//           (3) the rest of vertices labelled with more indices.
//          Description of the root. Vertex type (1)
//           There is a special vertex (the first one) whose content is
//           the following:
//             [3] lpp of the given ideal
//             [4] the given ideal
//             [5] the red-representation  of the (optional) given null and non-null
//                 conditions (see redspec for the description).
//             [6] MRCGS (to remember which algorithm has been used). If the
//                 algorithm used is rcgs of crcgs then this will be stated
//                 at this vertex (RCGS or CRCGS).
//           Description of vertices type (2). These are the vertices that
//           initiate a segment, and are labelled with a single integer.
//             [3] lpp (ideal) of the reduced basis. If they are repeated lpp's this
//                 will correspond to a sheaf.
//             [4] the reduced basis (ideal) of the segment.
//           Description of vertices type (3). These vertices have as first
//           label i and descend form vertex i in the position of the label
//           (i,...). They contain moreover a unique prime ideal in the parameters
//           and form ascending chains of ideals.
//          How is to be read the mrcgs tree? The vertices with an even number of
//          integers in the label are to be considered as additive and those
//          with an odd number of integers in the label are to be considered as
//          substraction. As an example consider the following vertices:
//          v1=((i),2,lpp,B),
//          v2=((i,1),2,P_(i,1)),
//          v3=((i,1,1),2,P_(i,1,1)),
//          v4=((i,1,1,1),1,P_(i,1,1,1)),
//          v5=((i,1,1,1,1),0,P_(i,1,1,1,1)),
//          v6=((i,1,1,2),1,P_(i,1,1,2)),
//          v7=((i,1,1,2,1),0,P_(i,1,1,2,1)),
//          v8=((i,1,2),0,P_(i,1,2)),
//          v9=((i,2),1,P_(i,2)),
//          v10=((i,2,1),0,P_(i,2,1)),
//          They represent the segment:
//          (V(i,1)\(((V(i,1,1) \ ((V(i,1,1,1) \ V(i,1,1,1,1)) u (V(i,1,1,2) \ V(i,1,1,2,1)))))
//          u V(i,1,2))) u (V(i,2) \ V(i,2,1))
//          and can also be represented by
//          (V(i,1) \ (V(i,1,1) u V(i,1,2))) u
//          (V(i,1,1,1) \ V(i,1,1,1)) u
//          (V(i,1,1,2) \ V(i,1,1,2,1)) u
//          (V(i,2) \ V(i,2,1))
//          where V(i,j,..) = V(P_(i,j,..))
//NOTE:     There are three fundamental routines in the old library redcgs.lib:
//          mrcgs, rcgs and crcgs.
//          mrcgs (Minimal Reduced CGS) is an algorithm that packs so much as it
//          is able to do (using algorithms adhoc) the segments with the same lpp,
//          obtaining the minimal number of segments. The hypothesis is that this
//          is very close to be canonical, but there is no proof of the uniqueness
//          of this minimal packing. Moreover, the segments obtained are not
//          locally closed, i.e. there are not always the difference of two varieties,
//          but are a union of differences of varieties.
//          The output can be visualized using cantreetoMaple, that will
//          write a file with the content of mrcgs that can be read in Maple
//          and plotted using the Maple plotcantree routine of the Monte's dpgb library
//KEYWORDS: rcgs, crcgs, buildtree, cantreetoMaple,
//EXAMPLE:  mrcgs; shows an example"
{
  int i=1;
  int @ish=1;
  exportto(Top,@ish);
  while((@ish) and (i<=size(F)))
  {
    @ish=ishomog(F[i]);
    i++;
  }
  int comment=0;
  def N=ideal(0);
  def W=ideal(1);
  list L=#;
  for(i=1;i<=size(L)/2;i++)
  {
    if(L[2*i-1]=="null"){N=L[2*i];}
    else
    {
      if(L[2*i-1]=="nonnull"){W=L[2*i];}
      else
      {
        if(L[2*i-1]=="comment"){comment=L[2*i];}
      }
    }
  }
  def RR=basering;
  list LL=buildtree(F, #);
  setglobalrings();
  list S=selectcases(LL);
  list T=cantree(S);
  if(equalideals(N,ideal(0))==0)
  {
    T=reduceconds(T,N,W);
  }
  T[1][6]="MRCGS";
  T[1][4]=F;
  for (i=1;i<=size(F);i++)
  {
    T[1][3][i]=leadmonom(F[i]);
  }
  kill @ish;
  kill @P; kill @RP; kill @R;
  return(T);
}
//example
//{ "EXAMPLE:"; echo = 2;
//  ring R=(0,a1,a2,a3,a4),(x1,x2,x3,x4),dp;
//  ideal F=x4-a4+a2, x1+x2+x3+x4-a1-a3-a4, x1*x3*x4-a1*a3*a4, x1*x3+x1*x4+x2*x3+x3*x4-a1*a4-a1*a3-a3*a4;
//  "System="; F;
//  def T=mrcgs(F);
//  setglobalrings();
//  "mrcgs(F)="; T;
//  cantreetoMaple(T,"Tm","Tm.txt");
//  "cantodiffcgs(T)="; cantodiffcgs(T);
//  kill R;
//  ring R=(0,b,c,d,e,f),(x,y),dp;
//  ideal F1=x^2+b*y^2+2*c*x*y+2*d*x+2*e*y+f, 2*x+2*c*y+2*d, 2*b*y+2*c*x+2*e;
//  "System="; F1;
//  def T1=mrcgs(F1);
//  setglobalrings();
//  "mrcgs(F1)="; T1;
//  cantreetoMaple(T1,"T1m","T1m.txt");
//}

// reduceconds: when null and nonnull conditions are specified it
//              takes the output of cantree and reduces the tree
//              assuming the null and nonnull conditions
// input: list T (the output of cantree computed with null and nonull conditions
//        ideal N: null conditions
//        ideal W: non-null conditions
// output: the list T assuming the null and non-null conditions
static proc reduceconds(list T,ideal N,ideal W)
{
  int i; intvec lab; intvec labfu; list fu; int j; int t;
  list @L=T;
  exportto(Top,@L);
  int n=size(W);
  for (i=2;i<=size(@L);i++)
  {
    t=0; j=0;
    while ((not(t)) and (j<n))
    {
      j++;
      if (size(@L[i][1])>1)
      {
        if (memberpos(W[j],@L[i][3])[1])
        {
          t=1;
          @L[i][3]=ideal(1);
        }
      }
    }
  }
  for (i=2;i<=size(@L);i++)
  {
    if (size(@L[i][1])>1)
    {
      @L[i][3]=delidfromid(N,@L[i][3]);
    }
  }
  for (i=2;i<=size(@L);i++)
  {
    if ((size(@L[i][1])>1) and (size(@L[i][1]) mod 2==1) and (equalideals(@L[i][3],ideal(0))))
    {
      lab=@L[i][1];
      labfu=delintvec(lab,size(lab));
      fu=tree(labfu,@L);
      @L[fu[2]][2]=@L[fu[2]][2]-1;
      deleteverts(lab);
    }
  }
  for (j=2; j<=size(@L); j++)
  {
    if (@L[j][2]>0)
    {
      deletebrotherscontaining(@L[j][1]);
    }
  }
  for (i=1;i<=@L[1][2];i++)
  {
    relabelingindices(intvec(i),intvec(i));
  }
  list TT=@L;
  kill @L;
  return(TT);
}

//**************End of cantree******************************

//**************Begin of CanTreeTo Maple********************

// cantreetoMaple
// input:  list L: the output of cantree
//         string T: the name of the table of Maple that represents L
//                   in Maple
//         string writefile: the name of the file where the table T
//                           is written
proc cantreetoMaple(list L, string T, string writefile)
"USAGE:   cantreetoMaple(T, TM, writefile);
          T: is the list provided by grobcovold with option ("out",1),
          TM: is the name (string) of the table variable in Maple that will
             represent the output of the fundamental routines,
          writefile: is the name (string) of the file where to write the content.
RETURN:   writes the list provided by grobcovold to a file
          containing the table representing it in Maple.
NOTE:     It can be called from the output of grobcovold with option ("out",1)
KEYWORDS: grobcovold, Maple
EXAMPLE:  cantreetoMaple; shows an example"
{
  short=0;
  if(size(L[1])!=6)
  {
    "  'Warning!' grobcovold must be called with option 'out' set to 1 to be operative";
    return();
  }
  int i;
  def R=basering;
  list L0=L[1];
  int numcases=L0[2];
  link LLw=":w "+writefile;
  string La=string("table(",T,");");
  write(LLw, La);
  close(LLw);
  link LLa=":a "+writefile;
  def RL=ringlist(R);
  list p=RL[1][2];
  string param=string(p[1]);
  if (size(p)>1)
  {
    for(i=2;i<=size(p);i++){param=string(param,",",p[i]);}
  }
  list v=RL[2];
  string vars=string(v[1]);
  if (size(v)>1)
  {
    for(i=2;i<=size(v);i++){vars=string(vars,",",v[i]);}
  }
  list xord;
  list pord;
  if (RL[1][3][1][1]=="dp"){pord=string("tdeg(",param);}
  else
  {
    if (RL[1][3][1][1]=="lp"){pord=string("plex(",param);}
  }
  if (RL[3][1][1]=="dp"){xord=string("tdeg(",vars);}
  else
  {
    if (RL[3][1][1]=="lp"){xord=string("plex(",vars);}
  }
  write(LLa,string(T,"[[___xord]]:=",xord,");"));
  write(LLa,string(T,"[[___pord]]:=",pord,");"));
  //write(LLa,string(T,"[[11]]:=true; "));
  list S;
  S=string(T,"[[0]]:=",numcases,";");
  write(LLa,S);
  S=string(T,"[[___method]]:=",L[1][6],";");
  // Method L[1][6];
  write(LLa,S);
  S=string(T,"[[___basis]]:=[",L0[4],"];");
  write(LLa,S);
  S=string(T,"[[___nullcond]]:=[",L0[5][1],"];");
  write(LLa,S);
  S=string(T,"[[___notnullcond]]:={",L0[5][2],"};");
  write(LLa,S);
  for (i=1;i<=numcases;i++)
  {
    S=ctlppbasis(T,L,intvec(i));
    write(LLa,S[1]);
    write(LLa,S[2]);
    write(LLa,S[3]);
    //write(LLa,S[4]);
    ctrecwrite(LLa, L, T, intvec(i),S[4]);
  }
  close(LLa);
}
example
{ "EXAMPLE:"; echo = 2;
  ring R=(0,b,c,d,e,f),(x,y),dp;
  ideal F=x^2+b*y^2+2*c*x*y+2*d*x+2*e*y+f, 2*x+2*c*y+2*d, 2*b*y+2*c*x+2*e;
  def T=grobcovold(F,"out",1);
  T;
  cantreetoMaple(T,"Tm","Tm.txt");
}

// ctlppbasis: auxiliary cantreetoMaple routine
// input:
//   string T: the name of the table in Maple
//   intvec lab: the label of the case
//   ideal B: the basis of the case
// output:
//   the string of T[[lab]] (basis); in Maple
static proc ctlppbasis(string T, list L, intvec lab)
{
  list u;
  intvec lab0=lab,0;
  u=tree(lab,L);
  list Li;
  Li[1]=string(T,"[[",lab,",___lpp]]:=[",u[1][3],"]; ");
  Li[2]=string(T,"[[",lab,"]]:=[",u[1][4],"]; ");
  Li[3]=string(T,"[[",lab0,"]]:=",u[1][2],"; ");
  Li[4]=u[1][2];
  return(Li);
}

// ctlppbasis: auxiliary cantreetoMaple routine
// recursive routine to write all elements
static proc ctrecwrite(LLa, list L, string T, intvec lab, int n)
{
  int i;
  intvec labi; intvec labi0;
  string S;
  list u;
  for (i=1;i<=n;i++)
  {
    labi=lab,i;
    u=tree(labi,L);
    S=string(T,"[[",labi,"]]:=[",u[1][3],"];");
    write(LLa,S);
    labi0=labi,0;
    S=string(T,"[[",labi0,"]]:=",u[1][2],";");
    write(LLa,S);
    ctrecwrite(LLa, L, T, labi, u[1][2]);
  }
}

//**************End of CanTreeTo Maple********************

//**************Begin homogenizing************************

// ishomog:
// Purpose: test if a polynomial is homogeneous in the variables or not
// input:  poly f
// output  1 if f is homogeneous, 0 if not
static proc ishomog(f)
{
  int i; poly r; int d; int dr;
  if (f==0){return(1);}
  d=deg(f); dr=d; r=f;
  while ((d==dr) and (r!=0))
  {
    r=r-lead(r);
    dr=deg(r);
  }
  if (r==0){return(1);}
  else{return(0);}
}

static proc rcgs(ideal F, list #)
//"USAGE:   rcgs(F);
//          F is the ideal from which to obtain the Reduced CGS.
//          From the old library redcgs.lib.
//          Alternatively, as option:
//          rcgs(F,L);
//          Options: We can give a list of options in the list L
//          of the form
//          ("null",ideal N,"nonnull",ideal W,"comment",int comment).
//          One can give none till 3 of these options by giving the
//          name of the option and the content.
//          When options "null" and/or "nonnull" are given, then the
//          parameter space is restricted to V(N)\V(h), where h is the product of
//          the non null polynomials in W. If the option "comment" is set to 1,
//          then information about the total number of segments of the
//          output is printed.
//          By default N=ideal(0) and W=ideal(1).
//          rcgs is the a routine whose output segments are always
//          locally closed and correspond to homogenizing the basis
//          compute its mrcgs and then reduce and de-homogenizing the result.
//          The result is a Reduced Comprehensive Groebner System.
//RETURN:   The list T representing the Reduced CGS.
//          The description given here is identical for mrcgs and crcgs.
//          The elements of the list T computed by rcgs are lists representing
//          a rooted tree.
//          Each element of the list T has the two first entries with the following content:
//           [1]: The label (intvec) representing the position in the rooted
//                tree:  0 for the root (and this is a special element)
//                       i for the root of the segment i
//                       (i,...) for the children of the segment i
//           [2]: the number of children (int) of the vertex.
//          There thus three kind of vertices:
//           (1) the root (first element labelled 0),
//           (2) the vertices labelled with a single integer i,
//           (3) the rest of vertices labelled with more indices.
//          Description of the root. Vertex type (1)
//           There is a special vertex (the first one) whose content is
//           the following:
//             [3] lpp of the given ideal
//             [4] the given ideal
//             [5] the red-representation  of the (optional) given null and non-null conditions
//                 (see redspec for the description)
//             [6] RCGS (to remember which algorithm has been used). If the
//                 algorithm used is mrcgs of crcgs then this will be stated
//                 at this vertex (MRCGS or CRCGS).
//           Description of vertices type (2). These are the vertices that
//           initiate a segment, and are labelled with a single integer.
//             [3] lpp (ideal) of the reduced basis. If they are repeated lpp's this
//                 will correspond to a sheaf.
//             [4] the reduced basis (ideal) of the segment.
//           Description of vertices type (3). These vertices have as first
//           label i and descend form vertex i in the position of the label
//           (i,...). They contain moreover a unique prime ideal in the parameters
//           and form ascending chains of ideals.
//          How is to be read the rcgs tree? The vertices with an even number of
//          integers in the label are to be considered as additive and those
//          with an odd number of integers in the label are to be considered as
//          substraction. As an example consider the following vertices:
//          v1=((i),2,lpp,B),
//          v2=((i,1),2,P_(i,1)),
//          v3=((i,1,1),2,P_(i,1,1)),
//          v4=((i,1,1,1),1,P_(i,1,1,1)),
//          v5=((i,1,1,1,1),0,P_(i,1,1,1,1)),
//          v6=((i,1,1,2),1,P_(i,1,1,2)),
//          v7=((i,1,1,2,1),0,P_(i,1,1,2,1)),
//          v8=((i,1,2),0,P_(i,1,2)),
//          v9=((i,2),1,P_(i,2)),
//          v10=((i,2,1),0,P_(i,2,1)),
//          They represent the segment:
//          (V(i,1)\(((V(i,1,1) \ ((V(i,1,1,1) \ V(i,1,1,1,1)) u (V(i,1,1,2) \ V(i,1,1,2,1)))))
//          u V(i,1,2))) u (V(i,2) \ V(i,2,1))
//          and can also be represented by
//          (V(i,1) \ (V(i,1,1) u V(i,1,2))) u
//          (V(i,1,1,1) \ V(i,1,1,1)) u
//          (V(i,1,1,2) \ V(i,1,1,2,1)) u
//          (V(i,2) \ V(i,2,1))
//          where V(i,j,..) = V(P_(i,j,..))
//NOTE:     There are three fundamental routines in the old library redcgs.lib:
//          mrcgs, rcgs and crcgs.
//          The output can be visualized using cantreetoMaple, that will
//          write a file with the content of rcgs that can be read in Maple
//          and plotted using the Maple plotcantree routine of the Monte's dpgb library
//KEYWORDS: mrcgs, crcgs, buildtree, cantreetoMaple,
//EXAMPLE:  rcgs; shows an example"
{
  int j; int i;
  poly f;
  int comment=0;
  def N=ideal(0);
  def W=ideal(1);
  list L=#;
  for(i=1;i<=size(L)/2;i++)
  {
    if(L[2*i-1]=="null"){N=L[2*i];}
    else
    {
      if(L[2*i-1]=="nonnull"){W=L[2*i];}
      else
      {
        if(L[2*i-1]=="comment"){comment=L[2*i];}
      }
    }
  }
  i=1; int postred=0;
  int ish=1;
  while ((ish) and (i<=size(F)))
  {
    ish=ishomog(F[i]);
    i++;
  }
  if (ish){return(mrcgs(F, #));}
  def RR=basering;
  list RRL=ringlist(RR);
  //if (RRL[3][1][1]!="dp"){ERROR("the order must be dp");}
  poly @t;
  ring H=0,@t,dp;
  def RH=RR+H;
  setring(RH);
  def FH=imap(RR,F);
  list u; ideal B; ideal lpp; intvec lab;
  FH=homog(FH,@t);
  def Nh=imap(RR,N);
  def Wh=imap(RR,W);
  list LL;
  if ((size(Nh)>0) or (size(Wh)>0))
  {
    LL=mrcgs(FH,list("null",Nh,"nonnull",Wh));
  }
  else
  {
    LL=mrcgs(FH);
  }
  setglobalrings();
  LL[1][3]=subst(LL[1][3],@t,1);
  LL[1][4]=subst(LL[1][4],@t,1);
  for (i=1; i<=LL[1][2]; i++)
  {
    lab=intvec(i);
    u=tree(lab,LL);
    postred=difflpp(u[1][3]);
    B=sortideal(subst(LL[u[2]][4],@t,1));
    lpp=sortideal(subst(LL[u[2]][3],@t,1));
    if (memberpos(1,B)[1]){B=ideal(1); lpp=ideal(1);}
    if (postred)
    {
      lpp=ideal(0);
      B=postredgb(mingb(B));
      for (j=1;j<=size(B);j++){lpp[j]=leadmonom(B[j]);}
    }
    else{"Sheaves present, not reduced bases in the case lpp = ";lpp;}
    LL[u[2]][4]=B;
    LL[u[2]][3]=lpp;
  }
  setring(RR);
  list LLL=imap(RH,LL);
  kill @P; kill @R; kill @RP;
  LLL[1][6]="RCGS";
  return(LLL);
}
//example
//{ "EXAMPLE:"; echo = 2;
//  ring R=(0,b,c,d,e,f),(x,y),dp;
//  ideal F=x^2+b*y^2+2*c*x*y+2*d*x+2*e*y+f, 2*x+2*c*y+2*d, 2*b*y+2*c*x+2*e;
//  def T=rcgs(F);
//  T;
//  cantreetoMaple(T,"Tr","Tr.txt");
//  cantodiffcgs(T);
//}

static proc difflpp(ideal lpp)
{
  int t=1; int i;
  poly lp1=lpp[1];
  poly lp;
  i=2;
  while ((i<=size(lpp)) and (t))
  {
    lp=lpp[i];
    if (lp==lp1){t=0;}
    lp1=lp;
    i++;
  }
  return(t);
}

// redgb: given a minimal bases (gb reducing) it
// reduces each polynomial wrt to the others
static proc postredgb(ideal F)
{
  ideal G;
  ideal H;
  int i;
  if (size(F)==0){return(ideal(0));}
  for (i=1;i<=size(F);i++)
  {
    H=delfromideal(F,i);
    G[i]=pdivi(F[i],H)[1];
  }
  return(G);
}

static proc crcgs(ideal F, list #)
//"USAGE:   crcgs(F);
//          F is the ideal from which to obtain the Canonical Reduced CGS.
//          From the old library redcgs.lib.
//          Alternatively, as option:
//          crcgs(F,L);
//          Options: We can give a list of options in the list L
//          of the form
//          ("null",ideal N,"nonnull",ideal W,"comment",int comment).
//          One can give none till 3 of these options by giving the
//          name of the option and the content.
//          When options "null" and/or "nonnull" are given, then the
//          parameter space is restricted to V(N)\V(h), where h is the product of
//          the non null polynomials in W. If the option "comment" is set to 1,
//          then information about the total number of segments of the
//          output is printed.
//          By default N=ideal(0) and W=ideal(1).
//          crcgs is a routine whose output segments are always
//          locally closed and correspond to homogenizing the ideal
//          compute its mrcgs and then reduce and de-homogenizing the result.
//          The result is in principle the Canonical Comprehensive Groebner System,
//          similar to the result obtained by the fundamental routine grobcov,
//          but the output is less friendly and not certified to be always
//          the canonical Groebner cover.
//RETURN:   The list T representing the canonical Reduced CGS.
//          The description given here is identical for mrcgs and rcgs.
//          The elements of the list T computed by crcgs are lists representing
//          a rooted tree.
//          Each element of the list T has the two first entries with the following content:
//           [1]: The label (intvec) representing the position in the rooted
//                tree:  0 for the root (and this is a special element)
//                       i for the root of the segment i
//                       (i,...) for the children of the segment i
//           [2]: the number of children (int) of the vertex.
//          There thus three kind of vertices:
//           (1) the root (first element labelled 0),
//           (2) the vertices labelled with a single integer i,
//           (3) the rest of vertices labelled with more indices.
//          Description of the root. Vertex type (1)
//           There is a special vertex (the first one) whose content is
//           the following:
//             [3] lpp of the given ideal
//             [4] the given ideal
//             [5] the red-representation  of the (optional) given null and non-null conditions
//                 (see redspec for the description)
//             [6] CRCGS (to remember which algorithm has been used). If the
//                 algorithm used is mrcgs of rcgs then this will be stated
//                 at this vertex (MRCGS or RCGS).
//           Description of vertices type (2). These are the vertices that
//           initiate a segment, and are labelled with a single integer.
//             [3] lpp (ideal) of the reduced basis. If they are repeated lpp's this
//                 will correspond to a sheaf.
//             [4] the reduced basis (ideal) of the segment.
//           Description of vertices type (3). These vertices have as first
//           label i and descend form vertex i in the position of the label
//           (i,...). They contain moreover a unique prime ideal in the parameters
//           and form ascending chains of ideals.
//          How is to be read the crcgs tree? The vertices with an even number of
//          integers in the label are to be considered as additive and those
//          with an odd number of integers in the label are to be considered as
//          substraction. As an example consider the following vertices:
//          v1=((i),2,lpp,B),
//          v2=((i,1),2,P_(i,1)),
//          v3=((i,1,1),2,P_(i,1,1)),
//          v4=((i,1,1,1),1,P_(i,1,1,1)),
//          v5=((i,1,1,1,1),0,P_(i,1,1,1,1)),
//          v6=((i,1,1,2),1,P_(i,1,1,2)),
//          v7=((i,1,1,2,1),0,P_(i,1,1,2,1)),
//          v8=((i,1,2),0,P_(i,1,2)),
//          v9=((i,2),1,P_(i,2)),
//          v10=((i,2,1),0,P_(i,2,1)),
//          They represent the segment:
//          (V(i,1)\(((V(i,1,1) \ ((V(i,1,1,1) \ V(i,1,1,1,1)) u (V(i,1,1,2) \ V(i,1,1,2,1)))))
//          u V(i,1,2))) u (V(i,2) \ V(i,2,1))
//          and can also be represented by
//          (V(i,1) \ (V(i,1,1) u V(i,1,2))) u
//          (V(i,1,1,1) \ V(i,1,1,1)) u
//          (V(i,1,1,2) \ V(i,1,1,2,1)) u
//          (V(i,2) \ V(i,2,1))
//          where V(i,j,..) = V(P_(i,j,..))
//NOTE:     There are three fundamental routines in the old library redcgs.lib:
//          mrcgs, rcgs and crcgs.
//          The output can be visualized using cantreetoMaple, that will
//          write a file with the content of rcgs that can be read in Maple
//          and plotted using the Maple plotcantree routine of the Monte's dpgb library
//KEYWORDS: mrcgs, crcgs, buildtree, cantreetoMaple,
//EXAMPLE:  rcgs; shows an example"
{
  int ish=1; int i=1;
  while ((ish) and (i<=size(F)))
  {
    ish=ishomog(F[i]);
    i++;
  }
  if (ish){return(mrcgs(F, #));}
  def RR=basering;
//  int comment=0;
//  def N=ideal(0);
//  def W=ideal(1);
//  list L=#;
//  for(i=1;i<=size(L)/2;i++)
//  {
//    if(L[2*i-1]=="null"){N=L[2*i];}
//    else
//    {
//      if(L[2*i-1]=="nonnull"){W=L[2*i];}
//      else
//      {
//        if(L[2*i-1]=="comment"){comment=L[2*i];}
//      }
//    }
//  }
  setglobalrings();
  setring(@RP);
  ideal FP=imap(RR,F);
  option(redSB);
  def G=std(FP);
  setring(RR);
  def GR=imap(@RP,G);
  kill @P; kill @RP; kill @R;
  list LL;
  LL=rcgs(GR, #);
  LL[1][6]="CRCGS";
  return(LL);
}
//example
//{ "EXAMPLE:"; echo = 2;
//  ring R=(0,b,c,d,e,f),(x,y),dp;
//  ideal F=x^2+b*y^2+2*c*x*y+2*d*x+2*e*y+f, 2*x+2*c*y+2*d, 2*b*y+2*c*x+2*e;
//  def T=crcgs(F);
//  T;
//  cantreetoMaple(T,"Tc","Tc.txt");
//  cantodiffcgs(T);
//}

//purpose ideal intersection called in @R and computed in @P
static proc idintR(ideal N, ideal M)
{
  def RR=basering;
  setring(@P);
  def Np=imap(RR,N);
  def Mp=imap(RR,M);
  def Jp=idint(Np,Mp);
  setring(RR);
  return(imap(@P,Jp));
}

//purpose reduced groebner basis called in @R and computed in @P
static proc gbR(ideal N)
{
  def RR=basering;
  setring(@P);
  def Np=imap(RR,N);
  option(redSB);
  Np=std(Np);
  setring(RR);
  return(imap(@P,Np));
}

// purpose: given the output of a locally closed CGS (i.e. from rcgs or crcgs)
//          it returns the segments as difference of varieties.
static proc cantodiffcgs(list L)
//"USAGE:   canttodiffcgs(T);
//          T: is the list provided by mrcgs or crcgs or crcgs,
//RETURN:   The list transforming the content of these routines to a simpler
//          output where each segment corresponds to a single element of the list
//          that is described as difference of two varieties.
//
//          The first element of the list is identical to the first element
//          of the list provided by the corresponding cgs algorithm, and
//          contains general information on the call (see mrcgs).
//          The remaining elements are lists of 4 elements,
//          representing segments. These elements are
//           [1]: the lpp of the segment
//           [2]: the basis of the segment
//           [3]; the ideal of the first variety (radical)
//           [4]; the ideal of the second variety (radical)
//          The segment is V([3]) \ V([4]).
//
//NOTE:     It can be called from the output of mrcgs or rcgs of crcgs
//KEYWORDS: mrcgs, rcgs, crcgs, Maple
//EXAMPLE:  cantodiffcgs; shows an example"
{
  int i; int j; int k; int depth; list LL; list u; list v; list w;
  ideal N; ideal Nn; ideal M; ideal Mn; ideal N0; ideal W0;
  LL[1]=L[1];
  N0=L[1][5][1];
  W0=L[1][5][2];
  def RR=basering;
  setring(@P);
  def N0P=imap(RR,N0);
  def W0P=imap(RR,N0);
  ideal NP;
  ideal MP;
  setring(RR);
  for (i=2;i<=size(L);i++)
  {
    depth=size(L[i][1]);
    if (depth>3){ERROR("the given CGS has non locally closed segments");}
  }
  for (i=1;i<=L[1][2];i++)
  {
    N=ideal(1);
    M=ideal(1);
    u=tree(intvec(i),L);
    for (j=1;j<=u[1][2];j++)
    {
      v=tree(intvec(i,j),L);
      Nn=v[1][3];
      N=idintR(N,Nn);
      for (k=1;k<=v[1][2];k++)
      {
        w=tree(intvec(i,j,k),L);
        Mn=w[1][3];
        M=idintR(M,Mn);
      }
    }
    setring(@P);
    NP=imap(RR,N);
    MP=imap(RR,M);
    MP=MP+N0P;
    for (j=1;j<=size(W0P);j++){MP=MP+ideal(W0P[j]);}
    NP=NP+N0P;
    NP=gbR(NP);
    MP=gbR(MP);
    setring(RR);
    N=imap(@P,NP);
    M=imap(@P,MP);
    LL[i+1]=list(u[1][3],u[1][4],N,M);
  }
  return(LL);
}
//example
//{ "EXAMPLE:"; echo = 2;
//  ring R=(0,b,c,d,e,f),(x,y),dp;
//  ideal F=x^2+b*y^2+2*c*x*y+2*d*x+2*e*y+f, 2*x+2*c*y+2*d, 2*b*y+2*c*x+2*e;
//  def T=crcgs(F);
//  T;
//  cantreetoMaple(T,"Tc","Tc.txt");
//  cantodiffcgs(T);
//}

//**************End homogenizing************************

//**************End of redcgs************************

//**************Begin of Groebner Cover*****************

// incquotient
// incremental quotient
// Input: ideal N: a Groebner basis of an ideal
//        poly f:
// Output: Na = N:<f>
static proc incquotient(ideal N, poly f)
{
  poly g; int i;
  ideal Nb; ideal Na=N;

  // begins incquotient
  if (size(Na)==1)
  {
    g=gcd(Na[1],f);
    if (g!=1)
    {
      Na[1]=Na[1]/g;
    }
    attrib(Na,"IsSB",1);
    return(Na);
  }
  def P=basering;
  poly @t;
  ring H=0,@t,lp;
  def HP=H+P;
  setring(HP);
  def fh=imap(P,f);
  def Nh=imap(P,N);
  ideal Nht;
  for (i=1;i<=size(Nh);i++)
  {
    Nht[i]=Nh[i]*@t;
  }
  attrib(Nht,"isSB",1);
  def fht=(1-@t)*fh;
  option(redSB);
  Nht=std(Nht,fht);
  ideal Nc; ideal v;
  for (i=1;i<=size(Nht);i++)
  {
    v=variables(Nht[i]);
    if(memberpos(@t,v)[1]==0)
    {
      Nc[size(Nc)+1]=Nht[i]/fh;
    }
  }
  setring(P);
  ideal HH;
  def Nd=imap(HP,Nc); Nb=Nd;
  option(redSB);
  Nb=std(Nd);
  return(Nb);
}

// RrepNN: given a red-representation of a locally closed set and a new
//         assumed non-null polynomial f, it returns the new R-representation.
//         Called in any @P
//         13/09/2010
// input:
//   ideal N : the ideal of null-conditions
//   ideal W : non-null set of polynomials. (N,W) is a R-representation of the
//             initial locally closed set.
//   poly f  : A new assumed non-null polynomial
// returns: list (N1,W1), the new R-representation:
//   N1 = new radical of the null conditions of the R-representation
//   W1 = non-null list of polynomials of the new R-representation.
//   If the given conditions are not compatible, then N1=ideal(1). This should not
//     happen, because this has to be tested before using RrepNN.

static proc RrepNN(ideal N, ideal W, poly f)
//"USAGE:   RrepNN(N,W,f);
//          N: null conditions ideal of the initial R-representation
//          W: non-null list of polynomials of the initial R-representation
//          f: new assumed non-null polynomial
//RETURN:   a list (N1,W1) containing the new R-representation of the segment
//          (N,W) adding the new non-null condition f.
//NOTE:     Called from parameter ring (@P).
//KEYWORDS: representation
//EXAMPLE:  RrepNN; shows an example"
{
  ideal F=f; ideal W1=W;
  def N1=incquotient(N,f);
  option(redSB);
  N1=std(N1);
  //attrib(N1,"IsSB",1);
  def H=sqrfree(f);
  int i;
  for(i=1;i<=size(H);i++){W1[size(W1)+1]=reduce(H[i],N1);}

  W1=facvar(W1);
  if (size(W1)==0){W1=1;}
  return(list(N1,W1));
}
//example
//{ "EXAMPLE:"; echo = 2;
//  ring r=(0,a,b,c),(x,y),dp;
//  setglobalrings();
//  ideal N=(ab-c)*(a-b),(a-bc)*(a-b);
//  poly  h=(a+b)bc;
//  poly f=a-b;
//}

// RrepN:  given a red-representation of a locally closed set and a new
//         assumed null polynomial f, that is not identically null, it returns
//         the new red-representation.
//         Called in ring @P
//         13/09/2010
// input:
//   ideal N : the ideal of null-conditions
//   ideal W : non-null list of polynomials. (N,W) is a R-representation of the
//     initial locally closed set.
//   poly f  : A new assumed null polynomial
// returns: list (N1,W1), the new R-representation:
//   N1 = new radical of the null conditions of the R-representation
//   W1 = non-null list of polynomials of the new R-representation.
//   If the given conditions are not compatible, then N1=ideal(1).
static proc RrepN(ideal N, ideal W, poly f)
//"USAGE:   RrepN(N,W,f);
//          N: null conditions ideal of the initial R-representation
//          W: non-null list of polynomials of the initial R-representation
//          f: new assumed null polynomial
//RETURN:   a list (N1,W1) containing the new R-representation of the segment
//          (N,W) adding the new non-null condition f.
//NOTE:     Called from parameter ring (@P).
//KEYWORDS: representation
//EXAMPLE:  RrepN; shows an example"
{
  attrib(N,"isSB",1);
  def N1=std(N,f);
  option(redSB);
  N1=std(radical(N1));
  int i;
  poly h;
  for (i=1;i<=size(W);i++)
  {
    h=W[i];
    N1=incquotient(N1,h);
  }
  option(redSB);
  N1=std(N1);
  def W1=W;
  if (size(W1)==0){W1=1;}
  return(list(N1,W1));
}
//example
//{ "EXAMPLE:"; echo = 2;
//  ring r=(0,a,b,c),(x,y),dp;
//  setglobalrings();
//  ideal N=(ab-c)*(a-b),(a-bc)*(a-b);
//  poly  h=(a+b)bc;
//  poly f=a-b;
//  RrepN(N,h,f);
//}

// Rrep: generates a R-representation
//       called from any ring
//       it uses ring @P, thus the globalrings @P, @RP, @R must be
//       active by a previous call to setglobalrings();
//       13/09/2010
// input:
//   ideal N : the ideal of null-conditions (not necessarily radical nor canonical)
//   ideal W : set of non-null polynomials: if W corresponds to no non null
//             conditions then W=ideal(0)
//             otherwise it should be given as an ideal.
// returns: list (Na,Wa)
//   the R-representation of (N,W):
//   ideal Na = radical of the R-representation (canonical)
//   ideal Wa = set of non-null polynomials in the R-representation.
//             if it corresponds to no non null conditions then it is ideal(0)
//             otherwise the ideal is returned.
//   If the given conditions are not compatible, then N=ideal(1).
static proc Rrep(ideal Ni, ideal Wi)
//"USAGE:   Rrep(N,W);
//          N: null conditions ideal
//          W: set of non-null polynomials (ideal)
//RETURN:   a list (N1,W1) containing the R-representation of the segment (N,W).
//          N1 is the radical reduced ideal characterizing the segment.
//          V(N1) is the Zarisky closure of the segment (N,W).
//          The segment S=V(N1) \ V(h), where h=prod(w in W1)
//          N1 is uniquely determined and no prime component of N1 contains none of
//          the polynomials in W1.
//NOTE:     Can be called from ring @R but it works in ring @P. Thus
//          the globalrings @P, @RP, @R must be active by a previous call
//          to setglobalrings();
//KEYWORDS: R-representation
//EXAMPLE:  Rrep shows an example"
{
  def RR=basering;
  setring(@P);
  def N=imap(RR,Ni);
  option(redSB);
  N=std(radical(N));
  def W=imap(RR,Wi);
  if(size(W)==0){W=ideal(0);}
     //when there are no non-null conditions then W=ideal(1)
  else
  {
    W=facvar(W);
  }
  if (size(W)==0)
  {
    setring(RR);
    //def Wb=imap(@P,W);
    return(list(imap(@P,N), ideal(1)));
  }
  else
  {
    int i; //ideal F;
    for (i=1;i<=size(W);i++)
    {
      //F=W[i];
      N=incquotient(N,W[i]);
    }
    option(redSB);
    N=std(N);
    setring(RR);
    def Nb=imap(@P,N);
    def Wb=imap(@P,W);
    if (equalideals(Wb,ideal(0))){Wb=ideal(1);}
    return(list(Nb,Wb));
  }
}
//example
//{ "EXAMPLE:"; echo = 2;
//  ring R=(0,a,b,c),(x,y),dp;
//  setglobalrings();
//  ideal N=(ab-c)*(a-b),(a-bc)*(a-b);
//  ideal W=a^2-b^2,bc;
//  Rrep(N,W);
//}


// eliminate the ith element from a list
static proc elimfromlist(list l, int i)
{
  list L; int j;
  for(j=1;j<=i-1;j++)
  {L[j]=l[j];}
  for(j=i+1;j<=size(l);j++)
  {L[j-1]=l[j];}
  return(L);
}

static proc idbefid(ideal a, ideal b)
{
  poly fa; poly fb; poly la; poly lb;
  int te=1; int i; int j;
  int na=size(a);
  int nb=size(b);
  int nm;
  if (na<=nb){nm=na;} else{nm=nb;}
  for (i=1;i<=nm; i++)
  {
    fa=a[i]; fb=b[i];
    while((fa!=0) or (fb!=0))
    {
      la=lead(fa);
      lb=lead(fb);
      fa=fa-la;
      fb=fb-lb;
      la=leadmonom(la);
      lb=leadmonom(lb);
      if(leadmonom(la+lb)!=la){return(1);}
      else{if(leadmonom(la+lb)!=lb){return(2);}}
    }
  }
  if(na<nb){return(1);} else{if(na>nb){return(2);} else{return(0);}}
}

static proc sortlistideals(list L)
{
  int i; int j; int n;
  ideal a; ideal b;
  list LL=L;
  list NL;
  int k; int te;
  i=1;
  while(size(LL)>0)
  {
    k=1;
    for(j=2;j<=size(LL);j++)
    {
      te=idbefid(LL[k],LL[j]);
      if (te==2){k=j;}
    }
    NL[size(NL)+1]=LL[k];
    n=size(LL);
    if (n>1){LL=elimfromlist(LL,k);} else{LL=list();}
  }
  return(NL);
}

// returns 1 if the two lists of ideals are equal and 0 if not
static proc equallistideals(list L, list M)
{
  int t; int i;
  if (size(L)!=size(M)){return(0);}
  else
  {
    t=1;
    if (size(L)>0)
    {
      i=1;
      while ((t==1) and (i<=size(L)))
      {
        if (equalideals(L[i],M[i])==0){t=0;}
        i++;
      }
    }
    return(t);
  }
}

// RtoPrepNew
// Computes the P-representaion of a R-representaion (N,W) of a set
// input:
//    ideal N (null conditions, must be radical)
//    ideal W (non-null conditions ideal)
//    list L  must contain the radical decomposition of N.
// output:
//    the ((p_1,(p_11,..,p_1k_1)),..,(p_r,(p_r1,..,p_rk_r)));
//    the Prep of V(N) \ V(h), where h=prod(w in W).
static proc RtoPrepNew(ideal N, ideal W)
{
  int i; int j; list L0;
  if (N[1]==1)
  {
    L0[1]=list(ideal(1),list(ideal(1)));
    return(L0);
  }
  def RR=basering;
  setring(@P);
  ideal Np=imap(RR,N);
  ideal Wp=imap(RR,W);
  list Lp=minGTZ(Np);
  for(i=1;i<=size(Lp);i++)
  {
    option(redSB);
    Lp[i]=std(Lp[i]);
  }
  //list Lp=imap(RR,L);
  poly h=1;
  for (i=1;i<=size(Wp);i++){h=h*Wp[i];}
  list r; list Ti; list LL;
  for (i=1;i<=size(Lp);i++)
  {
    Ti=minGTZ(Lp[i]+h);
    for(j=1;j<=size(Ti);j++)
    {
      option(redSB);
      Ti[j]=std(Ti[j]);
    }
    //list LL[i];
    LL[i]=list(Lp[i],Ti);
  }
  setring(RR);
  return(imap(@P,LL));
}

// splitR: a new leading coefficient f is given to a R-representation
//        then splitR computes the two new R-representation by
//        considering it null, and non null.
//        Can be called from any ring but it works in ring @P
//        14/09/2010
// given the R-representation (N,W) and a new`polynomial f,
//        it outputs the null and the non-null R-representations adding f.
//        if the output R-representation (N0,W0) has N0==ideal(1) then
//        there must be no split and recbtcgs must continue on
//        the compatible (N1,W1) R-representation.
// input:
//    ideal N: null-ideal of the R-representation
//    ideal W: non-null list of polynomials of the R-representation
//    poly f coefficient to split if needed
// output:
//    list L = (list(ideal N0, ideal W0), list(ideal N1, ideal W1))
static proc splitR(ideal Ni, ideal Wi, poly fi)
{
  def RR=basering;
  setring(@P);
  def f=imap(RR,fi);
  def N=imap(RR,Ni);
  def W=imap(RR,Wi);
  def L0=RrepN(N,W,f);
  if(L0[1][1]==1)
  {
    setring(RR);
    def LL0=list(ideal(1),ideal(1));
    list LL1=list(Ni,Wi);
    return(list(LL0,LL1));
  }
  else
  {
    def L1=RrepNN(N,W,f);
    setring(RR);
    def LL0=imap(@P,L0);
    def LL1=imap(@P,L1);
    return(list(LL0,LL1));
  }
}

// Prep
// Computes the P-representation of V(N) \ V(M).
// input:
//    ideal N (null ideal) (not necessarily radical nor maximal)
//    ideal M (hole ideal) (not necessarily containing N)
// output:
//    the ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r)));
//    the Prep of V(N)\V(M)
// Assumed to work in the ring @P of the parameters
static proc Prep(ideal N, ideal M)
{
  if (N[1]==1)
  {
    //L0=list(list(ideal(1),list(ideal(1))));
    return(list(list(ideal(1),list(ideal(1)))));
  }
  def RR=basering;
  setring(@P);
  ideal Np=imap(RR,N);
  ideal Mp=imap(RR,M);
  int i; int j; list L0;

  list Ni=minGTZ(Np);
  list prep;
  for(j=1;j<=size(Ni);j++)
  {
    option(redSB);
    Ni[j]=std(Ni[j]);
  }
  list Mij;
  for (i=1;i<=size(Ni);i++)
  {
    Mij=minGTZ(Ni[i]+Mp);
    for(j=1;j<=size(Mij);j++)
    {
      option(redSB);
      Mij[j]=std(Mij[j]);
    }
    if ((size(Mij)==1) and (equalideals(Ni[i],Mij[1])==1)){;}
    else
    {
        prep[size(prep)+1]=list(Ni[i],Mij);
    }
  }
  if (size(prep)==0){prep=list(list(ideal(1),list(ideal(1))));}
  setring(RR);
  return(imap(@P,prep));
}

// PtoCrep
// Computes the C-representation from the P-representation.
// input:
//    list ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r)));
//         the P-representation of V(N)\V(M)
// output:
//    list (ideal ida, ideal idb)
//    the C-represen taion of V(N)\V(M) = V(ida)\V(idb)
// Assumed to work in the ring @P of the parameters
static proc PtoCrep(list L)
{
  def RR=basering;
  setring(@P);
  def Lp=imap(RR,L);
  int i; int j;
  ideal ida=ideal(1); ideal idb=ideal(1); list Lb; ideal N;
  for (i=1;i<=size(Lp);i++)
  {
    option(returnSB);
    N=Lp[i][1];
    ida=intersect(ida,N);
    Lb=Lp[i][2];
    for(j=1;j<=size(Lb);j++)
    {
      idb=intersect(idb,Lb[j]);
    }
  }
  def La=list(ida,idb);
  setring(RR);
  return(imap(@P,La));
}

// addnewpairs:
// 14/09/2010
// input:
//    ideal F, the given ideal
//    list P: the list of existing pairs to be computed
//    int l (the new index to add S-pols)
// output: list of ordered pairs (i,j,lcmij) of F in ascending order of lcmij
//         adding the new (i,l,lcmil) and placing them in order of ascending lcm
//         if a pair verifies Buchberger 1st criterion it is not stored
// ring @R
static proc addnewpairs(ideal F, list P, int l)
{
  int i;
  poly lm;
  poly lpf;
  poly lpg;
  list P1=P;
  list pair;
  if (size(F)<=1){return(P);}
  for (i=1;i<l;i++)
  {
    lm=lcmlmonoms(F[i],F[l]);
    // Buchberger 1st criterion
    lpf=leadmonom(F[i]);
    lpg=leadmonom(F[l]);
    if (lpf*lpg!=lm)
    {
      pair=(i,l,lm);
      P1=placepairinlist(pair,P1);
    }
  }
  return(P1);
}

// DiscussPolys: given the data in a vertex of btcgs (BuildTree), it analyzes the
//               leadcoef of the polynomials in B until it finds
//               that one of them can be either null or non null.
//               In that case, recbtcgs has to split into two branches, and then
//               l < size(B)
//               If not, and at the end only the non null option is compatible
//               then the reduced B has all the leadcoef non null, and then l=size(B).
//               15/09/2010
// ring @R
// input:
//    B:   (ideal) the actual basis
//    N:   (ideal) null conditions (R-rep)
//    W:   (ideal) non-null conditions set (R-rep)
//    P:   (list) of pairs of indices of S-polynomials that can and must be computed
//         (its leading coefficients are non-null, and using Buchberger's
//          criterions they are to be computed)
//    l:   (integer) representing the last polynomial in B for which the leading
//         coefficient is already assumed non-null.
// output: list of (cond,lpp,B,N0,W0,P0,l0,N1,W1,P1,l1)
//         cond (poly) is the polynomial responsible of the branch
//         B is the new discussed basis. (It can contain less polynomials when
//         some polynomial has been reduced to 0 by previous null-assumptions.
//         (N0,W0,P0,l0) and (N1,W1,P1,l1) are respectively the R-representation,
//         list of S-polys to be computed, and the last poly with assumed non-null
//         coefficient in both the null side and the non-null side.
static proc DiscussPolys(ideal B, ideal N, ideal W, list P, int l)
{
  list Pn=P; ideal Bn=B; int ln=l; ideal Nn=N; ideal Wn=W;
  int testsplit=0;
  poly f; poly lc; list L; int j;
  int l0; int l1; list P0; list P1; ideal N0; ideal W0; ideal N1; ideal W1;
  while((testsplit==0) and (ln<size(Bn)))
  {
    //f=redcoefs(Bn[ln+1],Nn);
    f=pnormalform(Bn[ln+1],Nn,Wn);
    if (f==0)
    {
      Bn=delfromideal(Bn,ln+1); //lppn=delfromideal(lppn,ln+1);
    }
    else
    {
      Bn[ln+1]=f;
      lc=leadcoef(f);
      L=splitR(Nn,Wn,lc);
      N0=L[1][1];
      W0=L[1][2];
      N1=L[2][1];
      W1=L[2][2];
      P1=addnewpairs(Bn,Pn,ln+1); // uses Buchberger pair selection and standard order
      if(N0[1]<>1)
      {
        testsplit=1;
        l0=ln; l1=ln+1;
        P0=Pn;
      }
      else
      {
        Pn=P1; P0=list(); ln=ln+1; Nn=N1; Wn=W1; l1=ln;
      }
    }
  }
  if(testsplit==0)
  {
    N1=Nn; W1=Wn; N0=ideal(1); W0=ideal(0); P0=list();
    l0=size(Bn); l1=size(Bn); P1=Pn;
  }
  return(list(lc,Bn,N0,W0,P0,l0,N1,W1,P1,l1));
}

// DiscussSPolys: given the data in a vertex of btcgs (BuildTree),
//               and when DiscussPolys has already built a vertex where
//               all the leadcoef are non-null in the R-representation,
//               it computes and reduces the S-polys in the list P in order
//               until it finds some non-reducing one. Then adds it to the
//               basis and modifies the list P.
//               Then it calls splitR and if the leadcoef non-null is, it
//               continues with the next S-poly in the list.
//               Else it finishes and recbtcgs will need to split.
//               15/09/2010
// ring @R
// input:
//    B:   (ideal) the actual basis
//    N:   (ideal) null conditions (R-rep)
//    W:   (ideal) non-null conditions set (R-rep)
//    P:   (list) of pairs of indices of S-polynomials that can and must be computed
//         (its leading coefficients are non-null, and using Buchberger's
//          criterions they are to be computed)
//    l:   (integer) representing the last polynomial in B for which the leading
//         coefficient is already assumed non-null.
// output: list of (cond,lpp,B,N0,W0,P0,l0,N1,W1,P1,l1)
//         cond (poly) is the polynomial responsible of the branch
//         B is the new discussed basis. (It can contain less polynomials when
//         some polynomial has been reduced to 0 by previous null-assumptions.
//         (N0,W0,P0,l0) and (N1,W1,P1,l1) are respectively the R-representation,
//         list of S-polys to be computed, and the last poly with assumed non-null
//         coefficient in both the null side and the non-null side.
static proc DiscussSPolys(ideal B,ideal N,ideal W,list P,int l)
{
  def RR=basering;
  list Pn=P; ideal Bn=B; int ln=l; ideal Nn=N; ideal Wn=W;
  int testsplit=0;
  poly lc; list L; int i; int j; poly S; list pair;
  int l0; int l1; list P0; list P1; ideal N0; ideal W0; ideal N1; ideal W1;
//  poly lc0;
  while((testsplit==0) and (size(Pn)<>0))
  {
    pair=Pn[1];
    i=pair[1]; j=pair[2];
    Pn=delete(Pn,1);
    lc=1; N1=Nn; W1=Wn;
    S=pspol(Bn[i],Bn[j]);
    S=pdivi(S,Bn)[1];
    //S=redcoefs(S,Nn);
    S=pnormalform(S,Nn,Wn);
    if (S<>0)
    {
      Bn[size(Bn)+1]=S;
      lc=leadcoef(S);
      ln=ln+1;
      L=splitR(Nn,Wn,lc);
      N0=L[1][1];
      W0=L[1][2];
      N1=L[2][1];
      W1=L[2][2];
      P1=addnewpairs(Bn,Pn,ln); // uses Buchberger pair selection and standard order
      if(N0[1]<>1)
      {
        testsplit=1;
        l0=ln-1; l1=ln;
        P0=Pn;
      }
      else
      {
        Pn=P1;  Nn=N1; Wn=W1; P0=list(); W0=ideal(0);
      }
    }
  }
  if(testsplit==0)
  {
    N0=ideal(1); W0=ideal(0); P0=list(); l0=0;  N1=Nn; W1=Wn;
    l1=size(Bn);
  }
  return(list(lc,Bn,N0,W0,P0,l0,N1,W1,P1,l1));
}

// cgsdr
// 20/09/2010
proc cgsdr(ideal F, list #)
"USAGE:   cgsdr(F); To compute a disjoint, reduced CGS.
          cgsdr is the starting point of the fundamental routine grobcov.
          It is to be used if only a disjoint reduced CGS is required.
          F: ideal in Q[a][x] (parameters and variables) to be discussed.

          Options: To modify the default options, pairs of arguments
          -option name, value- of valid options must be added to the call.

          Options:
            "null",ideal N: The default is "null",ideal(0).
            "nonnull",ideal W: The default "nonnull",ideal(1).
                When options "null" and/or "nonnull" are given, then
                the parameter space is restricted to V(N) \ V(h), where
                h is the product of the polynomials w in W.
            "comment",0-1: The default is "comment",0. Setting "comments",1
                will provide information about the development of the
                computation.
          One can give none till 3 of these options.
RETURN:   Returns a list T describing a reduced and disjoint comprehensive
          Groebner system (CGS), and whose segments correspond to
          constant leading power products (lpp) of the reduced Groebner
          basis. The returned list is of the form:
          (
            (lpp, (basis,segment),...,(basis,segment)),
            ..,,
            (lpp, (basis,segment),...,(basis,segment))
          )
          The bases are the reduced Groebner bases (after normalization)
          for each point of the corresponding segment.
          Each segment is given by a reduced representation (Ni,Wi), with
          Ni radical and V(Ni)=Zariski closure of the segment Si=V(Ni)\V(hi),
          where hi is the product of the polynomials w in Wi.
NOTE:     The basering R, must be of the form Q[a][x], a=parameters,
          x=variables, and should be defined previously, and the ideal
          defined on R.
KEYWORDS: CGS, disjoint, reduced, comprehensive Groebner system
EXAMPLE:  cgsdr; shows an example"
{
  list @T;
  exportto(Top,@T);
  setglobalrings();
  int i;
  ideal B;
  poly f;
  def N=ideal(0);
  def W=ideal(1);
  int comment=0;
  list L=#;
  for(i=1;i<=size(L)/2;i++)
  {
    if(L[2*i-1]=="null"){N=L[2*i];}
    else
    {
      if(L[2*i-1]=="nonnull"){W=L[2*i];}
      else
      {
        if(L[2*i-1]=="comment"){comment=L[2*i];}
      }
    }
  }
  if(N!=0)
  {
    def LL=Rrep(N,W);
    N=LL[1];
    W=LL[2];
    for (i=1;i<=size(F);i++)
    {
      f=pnormalform(F[i],N,W);
      if (f!=0){B[size(B)+1]=f;}
    }
  }
  else {B=F;}
  reccgsdr(B,N,W,list(),0);
  def T=@T;
  if (comment==1)
  {string("Number of segments in cgsdr (total) = ",size(T));}
  kill @T;
  kill @P; kill @RP; kill @R;
  return(grsegments(T));
}
example
{ "EXAMPLE:"; echo = 2;
  "Casas conjecture for degree 4";
  ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp;
  ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0),
          x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1),
          x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0),
          x2^2+(2*a3)*x2+(a2),
          x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0),
          x3+(a3);
  cgsdr(F);
}

//reccgsdr
// 20/09/2010
static proc reccgsdr(ideal B, ideal N, ideal W, list P, int l)
{
  ideal Bn=B; ideal Nn=N; ideal Wn=W; list Pn=P; int ln=l;  ideal lppn;
  list L; int i;
  poly lc; ideal N0; ideal W0; list P0; int l0;
  ideal N1=Nn; ideal W1=Wn; list P1=Pn; int l1=ln;
  if (l>0)
  {
    if (size(variables(B[l]))==0)
    {
      lppn=1; Bn=1;
      @T[size(@T)+1]=list(lppn,Bn,N,W);
      return();
    }
  }
  if (ln<size(Bn))
  {
    L=DiscussPolys(Bn, Nn, Wn, Pn, ln);
    lc=L[1]; Bn=L[2]; N0=L[3]; W0=L[4]; P0=L[5]; l0=L[6];
                                 N1=L[7]; W1=L[8]; P1=L[9]; l1=L[10];
    ln=l0;
  }
  if ((ln==size(Bn)) and (size(Bn)<>0))
  {
    L=DiscussSPolys(Bn, N1, W1, P1, l1);
    lc=L[1]; Bn=L[2]; N0=L[3]; W0=L[4]; P0=L[5]; l0=L[6];
                                 N1=L[7]; W1=L[8]; P1=L[9]; l1=L[10];
  }
  if (N0[1]<>1)
  {
    reccgsdr(Bn, N0,W0,P0,l0);
    reccgsdr(Bn, N1,W1,P1,l1);
  }
  else
  {
    if (equalideals(N1,ideal(1))==0)
    {
      Bn=mingb(Bn);
      Bn=redgb(Bn,N1,W1);
      lppn=ideal(0);
      for (i=1; i<=size(Bn);i++){lppn[i]=leadmonom(Bn[i]);}
      @T[size(@T)+1]=list(lppn,Bn,N1,W1);
    }
  }
}

// input:  internal routine called by cgsdr at the end to improve the output
// output: grouped segments by lpp obtained in cgsdr
static proc grsegments(list T)
{
  int i;
  list L;
  list lpp;
  list lp;
  list ls;
  int n=size(T);
  lpp[1]=T[n][1];
  L[1]=list(lpp[1],list(list(T[n][2],T[n][3],T[n][4])));
  if (n>1)
  {
    for (i=1;i<=size(T)-1;i++)
    {
      lp=memberpos(T[n-i][1],lpp);
      if(lp[1]==1)
      {
        ls=L[lp[2]][2];
        ls[size(ls)+1]=list(T[n-i][2],T[n-i][3],T[n-i][4]);
        L[lp[2]][2]=ls;
      }
      else
      {
        lpp[size(lpp)+1]=T[n-i][1];
        L[size(L)+1]=list(T[n-i][1],list(list(T[n-i][2],T[n-i][3],T[n-i][4])));
      }
    }
  }
  //"L in groupsegments="; L;
  return(L);
}

// grRtoPrep
// input:  L (list) is the output of cgsdr
// output: LL (list) the same list but the segments are expressed
//                   in canonical representations:
//  ( (lpp, (basis,
//             ((P_1),(P_{11},...,P_{1t1}))
//             ...
//             ((P_j),(P_{j1},...,P_{jtj}))
//          )
//          ...
//          (basis,
//             ((P_1),(P_{11},...,P_{1t1}))
//             ...
//             ((P_j),(P_{j1},...,P_{jtj}))
//          )
//    )
//    ...
//    (lpp, (basis,
//             ((P_1),(P_{11},...,P_{1t1}))
//             ...
//             ((P_j),(P_{j1},...,P_{jtj}))
//          )
//          ...
//          (basis,
//             ((P_1),(P_{11},...,P_{1t1}))
//             ...
//             ((P_j),(P_{j1},...,P_{jtj}))
//          )
//    )
//  )
static proc grRtoPrep(list L)
{
  int i; int j;
  list LL; list ct;
  // size(L)=number of lpp-segments
  for (i=1;i<=size(L);i++)
  {
    LL[i]=list();
    LL[i][1]=L[i][1];
    // L[i][1]=lpp
    LL[i][2]=list();
    for (j=1;j<=size(L[i][2]);j++)
    {
      ct=RtoPrepNew(L[i][2][j][2],L[i][2][j][3]); // ,L[i][2][j][5]
      LL[i][2][j]=list();
      LL[i][2][j][1]=L[i][2][j][1];
      // L[i][2][j][1]=label
      LL[i][2][j][2]=L[i][2][j][2];
      // L[i][2][j][2]=basis
      LL[i][2][j][3]=ct;
    }
  }
  return(LL);
}

// idcontains
// input: ideal p, ideal q
// output: 1 if p contains q,  0 otherwise
static proc idcontains(ideal p, ideal q)
{
  int t; int i;
  t=1; i=1;
  def RR=basering;
  setring @P;
  def P=imap(RR,p);
  def Q=imap(RR,q);
  attrib(P,"isSB",1);
  poly r;
  while ((t==1) and (i<=size(Q)))
  {
    r=reduce(Q[i],P);
    if (r!=0){t=0;}
    i++;
  }
  setring RR;
  return(t);
}

// selectminindeals
//   given a list of ideals returns the list of integers corresponding
//   to the minimal ideals in the list
// input: L (list of ideals)
// output: the list of integers corresponding to the minimal ideals in L
static proc selectminideals(list L)
{
  if (size(L)==0){return(L);}
  def RR=basering;
  setring @P;
  def Lp=imap(RR,L);
  int i; int j; int t; intvec notsel;
  list P;
  for (i=1;i<=size(Lp);i++)
  {
    if(memberpos(i,notsel)[1]==1)
    {
      i++;
      if(i>size(Lp)){break;}
    }
    t=1;
    j=1;
    while ((t==1) and (j<=size(Lp)))
    {
      if (i==j){j++;}
      if ((j<=size(Lp)) and (memberpos(j,notsel)[1]==0))
      {

        if (idcontains(Lp[i],Lp[j])==1)
        {
          notsel[size(notsel)+1]=i;
          t=0;
        }
      }
      j++;
    }
    if (t==1){P[size(P)+1]=i;}
  }
  setring(RR);
  return(P);
}

// LCUnion
// Given a list of the P-representations of locally closed segments
// for which we know that the union is also locally closed
// it returns the P-representation of its union
// input:  L list of segments in P-representation
//      ((p_j^i,(p_j1^i,...,p_jk_j^i | j=1..t_i)) | i=1..s )
//      where i represents a segment
// output: P-representation of the union
//       ((P_j,(P_j1,...,P_jk_j | j=1..t)))
static proc LCUnion(list LL)
{
  def RR=basering;
  setring(@P);
  def L=imap(RR,LL);
  int i; int j; int k; list H; list C; list T;
  list L0; list P0; list P; list Q0; list Q;
  for (i=1;i<=size(L);i++)
  {
    for (j=1;j<=size(L[i]);j++)
    {
      P0[size(P0)+1]=L[i][j][1];
      L0[size(L0)+1]=intvec(i,j);
    }
  }
  Q0=selectminideals(P0);
  for (i=1;i<=size(Q0);i++)
  {
    Q[i]=L0[Q0[i]];
    P[i]=L[Q[i][1]][Q[i][2]];
  }
  // P is the list of the maximal components of the union
  //   with the corresponding initial holes.
  // Q is the list of intvec positions in L of the first element of the P's
  //   Its elements give (num of segment, num of max component (=min ideal))
  for (k=1;k<=size(Q);k++)
  {
    H=P[k][2]; // holes of P[k][1]
    for (i=1;i<=size(L);i++)
    {
      if (i!=Q[k][1])
      {
        for (j=1;j<=size(L[i]);j++)
        {
          C[size(C)+1]=L[i][j];
        }
      }
    }
    T[size(T)+1]=list(Q[k],P[k][1],addpart(H,C));
  }
  setring(RR);
  def TT=imap(@P,T);
  return(TT);
}

// Called by LCUnion to modify the holes of a primepart of the union
// by the addition of the segments that do not correspond to that part
// Works on @P ring.
// Input:
//   H=(p_i1,..,p_is) the holes of a component to be transformed by the addition of
//        the segments C that do not correspond to that component
//   C=((q_1,(q_11,..,q_1l_1)),..,(q_k,(q_k1,..,q_kl_k)))
//        the list of segments to be added to the holes
static proc addpart(list H, list C)
{
  list Q; int i; int j; int k; int l; int t; int t1;
  Q=H; intvec notQ; list QQ; list addq;
  ideal q;
  i=1;
  while (i<=size(Q))
  {
    if (memberpos(i,notQ)[1]==0)
    {
      q=Q[i];
      t=1; j=1;
      while ((t==1) and (j<=size(C)))
      {
        if (equalideals(q,C[j][1])==1)
        {
          t=0;
          for (k=1;k<=size(C[j][2]);k++)
          {
            t1=1;
            //kill addq;
            //list addq;
            l=1;
            while((t1==1) and (l<=size(Q)))
            {
              if ((l!=i) and (memberpos(l,notQ)[1]==0))
              {
                if (idcontains(C[j][2][k],Q[l])==1)
                {
                  t1=0;
                }
              }
              l++;
            }
            if (t1==1)
            {
              addq[size(addq)+1]=C[j][2][k];
            }
          }
          if((size(notQ)==1) and (notQ[1]==0)){notQ[1]=i;}
          else {notQ[size(notQ)+1]=i;}
        }
        j++;
      }
      if (size(addq)>0)
      {
        for (k=1;k<=size(addq);k++)
        {
          Q[size(Q)+1]=addq[k];
        }
        kill addq;
        list addq;
      }
      //print("Q="); Q; print("notQ="); notQ;
    }
    i++;
  }
  for (i=1;i<=size(Q);i++)
  {
    if(memberpos(i,notQ)[1]==0)
    {
      QQ[size(QQ)+1]=Q[i];
    }
  }
  if (size(QQ)==0){QQ[1]=ideal(1);}
  return(addpartfine(QQ,C));
}

// Called by addpart to finish the modification of the holes of a primepart
// of the union by the addition of the segments that do not correspond to
// that part.
// Works on @P ring.
static proc addpartfine(list H, list C0)
{
  int i; int j; int k; int te; intvec notQ; int l; list sel; int used;
  intvec jtesC;
  if ((size(H)==1) and (equalideals(H[1],ideal(1)))){return(H);}
  if (size(C0)==0){return(H);}
  def RR=basering;
  setring(@P);
  list newQ; list nQ; list Q; list nQ1; list Q0;
  def Q1=imap(RR,H);
  //Q1=sortlistideals(Q1);
  def C=imap(RR,C0);
  while(equallistideals(Q0,Q1)==0)
  {
    Q0=Q1;
    i=0;
    Q=Q1;
    kill notQ; intvec notQ;
    while(i<size(Q))
    {
      i++;
      for(j=1;j<=size(C);j++)
      {
        te=idcontains(Q[i],C[j][1]);
        if(te==1)
        {
          for(k=1;k<=size(C[j][2]);k++)
          {
            if(idcontains(Q[i],C[j][2][k])==1)
            {
              te=0; break;
            }
          }
          if (te==1)
          {
            used++;
            if ((size(notQ)==1) and (notQ[1]==0)){notQ[1]=i;}
            else{notQ[size(notQ)+1]=i;}
            kill newQ; list newQ;
            for(k=1;k<=size(C[j][2]);k++)
            {
              nQ=minGTZ(Q[i]+C[j][2][k]);
              for(l=1;l<=size(nQ);l++)
              {
                option(redSB);
                nQ[l]=std(nQ[l]);
                newQ[size(newQ)+1]=nQ[l];
              }
            }
            sel=selectminideals(newQ);
            kill nQ1; list nQ1;
            for(l=1;l<=size(sel);l++)
            {
              nQ1[l]=newQ[sel[l]];
            }
            newQ=nQ1;
            for(l=1;l<=size(newQ);l++)
            {
              Q[size(Q)+1]=newQ[l];
            }
            break;
          }
        }
      }
    }
    kill Q1; list Q1;
    for(i=1;i<=size(Q);i++)
    {
      if(memberpos(i,notQ)[1]==0)
      {
        Q1[size(Q1)+1]=Q[i];
      }
    }
    sel=selectminideals(Q1);
    kill nQ1; list nQ1;
    for(l=1;l<=size(sel);l++)
    {
      nQ1[l]=Q1[sel[l]];
    }
    Q1=nQ1;
  }
  setring(RR);
  //if(used>0){string("addpartfine was ", used, " times used");}
  return(imap(@P,Q1));
}

//// specswell
//// used only in specswellonlpp (not used, can be deleted)
//// input:
////   given two corresponding polynomials g1 and g2 with the same lpp
////   g1 belonging to the basis in the segment N1,W1
////   g2 belonging to the basis in the segment N2,W2
//// output:
////   1 if g1 spezializes well to g2 on the whole (N2,W2) segment
////   0 if not
//proc specswell(poly g1, poly g2, ideal N2, ideal W2)
//{
//  poly S;
//  S=leadcoef(g2)*g1-leadcoef(g1)*g2;
//  def RR=basering;
//  setring(@RPt);
//  def SR=imap(RR,S);
//  def N2R=imap(RR,N2);
//  attrib(N2R,"isSB",1);
//  poly S2R=reduce(SR,N2R);
//  setring(RR);
//  def S2=imap(@RPt,S2R);
//  //if (S2==0)
//  //if (nonnull(leadcoef(g1),N2,W2)==1)
//  if ((S2==0) and (nonnull(leadcoef(g1),N2,W2)))
//  {return(1);}
//  else {return(0);}
//}
//
//// specswellonlpp
//// not used, can be deleted
//// input:
////   given a generic polynomial g with given lpp
////   and the list of tripets (p,N,W) of all the segments in
////   the same lpp-segment, where p is the correct image of g on (N,W)
//// output:
////   1 if g spezializes well to p on the whole (N,W) segment for all segments
////   0 if not
//proc specswellonlpp(poly g, list L)
//{
//  int i=1; int t=1;
//  while ((t==1) and (i<=size(L)))
//  {
//    t=specswell(g, L[i][1],L[i][2],L[i][3]);
//    i++;
//  }
//  return(t);
//}

// specswellCrep
// input:
//   given two corresponding polynomials g1 and g2 with the same lpp
//   g1 belonging to the basis in the segment ida1,idb1
//   g2 belonging to the basis in the segment ida2,idb2
// output:
//   1 if g1 spezializes well to g2 on the whole (ida2,idb2) segment
//   0 if not
static proc specswellCrep(poly g1, poly g2, ideal ida2)
{
  poly S;
  S=leadcoef(g2)*g1-leadcoef(g1)*g2;
  def RR=basering;
  setring(@RPt);
  def SR=imap(RR,S);
  def ida2R=imap(RR,ida2);
  attrib(ida2R,"isSB",1);
  poly S2R=reduce(SR,ida2R);
  setring(RR);
  def S2=imap(@RPt,S2R);
  if (S2==0){return(1);}   // and (nonnullCrep(leadcoef(g1),ida2,idb2))
  else {return(0);}
}


// gcover
// input: ideal F: a generating set of a homogeneous ideal in Q[a][x]
//    list GenCase: Containing the generic case with basis 1 if it exists
//    list #: optional
// output: the list
//   S=((lpp, generic basis, Rrep, Crep),..,(lpp, generic basis, Rrep, Crep))
//      where a Rrep is ( (p1,(p11,..,p1k_1)),..,(pj,(pj1,..,p1k_j)) )
//            a Crep is ( ida, idb )
static proc gcover(ideal F,list GenCase, list #)
{
  int i; int j; int k; ideal lpp; list GPi2; list pairspP; ideal B; int ti;
  int i1; int tes; int j1; int selind; int i2; int m;
  list prep; list crep; list LCU; poly p; poly lcp; list L; ideal FF;
  list NW=#;
  int CGS=NW[3];
  int comment=NW[4];
  NW=NW[1],NW[2];
  list GS; list GP;
  def RR=basering;
  int start=timer; int start0=start; int start1=start;
  if (CGS==0)
  {
    def BT=buildtree(F,list("null",NW[1],"nonnull",NW[2]));
    setglobalrings();
    def FC=finalcases(BT);
    GS=groupsegments(FC);
    if(comment==1)
    {
      string("Number of segments in buildtree (total) = ",size(FC));
      string("Number of lpp segments in groupsegments = ",size(GS));
      string("Time in buildtree = ",timer-start," sec");
    }
    start=timer;
    GP=groupRtoPrep(GS);
    if (comment==1){string("Time in groupRtoPrep = ",timer-start," sec");}
  }
  else
  {
    GS=cgsdr(F,list("null",NW[1],"nonnull",NW[2],"comment",comment));
    setglobalrings();
    if(comment==1)
    {
      string("Number of lpp segments in cgsdr = ",size(GS));
      string("Time in cgsdr = ",timer-start," sec");
    }
    start=timer;
    GP=grRtoPrep(GS);
    if(comment==1){string("Time in grRtoPrep = ",timer-start," sec");}
  }
  for(i=1;i<=size(GP);i++)
  {
    if(size(GP[i][2])>1){GP[i][3]=1;}
    else{GP[i][3]=0;}
  }
  int SizeGC=size(GenCase);
  if (SizeGC>0)
  {
    int te=0;
    list NewGen; list CH;
    for (i=1;i<=size(GP);i++)
    {
      if(equalideals(GP[i][1],ideal(1))==1)
      {
        te=1;
        NewGen[1]=GenCase;
        for(j=1;j<=size(GP[i][2]);j++)
        {
          NewGen[j+1]=GP[i][2][j];
        }
        GP[i][2]=NewGen;
        if(i!=1)
        { \\exchange cases i and 1
          CH=GP[i];
          GP[i]=GP[1];
          GP[1]=CH;
        }
        break;
      }
    }
    if (te==0) // add GenCase as a new case
    {
      CH[1]=GenCase;
      //CH[1]=list(ideal(1),list(GenCase));
      for (i=1;i<=size(GP);i++)
      {
        CH[i+1]=GP[i];
      }
      GP=CH;
    }
  }
  for(i=1;i<=size(GP);i++)
  {
    GP[i][3]=size(GP[i][2]);
  }
  list LL;
  list S;
  poly sp;
  ideal BB;
  start1=timer;
  for (i=1;i<=size(GP);i++)
  {
    kill LL;
    list LL;
    lpp=GP[i][1];
    GPi2=GP[i][2];
    kill pairspP; list pairspP;
    for(j=1;j<=size(GPi2);j++)
    {
      pairspP[size(pairspP)+1]=GPi2[j][3];
    }
    LCU=LCUnion(pairspP);
    kill prep; list prep;
    for(k=1;k<=size(LCU);k++)
    {
      prep[k]=list(LCU[k][2],LCU[k][3]);
      if (CGS==0)
      {
        B=GPi2[LCU[k][1][1]][2];
      }
      else
      {
        B=GPi2[LCU[k][1][1]][1];
      }
      LCU[k][1]=B;
    }
    // Deciding if combine is needed
    kill BB;
    ideal BB;
    tes=1; m=1;
    while((tes==1) and (m<=size(LCU[1][1])))
    {
      j=1;
      while((tes==1) and (j<=size(LCU)))
      {
        k=1;
        while((tes==1) and (k<=size(LCU)))
        {
          if(j!=k)
          {
            sp=pnormalform(pspol(LCU[j][1][m],LCU[k][1][m]),LCU[k][2],NW[2]);
            if(sp!=0){tes=0;}
          }
          k++;
        }
        if(tes==1)
        {
          BB[m]=LCU[j][1][m];
        }
        j++;
      }
      if(tes==0){break;}
      m++;
    }
    crep=PtoCrep(prep);
    if(tes==0)
    {
      // combine is needed
      kill B; ideal B;
      for (j=1;j<=size(LCU);j++)
      {
        LL[j]=LCU[j][2];
      }
      if (size(LCU)>1)
      {
        FF=precombint(LL);
      }
      for (k=1;k<=size(lpp);k++)
      {
        kill L; list L;
        for (j=1;j<=size(LCU);j++)
        {
          L[j]=list(LCU[j][2],LCU[j][1][k]);
        }
        if (size(LCU)>1)
        {
          B[k]=combine(L,FF);
        }
        else{B[k]=L[1][2];}
      }
    }
    else{B=BB;}
    for(j=1;j<=size(B);j++)
    {
      B[j]=pnormalform(B[j],crep[1],NW[2]);
    }
    S[i]=list(lpp,B,prep,crep,GP[i][3]);
  }
  if(comment==1)
  {
    string("Time in LCUnion + combine = ",timer-start1," sec");
  }
  kill @P; kill @RP; kill @R;
  return(S);
}

// grobcov
// input:
//    ideal F: a parametric ideal in Q[a][x], where a are the parameters
//             and x the variables
//    list #: (options) list("null",N,"nonnull",W,"can",Method,"cgs",CGS), where
//            N is the null conditions ideal (if desired)
//            W is the ideal of non-null conditions (if desired)
//            Method is 1 by default and can be set to 0 if we do not
//            need to obtain the canonical GC, but only a GC.
//            CGS is 1 by default and uses cgsdr. It can be set to 0 to
//            use the old buildtree instead.
// output:
//    list S: ((lpp,basis,(idp_1,(idp_11,..,idp_1s_1))), ..
//             (lpp,basis,(idp_r,(idp_r1,..,idp_rs_r))) ) where
//            each element of S corresponds to a lpp-segment
//            given by the lpp, the basis, and the P-representation of the segment
proc grobcov(ideal F,list #)
"USAGE:   grobcov(F); This is the fundamental routine of the
          library. It computes the Groebner cover of a parametric ideal
          (see (*) Montes A., Wibmer M., Groebner Bases for Polynomial
          Systems with parameters. JSC 45 (2010) 1391-1425.)
          The Groebner cover of a parametric ideal consist of a set of pairs
          (S_i,B_i), where the S_i are disjoint locally closed segments
          of the parameter space, and the B_i are the reduced Groebner
          bases of the ideal on every point of S_i.

          The ideal F must be defined on a parametric ring Q[a][x].
          Options: To modify the default options, pair of arguments
          -option name, value- of valid options must be added to the call.

          Options:
            "null",ideal N: The default is "null",ideal(0).
            "nonnull",ideal W: The default "nonnull",ideal(1).
                When options "null" and/or "nonnull" are given, then
                the parameter space is restricted to V(N) \ V(h), where
                h is the product of the polynomials w in W.
            "can",0-1: The default is "can",1. With the default option
                the homogenized ideal is computed before obtaining the
                Groebner cover, so that the result is the canonical
                Groebner cover. Setting "can",0 only homogenizes the basis
                so the result is not exactly canonical, but the computation
                is more efficient.
            "ext",0-1: The default is "ext",1. With the default option the
                full representation of the bases is computed (possible
                shaves) and often a simpler result is obtained. Setting
                "ext",0 only the generic representation is computed
                (single polynomials, but not specializing to non-zero at
                each point of the segment.
            "cgs",0-1: The default is "cgs",1. The default option uses the
                cgsdr routine of the actual library to compute the initial
                CGS (more efficient). Setting "cgs",0 it uses the routine
                cgsdrold of the old library redcgs.lib. This option can be
                tested if the default option does not terminate.
            "comment",0-1: The default is "comment",0. Setting "comments",1
                will provide information about the development of the
                computation.
          One can give none till 6 of these options.
RETURN:   The list
          (
           (lpp_1,basis_1,P-representation_1),
           ...
           (lpp_s,basis_s,P-represntation_s)
          )

          The lpp are constant over a segment and correspond to the
          set of lpp of the reduced Groebner basis for each point
          of the segment.

          Basis: to each element of lpp corresponds an I-regular function given          Groebner basis, and it is given in full representation (by
          in full representation (by default option "ext",1) or in
          generic representation (option "ext",0). The regular function is
          the corresponding element of the reduced Groebner basis for
          each point of the segment with the given lpp.
          For each point in the segment, the polynomial or the set of
          polynomials representing it, if they do not specialize to 0,
          then after normalization, specialize to the corresponding
          element of the reduced Groebner basis.

          The P-representation of a segment is of the form
          ((p_1,(p_11,..,p_1k1)),..,(p_r,(p_r1,..,p_rkr))
          representing the segment U_i (V(p_i) \ U_j (V(p_ij))), where the
          p's are prime ideals.

NOTE:     The basering R, must be of the form Q[a][x], a=parameters,
          x=variables, and should be defined previously. The ideal must
          be defined on R.
KEYWORDS: Groebner cover, parametric ideal, canonical, discussion of
          parametric ideal, multigrobcov, gencase1.
EXAMPLE:  grobcov; shows an example"
{
  list S; int i; int ish=1; list GBR; list BR; int j; int k;
  list NW; ideal idp; ideal idq; int s; ideal ext; list SS;
  ideal N; ideal W; int canop;  int extop; int CGS; int repop;
  int gradorder; int comment=0; int m;
  list L=#;
  // default options
  int start=timer;
  def RR=basering;
  list NW0;
  W=ideal(1);
  N=ideal(0);
  canop=1; // canop=0 for homogenizing the basis but not the ideal (not canonical)
           // canop=1 for working with the homogenized ideal
  repop=0; // repop=0 for representing the segments in Prep
           // repop=1 for representing the segments in Crep
           // repop=2 for representing the segments in Prep and Crep
  extop=1; // extop=1 if the full representation of the bases are to be computed
           // extop=0 if only generic representation of the bases are to be computed
  CGS=1;   // CGS=1 if cgsdr is to be used (default)
           // CGS=0 if buildtree is to be used instead
  for(i=1;i<=size(L)/2;i++)
  {
    if(L[2*i-1]=="can"){canop=L[2*i];}
    else
    {
      if(L[2*i-1]=="ext"){extop=L[2*i];}
      else
      {
        if(L[2*i-1]=="rep"){repop=L[2*i];}
        else
        {
          if(L[2*i-1]=="null"){N=L[2*i];}
          else
          {
            if(L[2*i-1]=="nonnull"){W=L[2*i];}
            else
            {
              if (L[2*i-1]=="cgs"){CGS=L[2*i];}
              else
              {
                if (L[2*i-1]=="comment"){comment=L[2*i];}
              }
            }
          }
        }
      }
    }
  }
  if (comment==1){string("Options: can = ",canop,", extend = ",extop,", cgs = ",CGS,", rep = ",repop);}
  for (i=1;i<=size(F);i++){ish=ishomog(F[i]); if(ish==0){break;}}
  NW0=list(N,W,CGS,comment);
  if (ish==1)
  {
    kill S;
    list gc;
    def S=gcover(F,gc,NW0);
    setglobalrings();
  }
  else
  {
    list RRL=ringlist(RR);
    if (RRL[3][1][1]=="dp"){gradorder=1;} else {gradorder=0;}
    RRL[3][1][1]="dp";
    //RRL[1][3][1][1]="dp"; // COMMENTED GIVES ERROR IN S53.
    def Pa=ring(RRL[1]);
    list Lx;
    Lx[1]=0;
    Lx[2]=RRL[2]+RRL[1][2];
    Lx[3]=RRL[1][3];
    Lx[4]=RRL[1][4];
    RRL[1]=0;
    def D=ring(RRL);
    def RP=D+Pa;
    setring(RP);
    def F1=imap(RR,F);
    def NW1=imap(RR,NW0);
    int gcyes=0;
    if (canop==1)
    {
      option(redSB);
      def F11=std(F1);
      setring(RR);
      list gc;
      def F2=imap(RP,F11);
      def NW2=imap(RP,NW1);
      if (size(NW2[1])==0)
      {
        gc=gencase1(F2,"compbas",0);
        if (size(gc)>0)
        {
          gcyes=1;
          NW2[1]=gc[4];
          //gc=delete(gc,4);
          list gcn;
          gcn[1]=ideal(1); // lpp
          gcn[2]=list(list(ideal(1),ideal(0),list(gc[3])));
          gc=gcn;
        }
      }
    }
    else
    {
      setring(RR);
      def NW2=NW0;
      def F2=imap(RP,F1);
    }
    //setglobalrings();
    setring RR; // ja hi es ?
    RRL=ringlist(RR);
    //if (RRL[3][1][1]!="dp"){ERROR("the order must be dp");}
    poly @t;
    ring H=0,@t,dp;
    def RH=RR+H;
    setring(RH);
    //kill @P;
    //kill @RP;
    //kill @R;
    //setglobalrings();
    //setring(@Rt);
    def FH=imap(RR,F2);
    list gcH;
    if (gcyes==1)
    {
      gcH=imap(RR,gc);
    }
    def NWH=imap(RR,NW2);
    for (i=1;i<=size(FH);i++)
    {
      FH[i]=homog(FH[i],@t);
    }
    def G=gcover(FH,gcH,NWH); // list(NWH[1],NWH[2],CGS,comment));
    for (i=1;i<=size(G);i++)
    {
      G[i][1]=subst(G[i][1],@t,1);
      G[i][2]=subst(G[i][2],@t,1);
    }
    setring(RR);
    setglobalrings();
    S=imap(RH,G);
    for (i=1;i<=size(S);i++)
    {
      S[i][2]=postredgb(mingb(S[i][2]));
      S[i][1]=postredgb(mingb(S[i][1]));
    }
  }
  // Now Extend;
  poly leadc;
  if (extop==1)
  {
    int start1=timer;
    for (i=1;i<=size(S);i++)
    {
      m=size(S[i][2]);
      for (j=1;j<=size(S[i][2]);j++)
      {
        idp=S[i][4][1];
        idq=S[i][4][2];
        if (size(idp)>0)
        {
          leadc=leadcoef(S[i][2][j]);
          kill ext;
          def ext=extend(S[i][2][j],idp,idq);
          if (typeof(ext)=="poly")
          {
            S[i][2][j]=pnormalform(ext,idp,W);
            //"T_Polynomial after extend="; S[i][2][j];
          }
          else
          {
            if(size(ext)==1)
            {
              S[i][2][j]=ext[1];
            }
            else
            {
              kill SS; list SS;
              for(s=1;s<=size(ext);s++)
              {
                ext[s]=pnormalform(ext[s],idp,W);
              }
              for(s=1;s<=size(S[i][2]);s++)
              {
                if(s!=j){SS[s]=S[i][2][s];}
                else{SS[s]=ext;}
              }
              S[i][2]=SS;
            }
          }
          //"T_ poly or ideal after extend="; S[i][2][j];
        }
      }
    }
    if(comment==1){string("Time in extend = ",timer-start1," sec");}
  }
  list Si; list nS;
  if (repop==0)
  {
    for(i=1;i<=size(S);i++)
    {
      Si=list(S[i][1],S[i][2],S[i][3]);
      nS[size(nS)+1]=Si;
    }
    S=nS;
  }
  else
  {
    if (repop==1)
    {
      for(i=1;i<=size(S);i++)
      {
        Si=list(S[i][1],S[i][2],S[i][4]);
        nS[size(nS)+1]=Si;
      }
      S=nS;
    }
  }
  kill @P; kill @RP; kill @R;
  if (comment==1)
  {
    string("Time for grobcov = ", timer-start," sec");
    string("Number of segments of grobcov = ", size(S));
  }
  return(S);
}
example
{ "EXAMPLE:"; echo = 2;
  "Casas conjecture for degree 4";
  ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp;
  ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0),
          x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1),
          x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0),
          x2^2+(2*a3)*x2+(a2),
          x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0),
          x3+(a3);
  grobcov(F);
}


// input:
//    poly g in K[a],
//    list P=(p_1,..p_r) representing a minimal prime decomposition
// output
//    poly f such taht f notin p_i forall i and
//           g-f in p_i forall i such that g notin p_i
static proc nonzerodivisor(poly gr, list Pr)
{
  def RR=basering;
  setring(@P);
  def g=imap(RR,gr);
  def P=imap(RR,Pr);
  int i; int k;  list J; ideal F;
  def f=g;
  ideal Pi;
  for (i=1;i<=size(P);i++)
  {
    option(redSB);
    Pi=std(P[i]);
    //attrib(Pi,"isST",1);
    if (reduce(g,Pi,1)==0){J[size(J)+1]=i;}
  }
  for (i=1;i<=size(J);i++)
  {
    F=ideal(1);
    for (k=1;k<=size(P);k++)
    {
      if (k!=J[i])
      {
        F=idint(F,P[k]);
      }
    }
    f=f+F[1];
  }
  setring(RR);
  def fr=imap(@P,f);
  return(fr);
}

// input:
//   int i:
//   list LPr: (p1,..,pr) of prime components of an ideal in K[a]
// output:
//   list (fr,fnr) of two polynomials that are equal on V(pi)
//       and fr=0 on V(P) \ V(pi), and fnr is nonzero on V(pj) for all j.
static proc deltai(int i, list LPr)
{
  def RR=basering;
  setring(@P);
  def LP=imap(RR,LPr);
  int j; poly p;
  def F=ideal(1);
  poly f;
  poly fn;
  ideal LPi;
  for (j=1;j<=size(LP);j++)
  {
    if (j!=i)
    {
      F=idint(F,LP[j]);
    }
  }
  p=0; j=1;
  while ((p==0) and (j<=size(F)))
  {
    LPi=LP[i];
    attrib(LPi,"isSB",1);
    p=reduce(F[j],LPi);
    j++;
  }
  f=F[j-1];
  fn=nonzerodivisor(f,LP);
  setring(RR);
  def fr=imap(@P,f);
  def fnr=imap(@P,fn);
  return(list(fr,fnr));
}

// input: a list of pairs ((p1,P1),..,(pr,Pr)) where
//    ideal pi is a prime component
//    poly Pi is the polynomial in K[a][x] on V(pi)\ V(Mi)
//    (p1,..,pr) are the prime decomposition of the lpp-segment
//    list crep =(ideal ida,ideal idb): the Crep of the segment.
//    list Pci of the intersecctions of all pj except the ith one
// output:
//    poly P on an open and dense set of V(p_1 int ... p_r)
static proc combine(list L, ideal F)
{
  // ATTENTION REVISE AND USE Pci and F
  int i; poly f;
  f=0;
  for(i=1;i<=size(L);i++)
  {
    f=f+F[i]*L[i][2];
  }
  f=elimconstfac(f);
  return(f);
}

// elimconstfac: eliminate the factors in the polynom f that are in K[a]
// input:
//   poly f:
//   list L: of components of the segment
// output:
//   poly f2  where the factors of f in K[a] that are non-null on any component
//   have been dropped from f
static proc elimconstfac(poly f)
{
  int cond; int i; int j; int t;
  if (f==0){return(f);}
  def RR=basering;
  setring(@R);
  poly ff=imap(RR,f);
  list l=factorize(ff,0);
  poly f1=1;
  for(i=2;i<=size(l[1]);i++)
  {
      f1=f1*(l[1][i])^(l[2][i]);
  }
  setring(RR);
  def f2=imap(@R,f1);
  return(f2);
}

// input:
//   poly f:  a polynomial in K[a]
//   ideal P: an ideal in K[a]
//   called from ring @R
// output:
//   t:  with value 1 if f reduces modulo P, 0 if not.
static proc nullin(poly f,ideal P)
{
  int t;
  def RR=basering;
  setring(@P);
  poly f0=imap(RR,f);
  ideal P0=imap(RR,P);
  attrib(P0,"isSB",1);
  if (reduce(f0,P0,1)==0){t=1;}
  else{t=0;}
  setring(RR);
  return(t);
}

static proc polyinparamsonly(poly f)
{
  int t;
  def RR=basering;
  setring @R;
  def f0=imap(RR,f);
  if (size(variables(f0))==0){t=1;}
  else{t=0;}
  setring(RR);
  return(t);
}

// monoms
static proc monoms(poly f)
{
  list L;
  if (f!=0) { L[size(f)]=list();}
  poly lm; poly lc; poly lp; poly Q; poly mQ;
  def p=f;
  int i=1;
  while (p!=0)
  {
    lm=lead(p);
    p=p-lm;
    lc=leadcoef(lm);
    lp=leadmonom(lm);
    L[i]=list(lc,lp);
    i++;
  }
  return(L);
}

// input:
//   poly f: a generic polynomial in the basis
//   ideal idp: such that ideal(S)=idp
//   ideal idq: such that S=V(idp)\V(idq)
////   NW the list of ((N1,W1),..,(Ns,Ws)) of red-rep of the grouped
////      segments in the lpp-segment  NO MORE USED
// output:
static proc extend(poly f, ideal idp, ideal idq)
{
  matrix CC; poly Q; list NewMonoms;
  int i;  int j;  poly fout; ideal idout;
  list L=monoms(f);
  int nummonoms=size(L)-1;
  Q=L[1][1];
  if (nummonoms==0){return(f);}
  for (i=2;i<=size(L);i++)
  {
    CC=matrix(extendcoef(L[i][1],Q,idp,idq));
    NewMonoms[i-1]=list(CC,L[i][2]);
  }
  if (nummonoms==1)
  {
    for(j=1;j<=ncols(NewMonoms[1][1]);j++)
    {
      fout=NewMonoms[1][1][2,j]*L[1][2]+NewMonoms[1][1][1,j]*NewMonoms[1][2];
      //fout=pnormalform(fout,idp,W);
      if(ncols(NewMonoms[1][1])>1){idout[j]=fout;}
    }
    if(ncols(NewMonoms[1][1])==1){return(fout);} else{return(idout);}
  }
  else
  {
    //int start=timer;
    list cfi;
    list coefs;
    for (i=1;i<=nummonoms;i++)
    {
      kill cfi; list cfi;
      for(j=1;j<=ncols(NewMonoms[i][1]);j++)
      {
        cfi[size(cfi)+1]=NewMonoms[i][1][2,j];
      }
      coefs[i]=cfi;
    }
    def indexpolys=findindexpolys(coefs);
    for(i=1;i<=size(indexpolys);i++)
    {
      fout=L[1][2];
      for(j=1;j<=nummonoms;j++)
      {
        fout=fout+(NewMonoms[j][1][1,indexpolys[i][j]])/(NewMonoms[j][1][2,indexpolys[i][j]])*NewMonoms[j][2];
      }
      fout=cleardenom(fout);
      if(size(indexpolys)>1){idout[i]=fout;}
    }
    if (size(indexpolys)==1){return(fout);} else{return(idout);}
  }
}

// input:
//   list coefs=( (q11,..,q1r_1),..,(qs1,..,qsr_1) )
//               of denominators of the monoms
// output:
//   list ind=(v_1,..,v_t) of intvec
//        each intvec v=(i_1,..,is) corresponds to a polynomial in the sheaf
//        that will be built from it in extend procedure.
static proc findindexpolys(list coefs)
{
  int i; int j; intvec numdens;
  for(i=1;i<=size(coefs);i++)
  {
    numdens[i]=size(coefs[i]);
  }
  def RR=basering;
  setring(@P);
  def coefsp=imap(RR,coefs);
  ideal cof; list combpolys; intvec v; int te; list mp;
  for(i=1;i<=size(coefsp);i++)
  {
    cof=ideal(0);
    for(j=1;j<=size(coefsp[i]);j++)
    {
      cof[j]=factorize(coefsp[i][j],3);
    }
    coefsp[i]=cof;
  }
  for(j=1;j<=size(coefsp[1]);j++)
  {
    v[1]=j;
    te=1;
    for (i=2;i<=size(coefsp);i++)
    {
      mp=memberpos(coefsp[1][j],coefsp[i]);
      if(mp[1])
      {
        v[i]=mp[2];
      }
      else{v[i]=0;}
    }
    combpolys[j]=v;
  }
  combpolys=reform(combpolys,numdens);
  setring RR;
  return(combpolys);
}


// extendcoef: given Q,P in K[a] where P/Q specializes on an open and dense subset
//      of the whole V(p1 int...int pr), it returns a basis of the module
//      of all syzygies equivalent to P/Q,
static proc extendcoef(poly P, poly Q, ideal idp, ideal idq)
{
  def RR=basering;
  setring(@P);
  def PL=ringlist(@P);
  PL[3][1][1]="dp";
  def P1=ring(PL);
  setring(P1);
  ideal idp0=imap(RR,idp);
  option(redSB);
  qring q=std(idp0);
  poly P0=imap(RR,P);
  poly Q0=imap(RR,Q);
  ideal PQ=Q0,-P0;
  module C=syz(PQ);
  setring @P;
  def idp1=imap(RR,idp);
  def idq1=imap(RR,idq);
  def C1=matrix(imap(q,C));
  def redC=selectregularfun(C1,idp1,idq1);
  setring(RR);
  def CC=imap(@P,redC);
  return(CC);
}

// input:
//   list L of the polynomials matrix CC
//      (we assume that one of them is non-null on V(N)\V(M))
//   ideal N, ideal M: ideals representing the locally closed set V(N)\V(M)
// assume to work in @P
static proc selectregularfun(matrix CC, ideal NN, ideal MM)
{
  int numcombused;
  def RR=basering;
  setring @P;
  def C=imap(RR,CC);
  def N=imap(RR,NN);
  def M=imap(RR,MM);
  if (ncols(C)==1){return(C);}

  int i; int j; int k; list c; intvec ci; intvec c0; intvec c1;
  list T; list T0; list T1; list LL; ideal N1;ideal M1; int te=0;
  for(i=1;i<=ncols(C);i++)
  {
    if((C[1,i]!=0) and (C[2,i]!=0))
    {
      if(c0==intvec(0)){c0[1]=i;}
      else{c0[size(c0)+1]=i;}
    }
  }
  def C1=submat(C,1..2,c0);
  for (i=1;i<=ncols(C1);i++)
  {
    c=comb(ncols(C1),i);
    for(j=1;j<=size(c);j++)
    {
      ci=c[j];
      numcombused++;
      if(i==1){N1=N+C1[2,j]; M1=M;}
      if(i>1)
      {
        kill c0; intvec c0 ; kill c1; intvec c1;
        c1=ci[size(ci)];
        for(k=1;k<size(ci);k++){c0[k]=ci[k];}
        T0=searchinlist(c0,LL);
        T1=searchinlist(c1,LL);
        N1=T0[1]+T1[1];
        M1=intersect(T0[2],T1[2]);
      }
      T=list(ci,PtoCrep(Prep(N1,M1)));
      LL[size(LL)+1]=T;
      if(equalideals(T[2][1],ideal(1))==1){te=1; break;}
    }
    if(te==1){break;}
  }
  ci=T[1];
  def Cs=submat(C1,1..2,ci);
  setring RR;
  return(imap(@P,Cs));
}

// input:
//   intvec c:
//   list L=( (c1,T1),..(ck,Tk) )
//      where the c's are assumed to be intvects
// output:
//   object T with index c
static proc searchinlist(intvec c,list L)
{
  int i; list T;
  for(i=1;i<=size(L);i++)
  {
    if (L[i][1]==c)
    {
      T=L[i][2];
      break;
    }
  }
  return(T);
}

// Input: C0 the matrtix of (P1,..,Pr)
//                          (Q1,..,Qr) of the regular function of a coefficient (P,Q)
//        NW0 the list of ((N1,W1),..(Ns,Ws)) of red-rep of the grouped
//        segments in the lpp-segment
// Output: (B, T) where
//        B is the submatrix of the selected minimal representants for the
//        regular function
//        T the matrix of ones and zeroes whose colums are associated
//        to the colums of B, with 1 in the segments where the representant
//        is nonnull and 0 if it can be.
static proc redext(matrix C0, list NW0)
{
  def RR=basering;
  setring(@P);
  def C=imap(RR,C0);
  def NW=imap(RR,NW0);
  int nc=ncols(C);
  int nr=size(NW);
  intmat T[nr][nc];
  int i; int j; int k; int t;
  for (i=1;i<=nc;i++)
  {
    for (j=1;j<=nr;j++)
    {
      t=nonnull(C[i][2],NW[j][1],NW[j][2]); // (Q,N,W)
      T[j,i]=t;
    }
  }
  int h; int tt=0;
  intvec c; intvec r;
  list cc;  int l;
  for (j=1;j<=2;j++){r[j]=j;}
  i=1;
  while((i<=nc) and (tt==0))
  {
    cc=comb(nc,i);
    tt=0;
    l=1;
    while((tt==0) and (l<=size(cc)))
    {
      tt=1;
      c=cc[l];
      j=1;
      while ((j<=nr) and (tt==1))
      {
        h=0;
        k=1;
        while ((h==0) and (k<=i))
        {
          if(T[j,c[k]]==1){h=1;}
          k++;
        }
        if (h==0){tt=0;}
        j++;
      }
      l++;
    }
    i++;
  }
  if (tt==0){"extendcoef does not extend to the whole S";}
  intvec rr;
  for (i=1;i<=nr;i++){rr[i]=i;}
  def B=submat(C,r,c);
  def TT=submat(T,rr,c);
  setring(RR);
  return(list(imap(@P,B),imap(@P,TT)));
}

// comb: the list of combinations of elements (1,..n) of order p
static proc comb(int n, int p)
{
  list L; list L0;
  intvec c; intvec d;
  int i; int j; int last;
  if ((n<0) or (n<p))
  {
    return(L);
  }
  if (p==1)
  {
    for (i=1;i<=n;i++)
    {
      c=i;
      L[size(L)+1]=c;
    }
    return(L);
  }
  else
  {
    L0=comb(n,p-1);
    for (i=1;i<=size(L0);i++)
    {
      c=L0[i]; d=c;
      last=c[size(c)];
      for (j=last+1;j<=n;j++)
      {
        d[size(c)+1]=j;
        L[size(L)+1]=d;
      }
    }
    return(L);
  }
}

// selectminsheaves
// Input: L=((v_11,..,v_1k_1),..,(v_s1,..,v_sk_s))
//    where:
//    The s lists correspond to the s coefficients of the polynomial f
//    (v_i1,..,v_ik_i) correspond to the k_i intvec v_ij of the
//    spezializations of the jth rekpresentant (Q,P) of the ith coefficient
//    v_ij is an intvec of size equal to the number of little segments
//    forming the lpp-segment of 0,1, where 1 represents that it specializes
//    to non-zedro an the whole little segment and 0 if not.
// Output: S=(w_1,..,w_j)
//    where the w_l=(n_l1,..,n_ls) are intvec of length size(L), where
//    n_lt fixes which element of (v_t1,..,v_tk_t) is to be
//    choosen to form the tth (Q,P) for the lth element of the sheaf
//    representing the I-regular function.
// The selection is done to obtian the minimal number of elements
//    of the sheaf that specializes to non-null everywhere.
static proc selectminsheaves(list L)
{
  list C=allsheaves(L);
  return(smsheaves(C[1],C[2]));
}

// Input:
//   list C of all the combrep
//   list L of the intvec that correesponds to each element of C
// Output:
//   list LL of the subsets of C that cover all the subsegments
//   (the union of the corresponding L(C) has all 1).
static proc smsheaves(list C, list L)
{
  int i; int i0; intvec W;
  int nor; int norn;
  intvec p;
  int sp=size(L[1]); int j0=1;
  for (i=1;i<=sp;i++){p[i]=1;}
  while (p!=0)
  {
    i0=0; nor=0;
    for (i=1; i<=size(L); i++)
    {
      norn=numones(L[i],pos(p));
      if (nor<norn){nor=norn; i0=i;}
    }
    W[j0]=i0;
    j0++;
    p=actualize(p,L[i0]);
  }
  list LL;
  for (i=1;i<=size(W);i++)
  {
    LL[size(LL)+1]=C[W[i]];
  }
  return(LL);
}

// allsheaves
// Input: L=((v_11,..,v_1k_1),..,(v_s1,..,v_sk_s))
//    where:
//    The s lists correspond to the s coefficients of the polynomial f
//    (v_i1,..,v_ik_i) correspond to the k_i intvec v_ij of the
//    spezializations of the jth rekpresentant (Q,P) of the ith coefficient
//    v_ij is an intvec of size equal to the number of little segments
//    forming the lpp-segment of 0,1, where 1 represents that it specializes
//    to non-zero on the whole little segment and 1 if not.
// Output:
//    (list LL, list LLS)  where
//    LL is the list of all combrep
//    LLS is the list of intvec of the corresponding elements of LL
static proc allsheaves(list L)
{
  intvec V; list LL; intvec W; int r; intvec U;
  int i; int j; int k;
  int s=size(L[1][1]); // s = number of little segments of the lpp-segment
  list LLS;
  for (i=1;i<=size(L);i++)
  {
    V[i]=size(L[i]);
  }
  LL=combrep(V);
  for (i=1;i<=size(LL);i++)
  {
    W=LL[i];   // size(W)= number of coefficients of the polynomial
    kill U; intvec U;
    for (j=1;j<=s;j++)
    {
      k=1; r=1; U[j]=1;
      while((r==1) and (k<=size(W)))
      {
        if(L[k][W[k]][j]==0){r=0; U[j]=0;}
        k++;
      }
    }
    LLS[i]=U;
  }
  return(list(LL,LLS));
}

// numones
// Input:
//   intvec v of (0,1) in each position
//   intvec pos: the positions to test
// Output:
//   int nor: the nuber of 1 of v in the positions given by pos.
static proc numones(intvec v, intvec pos)
{
  int i; int n;
  for (i=1;i<=size(pos);i++)
  {
    if (v[pos[i]]==1){n++;}
  }
  return(n);
}

// Input:  intvec p of zeros and ones
// Output: intvec W of the positions where p has ones.
static proc pos(intvec p)
{
  int i;
  intvec W; int j=1;
  for (i=1; i<=size(p); i++)
  {
    if (p[i]==1){W[j]=i; j++;}
  }
  return(W);
}

// actualize: actualizes zeroes of p
// Input:
//   intvec p: of zeroes and ones
//   intvec c: of zeroes and ones (of the same length)
// Output;
//   intvec pp: of zeroes and ones, where a 0 stays in pp[i] if either
//   already p[i]==0 or c[i]==1.
static proc actualize(intvec p, intvec c)
{
  int i; intvec pp=p;
  for (i=1;i<=size(p);i++)
  {
    if ((pp[i]==1) and (c[i]==1)){pp[i]=0;}
  }
  return(pp);
}

// combrep
// Input: V=(n_1,..,n_i)
// Output: L=(v_1,..,v_p) where p=prod_j=1^i (n_j)
//    is the list of all intvec v_j=(v_j1,..,v_ji) where 1<=v_jk<=n_i
static proc combrep(intvec V)
{
  list L; list LL;
  int i; int j; int k;  intvec W;
  if (size(V)==1)
  {
    for (i=1;i<=V[1];i++)
    {
      L[i]=intvec(i);
    }
    return(L);
  }
  for (i=1;i<size(V);i++)
  {
    W[i]=V[i];
  }
  LL=combrep(W);
  for (i=1;i<=size(LL);i++)
  {
    W=LL[i];
    for (j=1;j<=V[size(V)];j++)
    {
      W[size(V)]=j;
      L[size(L)+1]=W;
    }
  }
  return(L);
}

static proc reducemodN(poly f,ideal N)
{
  def RR=basering;
  setring(@RPt);
  def fa=imap(RR,f);
  def Na=imap(RR,N);
  attrib(Na,"isSB",1);
  // //option(redSB);
  // Na=std(Na);
  fa=reduce(fa,Na);
  setring(RR);
  def f1=imap(@RPt,fa);
  return(f1);
}

// computes the intersection of the ideals in S in @P
static proc intersp(list S)
{
  def RR=basering;
  setring(@P);
  def SP=imap(RR,S);
  option(returnSB);
  def NP=intersect(SP[1..size(SP)]);
  setring(RR);
  return(imap(@P,NP));
}

static proc radicalmember(poly f,ideal ida)
{
  int te;
  def RR=basering;
  setring(@P);
  poly fp=imap(RR,f);
  ideal idap=imap(RR,ida);
  poly @t;
  ring H=0,@t,dp;
  def PH=@P+H;
  setring(PH);
  def fH=imap(@P,fp);
  ideal idaH=imap(@P,idap);
  idaH[ncols(idaH)+1]=1-@t*fH;
  option(redSB);
  ideal G=std(idaH);
  //"G="; G;
  if (G==1){te=1;} else {te=0;}
  setring(RR);
  return(te);
}

// returns 1 if the poly f is nonnull on V(N)\V(M), 0 otherwise.
static proc NonNull(poly f, ideal N, ideal M)
{
  int te=1; int i;
  def RR=basering;
  setring(@P);
  poly fp=imap(RR,f);
  ideal Np=imap(RR,N);
  ideal Mp=imap(RR,M);
  ideal H;
  ideal Nf=Np+fp;
  for (i=1;i<=ncols(Mp);i++)
  {
    te=radicalmember(Mp[i],Nf);
    if (te==0) break;
  }
  setring RR;
  return(te);
}

// input:
//    matrix CC: CC=(p_a1 .. p_ar_a)
//                  (q_a1 .. q_ar_a)
//            the matrix of elements of a coefficient in oo[a].
//    (ideal ida, ideal idb): the canonical representation of the segment S.
// output:
//    list caout
//            the minimum set of elements of CC needed such that at least one
//            of the q's is non-null on S, as well as the C-rep of of the
//            points where the q's are null on S.
//            The elements of caout are of the form (p,q,prep);
static proc selectextendcoef(matrix CC, ideal ida, ideal idb)
{
  def RR=basering;
  setring(@P);
  def ca=imap(RR,CC);
  def N0=imap(RR,ida);
  ideal N;
  def M=imap(RR,idb);
  int r=ncols(ca);
  int i; int te=1; list com; int j; int k; intvec c; list prep;
  list cs; list caout;
  i=1;
  while ((i<=r) and (te==1))
  {
    com=comb(r,i);
    j=1;
    while((j<=size(com)) and (te==1))
    {
      N=N0;
      c=com[j];
      for (k=1;k<=i;k++)
      {
        N=N+ca[2,c[k]];
      }
      prep=Prep(N,M);
      if (i==1)
      {
        cs[j]=list(ca[1,j],ca[2,j],prep);
      }
      if ((size(prep)==1) and (equalideals(prep[1][1],ideal(1))))
      {
        te=0;
        for(k=1;k<=size(c);k++)
        {
          caout[k]=cs[c[k]];
        }
      }
      j++;
    }
    i++;
  }
  if (te==1){"error: extendcoef does not extend to the whole S";}
  setring(RR);
  return(imap(@P,caout));
}

// input:
//   ideal N1: in some basering (depends only on the parameters)
//   ideal N2: in some basering (depends only on the parameters)
// output:
//   ideal Np=N1+N2; computed in P
static proc plusP(ideal N1,ideal N2)
{
  def RR=basering;
  setring(@P);
  def N1p=imap(RR,N1);
  def N2p=imap(RR,N2);
  def Np=N1p+N2p;
  setring RR;
  return(imap(@P,Np));
}

// input:
//   list combpolys: (v1,..,vs)
//      where vi are intvec.
//   output outcomb: (w1,..,wt)
//      whre wi are intvec.
//      All the vi without zeroes are in outcomb, and those with zeroes are
//         combined to form new intvec with the rest
static proc reform(list combpolys, intvec numdens)
{
  list combp0; list combp1; int i; int j; int k; int l; list rest; intvec notfree;
  list free; intvec free1; int te; intvec v;  intvec w;
  int nummonoms=size(combpolys[1]);
  for(i=1;i<=size(combpolys);i++)
  {
    if(memberpos(0,combpolys[i])[1]==1)
    {
      combp0[size(combp0)+1]=combpolys[i];
    }
    else {combp1[size(combp1)+1]=combpolys[i];}
  }
  for(i=1;i<=nummonoms;i++)
  {
    kill notfree; intvec notfree;
    for(j=1;j<=size(combpolys);j++)
    {
      if(combpolys[j][i]<>0)
      {
        if(notfree[1]==0){notfree[1]=combpolys[j][i];}
        else{notfree[size(notfree)+1]=combpolys[j][i];}
      }
    }
    kill free1; intvec free1;
    for(j=1;j<=numdens[i];j++)
    {
      if(memberpos(j,notfree)[1]==0)
      {
        if(free1[1]==0){free1[1]=j;}
        else{free1[size(free1)+1]=j;}
      }
      free[i]=free1;
    }
  }
  list amplcombp; list aux;
  for(i=1;i<=size(combp0);i++)
  {
    v=combp0[i];
    kill amplcombp; list amplcombp;
    amplcombp[1]=intvec(v[1]);
    for(j=2;j<=size(v);j++)
    {
      if(v[j]!=0)
      {
        for(k=1;k<=size(amplcombp);k++)
        {
          w=amplcombp[k];
          w[size(w)+1]=v[j];
          amplcombp[k]=w;
        }
      }
      else
      {
        kill aux; list aux;
        for(k=1;k<=size(amplcombp);k++)
        {
          for(l=1;l<=size(free[j]);l++)
          {
            w=amplcombp[k];
            w[size(w)+1]=free[j][l];
            aux[size(aux)+1]=w;
          }
        }
        amplcombp=aux;
      }
    }
    for(j=1;j<=size(amplcombp);j++)
    {
      combp1[size(combp1)+1]=amplcombp[j];
    }
  }
  return(combp1);
}

static proc nonnullCrep(poly f0,ideal ida0,ideal idb0)
{
  int i;
  def RR=basering;
  setring(@P);
  def f=imap(RR,f0);
  def ida=imap(RR,ida0);
  def idb=imap(RR,idb0);
  def idaf=ida+f;
  int te=1;
  for(i=1;i<=size(idb);i++)
  {
    if(radicalmember(idb[i],idaf)==0)
    {
      te=0; break;
    }
  }
  setring(RR);
  return(te);
}

// input:  L: list of ideals (works in @P)
// output: F0: ideal of polys. F0[i] is a poly in the intersection of
//             all ideals in L except in the ith one, where it is not.
//             L=(p1,..,ps);  F0=(f1,..,fs);
//             F0[i] \in intersect_{j#i} p_i
static proc precombint(list L)
{
  int i; int j; int tes;
  def RR=basering;
  setring(@P);
  list L0; list L1; list L2; list L3; ideal F;
  L0=imap(RR,L);
  L1[1]=L0[1]; L2[1]=L0[size(L0)];
  for (i=2;i<=size(L0)-1;i++)
  {
    L1[i]=intersect(L1[i-1],L0[i]);
    L2[i]=intersect(L2[i-1],L0[size(L0)-i+1]);
  }
  L3[1]=L2[size(L2)];
  for (i=2;i<=size(L0)-1;i++)
  {
    L3[i]=intersect(L1[i-1],L2[size(L0)-i]);
  }
  L3[size(L0)]=L1[size(L1)];
  for (i=1;i<=size(L3);i++)
  {
    option(redSB); L3[i]=std(L3[i]);
  }
  for (i=1;i<=size(L3);i++)
  {
    tes=1; j=0;
    while((tes==1) and (j<size(L3[i])))
    {
      j++;
      option(redSB);
      L0[i]=std(L0[i]);
      if(reduce(L3[i][j],L0[i])!=0){tes=0; F[i]=L3[i][j];}
    }
    if (tes==1){"ERROR a polynomial in all p_j except p_i was not found";}
  }
  setring(RR);
  def F0=imap(@P,F);
  return(F0);
}

// precombinediscussion
// not used, can be deleted
// input:  list L: the LCU segment with bases for each pi component
// output: intvec vv:  vv[1]=(1 if the generic polynomial of the vv[2]
//                     component already specializes well,
//                     0 if combine is to be used)
//                     vv[2]=selind, the index for which the generic basis
//                     already specializes well if combine is not to be used (vv[1]=1).
static proc precombinediscussion(L,crep)
{
  int tes=1; int selind; int i1; int j1; poly p; poly lcp; intvec vv;
  if (size(L)==1){vv=1,1; return(vv);}
  for (i1=1;i1<=size(L);i1++)
  {
    tes=1;
    p=L[i1][2];
    lcp=leadcoef(p);


    if(nonnullCrep(lcp,crep[1],crep[2])==1)
    {
      for(j1=1;j1<=size(L);j1++)
      {
        if(i1!=j1)
        {
          if(specswellCrep(p,L[j1][2],L[j1][1])==0){tes=0; break;}
        }
      }
    }
    else{tes=0;}
    if(tes==1){selind=i1; break;}
  }
  vv=tes,selind;
  return(vv);
}

// only if N=0 and W=1
proc gencase1(ideal F, list #)
"USAGE:   gencase1(F); This routine determines the generic segment when
          the generic case has basis 1, and returns the empty list if not.
          It is useful, for example in automatic discovery of geometric
          theorems, to determine the prime varieties over which solutions exist.
          It can work, even if the complete grobcov does not finish.
          It serves to obtain a partial result that can be sometimes very useful.
          It is also used internally in the canonical computation grobcov,
          but can be called by the user. Only the basering Q[a][x] needs
          to be defined and the ideal given in this ring.
          Options: It allows an option list("compbas",0-1),
          If the routine is called with option
          ("compbas",0), then the given ideal must be the reduced
          Groebner basis of the ideal in the ring Q[x,a].
          If the routine is called by the user this option not to be used,
          and the algorithm will compute internally the reduced Groebner
          basis of the ideal in the ring Q[x,a].
RETURN:   The list of the generic case, when its basis is 1, or
          the empty list if not.
          The output is of the form
          (lpp=1,basis=1,(null ideal=0,(p1,..ps)),N)
          where (0,(p1,..,ps)) is the P-representation of the generic segment
          (the pi's are the prime components) and N is its intersection
NOTE:     The basering R, must be of the form Q[a][x], a=parameters,
          x=variables, and should be defined previously. The ideal must
          be defined on R.
KEYWORDS: generic segment, automatic discovery of geometric theorems,
EXAMPLE:  gencase1; shows an example"
{
  int compbas=1; list L=#;
     // compbas==1 the gbasis wrt vars+param must be computed now
     // compbas==0 the gbasis wrt vars+param is already computed
  def RR=basering; list empty; int i;
  setglobalrings();
  for(i=1;i<=size(L)/2;i++)
  {
    if(L[2*i-1]=="compbas"){compbas=L[2*i];}
  }
  if (compbas==1)
  {
    setring(@RP);
    def FP=imap(R,F);
    option(redSB);
    def G=std(FP);
    setring(RR);
    def F1=imap(@RP,G);
  }
  else {def F1=F;}
  ideal Zero;
  for(i=1;i<=size(F1);i++)
  {
    if (leadmonom(F1[i])==1)
    {
      Zero[size(Zero)+1]=F1[i];
    }
  }
  if (size(Zero)>0)
  {
    setring(@P);
    def ZeroP=imap(RR,Zero);
    //def N=radical(ZeroP);
    def holes=minGTZ(ZeroP);
    for(i=1;i<=size(holes);i++)
    {
      option(redSB);
      holes[i]=std(holes[i]);
    }
    def N=holes[1];
    for(i=2;i<=size(holes);i++)
    {
      N=intersect(N,holes[i]);
    }
    option(redSB);
    N=std(N);
    setring(RR);
    def hole=imap(@P,holes);
    def Nn=imap(@P,N);
    kill @P; kill @RP; kill @R;
    return(ideal(1),ideal(1),list(ideal(0),hole),Nn);
  }
  else
  {
    kill @P; kill @RP; kill @R;
    setring(RR);
    return(empty);
  }
}
example
{ "EXAMPLE:"; echo = 2;
  "Generic segment for the extended Steiner-Lehmus theorem";
  ring R=(0,x,y),(a,b,m,n,p,r),lp;
  ideal S=p^2-(x^2+y^2),
          -a*(y)+b*(x+p),
          -a*y+b*(x-1)+y,
          (r-1)^2-((x-1)^2+y^2),
          -m*(y)+n*(x+r-2) +y,
          -m*y+n*x,
          (a^2+b^2)-((m-1)^2+n^2);
  short=0;
  gencase1(S);
}

// minAssGTZ eliminating denominators
static proc minGTZ(ideal N);
{
  int i; int j;
  def L=minAssGTZ(N);
  for(i=1;i<=size(L);i++)
  {
    for(j=1;j<=size(L[i]);j++)
    {
      L[i][j]=cleardenom(L[i][j]);
    }
  }
  return(L);
}

proc multigrobcov(ideal F, list #)
"USAGE:   multigrobcov(F); This routine is to be used instead of grobcov
          when grobcov does not finish, and the generic case is expected
          to have basis 1. It can be useful for automating discovery of
          geometric theorems.
          The ideal F must be defined on a parametric ring Q[a][x].
          If the generic basis is not 1, then it returns the empty list,
          but if the generic basis is one then it computes the
          grobcov over each irreducible component of the complement of
          the generic segment and returns the generic segment and the
          different grobcov on each segment. From the result, the global
          grobcov can be deduced eliminating convenablement the inter-
          sections of the different grobcov computed over the components.
          Options: A list of options of the form
          ("comment",0-1,"can",0-1 can,"cgs",0-1,"ext",0-1), can be given.
          One can give none till 4 of these options by giving the
          name of the option and the value. Options "null" and "nonnull" are
          avoided.
          When option ("comment",1) is set, the routine provides information
          about the development of the computation. The default option
          is ("comment",0).
          When option ("can",0) is given, then the computation is
          done homogenizing the given basis but not computing the
          whole homogenized ideal. Thus in this case the result is not
          completely canonical but it is also useful. This option
          facilitates the computation. The default option is ("can",1).
          When option ("cgs",0) is set, then instead of using cgsdr
          for computing the initial reduced disjoint CGS, then
          cgsdrold is used. This can be tested when ("cgs",1) (the default
          option) fails. When option ("ext",0) is set, only the generic
          representation of the bases are computed instead of the
          full representation (the default option is ("ext",1)).
RETURN:   The list whose first element is the generic case, and the
          remaining elements are the grobcov over the different irreducible
          components in the complementary of the generic segment.
          the empty list if the generic case does not have basis 1.
NOTE:     The basering R, must be of the form Q[a][x], a=parameters,
          x=variables, and should be defined previously. The ideal must
          be defined on R.
KEYWORDS: grobcov, generic segment, automatic discovery of geometric theorems,
EXAMPLE:  multigrobcov; shows an example."
{
  int i; int comment=1; list L=#; ideal N; list gc; list GC; list GCA;
  int start=timer; int ni; int nw;
  for(i=1;i<=size(L)/2;i++)
  {
    if (L[2*i-1]=="comment"){comment=L[2*i];}
    else
    {
      if(L[2*i-1]=="null"){"multigrobcov does not allow null restriction"; ni=i;}
      else
      {
        if(L[2*i-1]=="nonnull"){"multigrobcov does not allow nonnull restriction"; nw=i;}
      }
    }
  }
  if (ni>0)
  {
    L=delete(L,2*ni-1); L=delete(L,2*ni-1);
    if(nw>0)
    {
      if(nw<ni)
      {
        L=delete(L,2*nw-1); L=delete(L,2*nw-1);
      }
      else
      {
        L=delete(L,2*nw-3); L=delete(L,2*nw-3);
      }
    }
  }
  else
  {
    if (nw>0){L=delete(L,2*nw-1);L=delete(L,2*nw-1);}
  }
  gc=gencase1(F);
  if(size(gc)==0)
  {
    string("The generic case is not 1, thus multigrobcov is not useful");
    return(gc);
  }
  else
  {
    if(comment==1){"Generic case ="; gc; " ";}
    def SS2=gc[3][2];
    GCA=list(list(list(gc[1],gc[2],list(gc[3]))));
    if(comment==1){"Components to study="; SS2;}
    for (i=1;i<=size(SS2);i++)
    {
      N=SS2[i];
      if(comment==1){" "; "Begin grobcov on the variety N ="; N;}
      L[size(L)+1]="null"; L[size(L)+1]=N;
      //"T_L=";L;
      GC=grobcov(F,L);
      GCA[size(GCA)+1]=GC;
    }
    if(comment==1){string("Time for multigrobcov = ",timer-start);}
    return(GCA);
  }
}
example
{
  "Generalization of the Steiner-Lehmus theorem";
  ring R=(0,x,y),(a,b,m,n,p,r),lp;
  ideal S=p^2-(x^2+y^2),
          -a*(y)+b*(x+p),
          -a*y+b*(x-1)+y,
          (r-1)^2-((x-1)^2+y^2),
          -m*(y)+n*(x+r-2) +y,
          -m*y+n*x,
          (a^2+b^2)-((m-1)^2+n^2);
  short=0;
  multigrobcov(S,list("can",0,"cgs",0,"comment",1));
}

proc cgsdrold(ideal F, list #)
"USAGE:   cgsdrold(F); To compute a disjoint, reduced CGS.
          From the old library redcgs.lib.
          cgsdrold is the starting point of the fundamental routine
          grobcovold of the library redcgs.lib.
          Use instead cgsdr. cgsdrold is only recommended for comparison
          with cgsdr or for didactic purposes to plot the tree (buildtree)
          using the routine buildtreetoMaple.
          F: ideal in Q[a][x] (parameters and variables) to be discussed.

          Options: To modify the default options, pairs of arguments
          -option name, value- of valid options must be added to the call.

          Options:
            "null",ideal N: The default is "null",ideal(0).
            "nonnull",ideal W: The default "nonnull",ideal(1).
                When options "null" and/or "nonnull" are given, then
                the parameter space is restricted to V(N) \ V(h), where
                h is the product of the polynomials w in W.
            "old",0-1: The default option is "old",1 that gives an output
                analogous to the one obtained by cgsdr. Setting "old",0
                provides an output representing a tree (buildtree), that
                can be plotted using the routine buildtreetoMaple.
            "comment",0-1: The default is "comment",0. Setting "comments",1
                will provide information about the development of the
                computation.
          One can give none till 4 of these options.
RETURN:   With the default option "old",1, it returns a list T describing
          a reduced and disjoint comprehensive Groebner system (CGS),
          whose segments correspond to constant leading power products (lpp)
          of the reduced Groebner basis. The returned list is of the form:
          (
            (lpp, (basis,segment),...,(basis,segment)),
            ..,,
            (lpp, (basis,segment),...,(basis,segment))
          )
          The bases are the reduced Groebner bases (after normalization)
          for each point of the corresponding segment.
          Each segment is given by a reduced representation (Ni,Wi), with
          Ni radical and V(Ni)=Zariski closure of the segment Si=V(Ni)\V(hi),
          where hi is the product of the polynomials w in Wi.
          Setting option "old",0 the output represents the tree and
          can then be transformed to a plot structure using the routine
          buildtreetoMaple.
          Its structure in this case is:
          The first element of the list is the root, and contains
            [1] label: intvec(-1)
            [2] number of children : int
            [3] the ideal F
            [4], [5], [6] the red-representation of the segment
                (null, non-null conditions, prime components of the null
                conditions) given (as option).
                ideal (0), ideal (1), list(ideal(0)) is assumed if
                no optional conditions are given.
            [7] the set of lpp of ideal F
            [8] condition that was taken to reach the vertex
                (poly 1, for the root).
          The remaining elements of the list represent vertices of the tree:
          with the same structure:
            [1] label: intvec (1,0,0,1,...) gives its position in the tree:
                first branch condition is taken non-null, second null,...
            [2] number of children (0 if it is a terminal vertex)
            [3] the specialized ideal with the previous assumed conditions
                to reach the vertex
            [4],[5],[6] the red-representation of the segment corresponding
                to the previous assumed conditions to reach the vertex
            [7] the set of lpp of the specialized ideal at this stage
            [8] condition that was taken to reach the vertex from the
                father's vertex (that was taken non-null if the last
                integer in the label is 1, and null if it is 0)
          The terminal vertices form a disjoint partition of the parameter
          space whose bases specialize to the reduced Groebner basis of the
          specialized ideal on each point of the segment and preserve
          the lpp. They form a disjoint reduced CGS, and is the only
          vertices grouped and ordered by lpp that is returned with the
          default option "old",1.

NOTE:     The basering R, must be of the form Q[a][x], a=parameters,
          x=variables, and should be defined previously, and the ideal
          defined on R.
KEYWORDS: CGS, cgsdr, buildtree, buildtreetoMaple, disjoint, reduced,
          comprehensive Groebner system
EXAMPLE:  cgsdrold; shows an example"
{
  int i; list L=#; int oldop=1;
  for(i=1;i<=size(L)/2;i++)
  {
    if(L[2*i-1]=="old"){oldop=L[2*i];}
  }
  def bt=buildtree(F, #);
  if (oldop==0){return(bt);}
  else
  {
    setglobalrings();
    def gs=groupsegments(finalcases(bt));
    int j;
    for (i=1;i<=size(gs);i++)
    {
      for (j=1;j<=size(gs[i][2]);j++)
      {
        gs[i][2][j]=delete(gs[i][2][j],1);
        gs[i][2][j]=delete(gs[i][2][j],4);
        if(equalideals(gs[i][2][j][3],ideal(0))){gs[i][2][j][3]=ideal(1);}
      }
    }
    kill @P; kill @R; kill @RP;
    return(gs);
  }
}
example
{ "EXAMPLE:"; echo = 2;
  ring R=(0,a1,a2,a3,a4),(x1,x2,x3,x4),dp;
  ideal F=x4-a4+a2,
          x1+x2+x3+x4-a1-a3-a4,
          x1*x3*x4-a1*a3*a4,
          x1*x3+x1*x4+x2*x3+x3*x4-a1*a4-a1*a3-a3*a4;
  cgsdrold(F);
  cgsdrold(F,"old",0);
}

proc grobcovold(ideal F,list #)
"USAGE:   grobcovold(F); This is the fundamental routine of the
          old library redcgs.lib. It is somewhat heuristic and does
          not certify the obtention of the canonical Groebner cover of
          a parametric ideal, as does grobcov, but usually it does or
          provides a warning if not. It allows different options, recalling
          all the different approaches of the old library redcgs.lib.
          Use grobcov instead. The use of grobcovold is only recommended
          to compare results or study alternatives.

          The ideal F must be defined on a parametric ring Q[a][x].
          Options: To modify the default options, pair of arguments
          -option name, value- of valid options must be added to the call.

          Options:
            "null",ideal N: The default is "null",ideal(0).
            "nonnull",ideal W: The default "nonnull",ideal(1).
                When options "null" and/or "nonnull" are given, then
                the parameter space is restricted to V(N) \ V(h), where
                h is the product of the polynomials w in W.
            "can",0-2: The default is "can",1. With the default option
                the homogenized ideal is computed before obtaining the
                Groebner cover, so that the result is the canonical
                Groebner cover. Setting "can",0 only homogenizes the basis
                so the result is not exactly canonical, but the computation
                is more efficient. Setting "can",2 no homogenization of
                the ideal is carried out, and the segments with same lpp
                are added so much as possible when a common basis is obtained.
                The result, in this case is not canonical nor the segments
                are always locally closed. Nevertheless it can have
                less segments as the canonical result.
            "out",0-1: The default is "out",0. With the default option the
                output is analogous to that of grobcov. If option "can",2
                is also set, then this representation can be somewhat
                confusing, because the segments are not always given in
                P-representation, as they are not always locally closed.
                With option "out",1 a representation in tree form is given
                providing a canonical representation of the segments, even if
                they are not locally closed. This representation can be transformed
                by the routine cantreetoMaple into a file that can be read
                in Maple and plotted with the plotcantree Maple routine of
                the old dpgb library, showing the tree.
            "comment",0-1: The default is "comment",0. Setting "comments",1
                will provide information about the development of the
                computation.
          One can give none till 5 of these options.
RETURN:   With the default option ("out",0), the list
          (
           (lpp_1,basis_1,P-representation_1)
           ...
           (lpp_s,basis_s,P-represntation_s)
          )
          With option "out",1, a list T representing a rooted tree.
          Each element of the list T has the two first entries with the
          following content:
           [1]: The label (intvec) representing the position in the rooted
                tree:  0 for the root (and this is a special element)
                       i for the root of the segment i
                       (i,...) for the children of the segment i
           [2]: the number of children (int) of the vertex.
          There are three kind of vertices:
           (1) the root (first element labelled 0),
           (2) the vertices labelled with a single integer i,
           (3) the rest of vertices labelled with more indices.
          Description of the root. Vertex type (1)
           There is a special vertex (the first one) whose content is
           the following:
             [3] lpp of the given ideal
             [4] the given ideal
             [5] the R-representation  of the (optional) given null and
                 non-null conditions.
             [6] CRCGS, RCGS, MRCGS depending on the "can" option (1,0,2).
           Description of vertices type (2). These are the vertices that
           initiate a segment, and are labelled with a single integer.
             [3] lpp (ideal) of the reduced basis. If they are repeated lpp's this
                 will correspond to a sheaf.
             [4] the reduced basis (ideal) of the segment.
           Description of vertices type (3). These vertices have as first
           label i and descend form vertex i in the position of the label
           (i,...). They contain moreover a unique prime ideal in the parameters
           and form ascending chains of ideals.
          How is to be read the mrcgs tree? The vertices with an even number of
          integers in the label are to be considered as additive and those
          with an odd number of integers in the label are to be considered as
          substraction. As an example consider the following vertices:
          v1=((i),2,lpp,B),
          v2=((i,1),2,P_(i,1)),
          v3=((i,1,1),2,P_(i,1,1)),
          v4=((i,1,1,1),1,P_(i,1,1,1)),
          v5=((i,1,1,1,1),0,P_(i,1,1,1,1)),
          v6=((i,1,1,2),1,P_(i,1,1,2)),
          v7=((i,1,1,2,1),0,P_(i,1,1,2,1)),
          v8=((i,1,2),0,P_(i,1,2)),
          v9=((i,2),1,P_(i,2)),
          v10=((i,2,1),0,P_(i,2,1)),
          They represent the segment:
          (V(i,1)\(((V(i,1,1) \ ((V(i,1,1,1) \ V(i,1,1,1,1)) u (V(i,1,1,2) \ V(i,1,1,2,1)))))
          u V(i,1,2))) u (V(i,2) \ V(i,2,1))
          and can also be represented by
          (V(i,1) \ (V(i,1,1) u V(i,1,2))) u
          (V(i,1,1,1) \ V(i,1,1,1)) u
          (V(i,1,1,2) \ V(i,1,1,2,1)) u
          (V(i,2) \ V(i,2,1))
          where V(i,j,..) = V(P_(i,j,..))

          The lpp are constant over a segment and correspond to the
          set of lpp of the reduced Groebner basis for each point
          of the segment.

          Basis: to each element of lpp corresponds an I-regular function given          Groebner basis, and it is given in full representation (by
          in full representation. The regular function is
          the corresponding element of the reduced Groebner basis for
          each point of the segment with the given lpp.
          For each point in the segment, the polynomial or the set of
          polynomials representing it, if they do not specialize to 0,
          then after normalization, specialize to the corresponding
          element of the reduced Groebner basis.

          The P-representation of a segment is of the form
          ((p_1,(p_11,..,p_1k1)),..,(p_r,(p_r1,..,p_rkr))
          representing the segment U_i (V(p_i) \ U_j (V(p_ij))), where the
          p's are prime ideals.

NOTE:     The basering R, must be of the form Q[a][x], a=parameters,
          x=variables, and should be defined previously. The ideal must
          be defined on R.
KEYWORDS: Groebner cover, grobcov, parametric ideal, canonical, discussion of
          parametric ideal.
EXAMPLE:  grobcovold; shows an example"
{
  int i;
  list LL=#;
  list T; list NT; list NTe;
  // default options
  int comment=0; int canop=1; int outop=0;
  int start=timer;
  ideal W=ideal(1);
  ideal N=ideal(0);
  canop=1; // canop=0 for homogenizing the basis but not the ideal (not canonical)
           //         (old rcgs)
           // canop=1 for homogenizing the ideal
           //         (old crcgs)
           // canop=2 for not homogenizing and try to minimize the segments
           //         (old mrcgs)
  outop=0; // outop=0 for an output analogous to grobcov (if canop<>2)
           // outop=1 for an output as in the old library redcgs.lib
           //         in form of tree that can be transformed into Maple.
  for(i=1;i<=size(LL)/2;i++)
  {
    if(LL[2*i-1]=="can"){canop=LL[2*i];}
    else
    {
      if(LL[2*i-1]=="out"){outop=LL[2*i];}
      else
      {
        if (LL[2*i-1]=="comment"){comment=LL[2*i];}
      }
    }
  }
  if (comment>=1){string("can = ",canop," out = ", outop," comment = ",comment);}
  if (canop==0){T=rcgs(F,LL);}
  else
  {
    if (canop==1){T=crcgs(F,LL);}
    else
    {
      if (canop==2){T=mrcgs(F,LL);}
    }
  }
  if (comment>=1){string("Time in grobcovold = ",timer-start," sec");}
  if (outop==0)
  {
    // transforming the output to the modern form
    i=2; list Cap; int indCap; list Cua; ideal idp; list idq;
    int tes;
    while(i<=size(T))
    {
      kill Cap; list Cap;
      if(size(T[i][1])==1)
      {
        Cap=list(T[i][3],T[i][4]);
        indCap=T[i][1][1];
        i++;
      }
      kill Cua; list Cua;
      while(T[i][1][1]==indCap)
      {
        if(size(T[i][1]) mod 2 ==0)
        {
          if(size(idq)!=0){Cua[size(Cua)+1]=list(idp,idq);}
          kill idq; list idq;
          idp=T[i][3];
        }
        else
        {
          idq[size(idq)+1]=T[i][3];
        }
        i++;
        if(i>size(T)){break;}
      }
      Cua[size(Cua)+1]=list(idp,idq);
      Cap[3]=Cua;
      NT[size(NT)+1]=Cap;
      kill idp; ideal idp; kill idq; list idq;
    }
    if (comment==2){"rcgs="; T;}
    return(NT);
  }
  else
  {
    return(T);
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  "Simple robot: A. Montes,";
  "New algorithm for discussing Groebner bases with parameters,";
  "JSC, 33: 183-208 (2002).";
  ring R=(0,r,z,l),(s1,c1,s2,c2), dp;
  ideal S10=c1^2+s1^2-1,
            c2^2+s2^2-1,
            r-c1-l*c1*c2+l*s1*s2,
            z-s1-l*c1*s2-l*s1*c2;
  grobcovold(S10,"comment",1);
  grobcovold(S10,"can",2,"comment",1);
}