//////////////////////////////////////////////////////////////////////////////
version="$Id$";
category="General purpose";
info="
LIBRARY:  grobcov.lib   Groebner Cover for parametric ideals.
PURPOSE:  Comprehensive Groebner Systems, Groebner Cover, Canonical Forms,
          Parametric Polynomial Systems.
          The library contains Montes-Wibmer'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.
          For this purpose, in this library, the implemented algorithm is
          Kapur-Sun-Wang algorithm, because it is the most efficient
          algorithm known for this purpose.

          D. Kapur, Y. Sun, and D.K. Wang.
          \"A New Algorithm for Computing Comprehensive Groebner Systems\".
          Proceedings of ISSAC'2010, ACM Press, (2010), 29-36.

          cgsdr can be called directly if only a disjoint reduced
          Comprehensive Groebner System (CGS) is required.

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,
@*         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 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 routine than buildtree
                   (the own routine that is no more used). It uses
                   now KSW algorithm.

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 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 killed by the call to any of the
                   basic routines.

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

pnormalf(f,E,N);   Reduces a parametric polynomial f over V(E) \ V(N)
                   E is the null ideal and N the non-null ideal
                   over the parameters.

extend(GC); When the grobcov of an ideal has been computed
                   with the default option ('ext',0) and the explicit
                   option ('rep',2) (which is not the default), then
                   one can call extend (GC) (and options) to obtain the
                   full representation of the bases. With the default
                   option ('ext',0) only the generic representation of
                   the bases are computed, and one can obtain the full
                   representation using extend.

locus2d:      Special routine for determining the locus of points
                   of a two dimensional object. Given an ideal J with
                   two parameters (a,b) and so many variables as
                   needed, representing the system determining
                   the locus of points (a,b) who verify certain
                   geometrical properties, computing the grobcov of
                   J and applying to it locus2d, determines the locus.

locus2dto:   Transforms the output of locus2d to a string that
                   can be reed from different computational systems.

SEE ALSO: compregb_lib
";

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

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

// Library grobcov.lib
// (Groebner cover):
// Release 1: (public)
// Initial data: 21-1-2008
// Final data: 3-7-2008
// Release 2: (private)
// Initial data: 6-9-2009
// Final data: 25-10-2011
// Release 3: (this release, public)
// Initial data: 1-7-2012
// Final data: 4-9-2012
// 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, extend, pdivi, pnormalf, locus2d, locus2dto,
          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))
  {
    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;
ringlist(R);
ringlist(@P);
ringlist(@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)
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){setring(RR); 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);
};

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);
//}

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
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
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 a parametric ideal F in Q[a][x].
// input:
//   poly  f (in the parametric ring @R)
//   ideal F (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
          coefficient 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 te=0;
  if (defined(@P)==1){te=1;}
  else{setglobalrings();}
  def R=basering;
  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;
  if(te==0){kill @P; kill @R; kill @RP;}
  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 then S reduces by Buchberger 1st criterion
//                (not used)
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.
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){setring(RR); 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);
//}

// Ered: eliminates the factors in the polynom f that are non-null.
//       In ring @R
// input:
//   poly f:
//   ideal E  of null-conditions
//   ideal N  of non-null conditions
//        (E,N) represents V(E)\V(N),
//        Ered eliminates the non-null factors of f in V(E)\V(N)
// output:
//   poly f2  where the non-null conditions have been dropped from f
proc Ered(poly f,ideal E, ideal N)
{
  def RR=basering;
  setring(@R);
  poly ff=imap(RR,f);
  ideal EE=imap(RR,E);
  ideal NN=imap(RR,N);
  if((ff==0) or (equalideals(NN,ideal(1)))){setring(RR); return(f);}
  def v=variables(ff);
  int i;
  poly X=1;
  for(i=1;i<=size(v);i++){X=X*v[i];}
  matrix M=coef(ff,X);
  setring(@P);
  def RPE=imap(@R,EE);
  def RPN=imap(@R,NN);
  matrix Mp=imap(@R,M);
  poly g=Mp[2,1];
  if (size(Mp)!=2)
  {
    for(i=2;i<=size(Mp) div 2;i++)
    {
      g=gcd(g,Mp[2,i]);
    }
  }
  if (g==1){setring(RR); return(f);}
  else
  {
    def wg=factorize(g,2);
    if (wg[1][1]==1){setring(RR); return(f);}
    else
    {
      poly simp=1;
      int te;
      for(i=1;i<=size(wg[1]);i++)
      {
        te=inconsistent(RPE+wg[1][i],RPN);
        if(te)
        {
          simp=simp*(wg[1][i])^(wg[2][i]);
        }
      }
    }
    if (simp==1){setring(RR); return(f);}
    else
    {
      setring(RR);
      def simp0=imap(@P,simp);
      def f2=f/simp0;
      return(f2);
    }
  }
}

// pnormalf: reduces a polynomial f wrt a V(E)\V(N)
//           dividing by E and eliminating factors in N.
//           called in the ring @R,
//           operates in the ring @RP.
// input:
//         poly  f
//         ideal E  (depends only on the parameters)
//         ideal N  (depends only on the parameters)
//                  (E,N) represents V(E)\V(N)
//         optional:
// output: poly f2 reduced wrt to V(E)\V(N)
proc pnormalf(poly f, ideal E, ideal N)
"USAGE:   pnormalf(f,E,N);
          f: the polynomial to be reduced modulo V(E)\V(N)
          of a segment in the parameters.
          E: the null conditions ideal
          N: the non-null conditions
RETURN:   a reduced polynomial g of f, whose coefficients are reduced
          modulo E and having no factor in N.
NOTE: Should be called from ring Q[a][x].
          Ideals E and N must be given by polynomials
          in the parameters.
KEYWORDS: division, pdivi, reduce
EXAMPLE:  pnormalf; shows an example"
{
    def RR=basering;
    int te=0;
    if (defined(@P)){te=1;}
    else{setglobalrings();}
    setring(@RP);
    def fa=imap(RR,f);
    def Ea=imap(RR,E);
    def Na=imap(RR,N);
    option(redSB);
    Ea=std(Ea);
    def r=cld(reduce(fa,Ea));
    poly f1=r[1];
    f1=Ered(r[1],Ea,Na);
    setring(RR);
    def f2=imap(@RP,f1);
    if(te==0){kill @R; kill @RP; kill @P;}
    return(f2)
};
example
{ "EXAMPLE:"; echo = 2;
  ring R=(0,a,b,c),(x,y),dp;
  poly f=(b^2-1)*x^3*y+(c^2-1)*x*y^2+(c^2*b-b)*x+(a-bc)*y;
  ideal E=(c-1);
  ideal N=a-b;
  pnormalf(f,E,N);
}

// 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)
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));
}

// 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
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
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
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
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
// eliminates the polynomials whose lpp are divisible by some
// smaller one
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);
}

// redgbn: given a minimal basis (gb reducing) it
// reduces each polynomial wrt to V(E) \ V(N)
proc redgbn(ideal F, ideal E, ideal N)
{
  int te=0;
  if (defined(@P)==1){te=1;}
  ideal G=F;
  ideal H;
  int i;
  if (size(G)==0){return(ideal(0));}
  for (i=1;i<=size(G);i++)
  {
    H=delfromideal(G,i);
    G[i]=pnormalf(pdivi(G[i],H)[1],E,N);
    G[i]=primepartZ(G[i]);
  }
  if(te==1){setglobalrings();}
  return(G);
};

// eliminates repeated elements form an ideal
proc elimrepeated(ideal F)
{
  int i;
  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);
}

// equalideals
// input: 2 ideals F and G;
// output: 1 if they are identical (the same polynomials in the same order)
//         0 else
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
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);
}

//**************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
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);}
}

// postredgb: given a minimal basis (gb reducing) it
// reduces each polynomial wrt to the others
proc postredgb(ideal F)
{
  int te=0;
  if(defined(@P)==1){te=1;}
  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];
  }
  if(te==1){setglobalrings();}
  return(G);
}

//purpose ideal intersection called in @R and computed in @P
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
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));
}

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

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

// incquotient
// incremental quotient
// Input: ideal N: a Groebner basis of an ideal
//        poly f:
// Output: Na = N:<f>
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);
}

// eliminates the ith element from a list
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);
}

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);}
  }
}

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
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) and (i<=size(L)))
      {
        if (equalideals(L[i],M[i])==0){t=0;}
        i++;
      }
    }
    return(t);
  }
}

// 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
proc Prep(ideal N, ideal M)
{
  if (N[1]==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-representaion of V(N)\V(M) = V(ida)\V(idb)
// Assumed to work in the ring @P of the parameters
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));
}

// input: F a parametric ideal in Q[a][x]
// output: a rComprehensive Groebner System disjoint and reduced.
//      It uses Kapur-Sun-Wang algorithm, and with the options
//      can compute the homogenization before  (('can',0) or ( 'can',1))
//      and dehomogenize the result.
proc cgsdr(ideal F, list #)
"USAGE: cgsdr(F); To compute a disjoint, reduced CGS.
          cgsdr is the starting point of the fundamental routine grobcov.
          Inside grobcov it is used only with options 'can' set to 0,1 and
          not with options ('can',2).
          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:
            "can",0-1-2: The default value is "can",2. In this case no
                homogenization is done. With option ("can",0) the given
                basis is homogenized, and with option ("can",1) the
                whole given ideal is homogenized before computing the
                cgs and dehomogenized after.
                with option ("can",0) the homogenized basis is used
                with option ("can",1) the homogenized ideal is used
                with option ("can",2) the given basis is used
            "null",ideal E: The default is ('null',ideal(0)).
            "nonnull",ideal N: The default (nonnull,ideal(1)).
                When options 'null' and/or 'nonnull' are given, then
                the parameter space is restricted to V(E)\V(N).
            "comment",0-1: The default is ('comment',0). Setting ('comment',1)
                will provide information about the development of the
                computation.
            "out",0-1: 1 (default) the output segments are given as
                as difference of varieties.
                0: the output segments are given in P-representation
                and the segments grouped by lpp
                With options ("can",0) and ("can",1) the option ("out",1)
                is set to ("out,0) because it is not compatible.
          One can give none or whatever of these options.
          With the default options ("can",2,"out",1), only the
          Kapur-Sun-Wang algorithm is computed. This is very effectif
          but is only the starting point for the grobcov computation.
          When grobcov is computed, the call to cgsdr inside uses
          specific options that are more expensive ("can",0-1,"out",0).
RETURN:   Returns a list T describing a reduced and disjoint
          Comprehensive Groebner System (CGS),
          With option ("out",0)
           the segments are grouped by
           leading power products (lpp) of the reduced Groebner
           basis and given in P-representation.
           The returned list is of the form:
           (
             (lpp, (num,basis,segment),...,(num,basis,segment),lpp),
             ..,,
             (lpp, (num,basis,segment),...,(num,basis,segment),lpp)
           )
           The bases are the reduced Groebner bases (after normalization)
           for each point of the corresponding segment.

           The third element of each lpp segment is the lpp of the
           used ideal in the CGS as a string:
            with option ("can",0) the homogenized basis is used
            with option ("can",1) the homogenized ideal is used
            with option ("can",2) the given basis is used

          With option ("out",1) (default)
           only KSW is applied and segments are given as
           difference of varieties and are not grouped
           The returned list is of the form:
           (
             (E,N,B),..(E,N,B)
           )
           E is the null variety
           N is the nonnull variety
           segment = V(E)\V(N)
           B is the reduced Groebner basis

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"
{
  def RR=basering;
  setglobalrings();
  // INITIALIZING OPTIONS
  int i; int j;
  int can=2;
  int out=1;
  poly f;
  ideal B;
  def E=ideal(0);
  def N=ideal(1);
  int comment=0;
  int start=timer;
  list L=#;
  for(i=1;i<=size(L) div 2;i++)
  {
    if(L[2*i-1]=="null"){E=L[2*i];}
    else
    {
      if(L[2*i-1]=="nonnull"){N=L[2*i];}
      else
      {
        if(L[2*i-1]=="comment"){comment=L[2*i];}
        else
        {
          if(L[2*i-1]=="can"){can=L[2*i];}
          else
          {
            if(L[2*i-1]=="out"){out=L[2*i];}
          }
        }
      }
    }
  }
  //if(can==2){out=1;}
  B=F;
  if ((printlevel) and (comment==0)){comment=printlevel;}
  if((can<2) and (out>0)){"Option out,1 is not compatible with can,0,1"; out=0;}
  // DEFINING OPTIONS
  list LL;
  LL[1]="can";     LL[2]=can;
  LL[3]="comment"; LL[4]=comment;
  LL[5]="out";     LL[6]=out;
  LL[7]="null";    LL[8]=E;
  LL[9]="nonnull"; LL[10]=N;
  if(comment>=1)
  {
    "Begin cgsdr with options: "+string(LL);
  }
  int ish;
  for (i=1;i<=size(B);i++){ish=ishomog(B[i]); if(ish==0){break;};}
  if (ish)
  {
     if(comment>0){"The given system is homogneous";}
    can=0;
  }
  // ACTING DEPENDING ON OPTIONS
  if(can==2)
  {
    // WITHOUT HOMOHGENIZING
    if(comment>0){"Option of cgsdr: do not homogenize";}
    def GS=KSW(B,LL);
    setglobalrings();
  }
  else
  {
    if(can==1)
    {
      // COMPUTING THE HOMOGOENIZED IDEAL
      if(comment>0){"Homogenizing the whole ideal: option can=1"; }
      list RRL=ringlist(RR);
      RRL[3][1][1]="dp";
      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 B1=imap(RR,B);
      option(redSB);
      B1=std(B1);
      setring(RR);
      def B2=imap(RP,B1);
    }
    else
    { // (can=0)
      if(comment>0){"Homogenizing the basis: option can=0";}
      def B2=B;
    }
    // COMPUTING HOMOGENIZED CGS
    poly @t;
    ring H=0,@t,dp;
    def RH=RR+H;
    setring(RH);
    setglobalrings();
    def BH=imap(RR,B2);
    def LH=imap(RR,LL);
    for (i=1;i<=size(BH);i++)
    {
      BH[i]=homog(BH[i],@t);
    }
    if (comment>=1){"Homogenized system = "; BH;}
    def GSH=KSW(BH,LH);
    setglobalrings();
    // DEHOMOGENIZING THE RESULT
    if(out==0)
    {
      for (i=1;i<=size(GSH);i++)
      {
        GSH[i][1]=subst(GSH[i][1],@t,1);
        for(j=1;j<=size(GSH[i][2]);j++)
        {
          GSH[i][2][j][2]=subst(GSH[i][2][j][2],@t,1);
        }
      }
    }
    else
    {
      for (i=1;i<=size(GSH);i++)
      {
        GSH[i][3]=subst(GSH[i][3],@t,1);
        GSH[i][7]=subst(GSH[i][7],@t,1);
      }
    }
    setring(RR);
    def GS=imap(RH,GSH);
    setglobalrings();
    if(out==0)
    {
      for (i=1;i<=size(GS);i++)
      {
        GS[i][1]=postredgb(mingb(GS[i][1]));
        for(j=1;j<=size(GS[i][2]);j++)
        {
          GS[i][2][j][2]=postredgb(mingb(GS[i][2][j][2]));
        }
      }
    }
    else
    {
      for (i=1;i<=size(GS);i++)
      {
        if(GS[i][2]==1)
        {
          GS[i][3]=postredgb(mingb(GS[i][3]));
          GS[i][7]=postredgb(mingb(GS[i][7]));
        }
      }
    }
  }
  if(defined(@P)){kill @P; kill @R; kill @RP;}
  return(GS);
}
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);
}

// input:  internal routine called by cgsdr at the end to group the
//         lpp segments and improve the output
// output: grouped segments by lpp obtained in cgsdr
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])));
      }
    }
  }
  return(L);
}

// idcontains
// input: ideal p, ideal q
// output: 1 if p contains q,  0 otherwise
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) 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
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])
    {
      i++;
      if(i>size(Lp)){break;}
    }
    t=1;
    j=1;
    while ((t) and (j<=size(Lp)))
    {
      if (i==j){j++;}
      if ((j<=size(Lp)) and (memberpos(j,notsel)[1]==0))
      {

        if (idcontains(Lp[i],Lp[j]))
        {
          notsel[size(notsel)+1]=i;
          t=0;
        }
      }
      j++;
    }
    if (t){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)))
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
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) and (j<=size(C)))
      {
        if (equalideals(q,C[j][1]))
        {
          t=0;
          for (k=1;k<=size(C[j][2]);k++)
          {
            t1=1;
            //kill addq;
            //list addq;
            l=1;
            while((t1) and (l<=size(Q)))
            {
              if ((l!=i) and (memberpos(l,notQ)[1]==0))
              {
                if (idcontains(C[j][2][k],Q[l]))
                {
                  t1=0;
                }
              }
              l++;
            }
            if (t1)
            {
              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.
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)
        {
          for(k=1;k<=size(C[j][2]);k++)
          {
            if(idcontains(Q[i],C[j][2][k]))
            {
              te=0; break;
            }
          }
          if (te)
          {
            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){"addpartfine was ", used, " times used";}
  return(imap(@P,Q1));
}

// 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
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 #: optional
// output: the list
//   S=((lpp, generic basis, Prep, Crep),..,(lpp, generic basis, Prep, Crep))
//      where a Prep is ( (p1,(p11,..,p1k_1)),..,(pj,(pj1,..,p1k_j)) )
//            a Crep is ( ida, idb )
proc gcover(ideal F,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; ideal FF;
  list lpi;
  string lpph;
  list L=#;
  int canop=1;
  int extop=1;
  int repop=0;
  ideal E=ideal(0);;
  ideal N=ideal(1);;
  int comment;
  for(i=1;i<=size(L) div 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"){E=L[2*i];}
          else
          {
            if(L[2*i-1]=="nonnull"){N=L[2*i];}
            else
            {
              if (L[2*i-1]=="comment"){comment=L[2*i];}
            }
          }
        }
      }
    }
  }
  list GS; list GP;
  def RR=basering;
  GS=cgsdr(F,L); // "null",NW[1],"nonnull",NW[2],"cgs",CGS,"comment",comment);
  setglobalrings();
  int start=timer;
  GP=GS;
  ideal lppr;
  list LL;
  list S;
  poly sp;
  ideal BB;
  for (i=1;i<=size(GP);i++)
  {
    kill LL;
    list LL;
    lpp=GP[i][1];
    GPi2=GP[i][2];
    lpph=GP[i][3];
    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]);
      B=GPi2[LCU[k][1][1]][2]; // ATENTION last 1 has been changed to [2]
      LCU[k][1]=B;
    }
    //"Deciding if combine is needed";
    kill BB;
    ideal BB;
    tes=1; m=1;
    while((tes) and (m<=size(LCU[1][1])))
    {
      j=1;
      while((tes) and (j<=size(LCU)))
      {
        k=1;
        while((tes) and (k<=size(LCU)))
        {
          if(j!=k)
          {
            sp=pnormalf(pspol(LCU[j][1][m],LCU[k][1][m]),LCU[k][2],N);
            if(sp!=0){tes=0;}
          }
          k++;
        }        //setglobalrings();
        if(tes)
        {
          BB[m]=LCU[j][1][m];
        }
        j++;
      }
      if(tes==0){break;}
      m++;
    }    //"T_BB="; BB;
    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]=pnormalf(B[j],crep[1],N);
    }
    S[i]=list(lpp,B,prep,crep,lpph);
    if(comment>=1)
    {
      lpi[size(lpi)+1]=string("[",i,"]");
      lpi[size(lpi)+1]=S[i][1];
    }
  }
  if(comment>=1)
  {
    "Time in LCUnion + combine = ",timer-start;
    if(comment>=2){"lpp=",lpi};
  }
  if(defined(@P)==1){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",0-1,ext",0-1, "rep",0-1-2)
//            where
//            N is the null conditions ideal (if desired)
//            W is the ideal of non-null conditions (if desired)
//            The value of "can"is 1 by default and can be set to 0 if we do not
//            need to obtain the canonical GC, but only a GC.
//            The value of "ext" is 0 by default and so the generic representation
//             of the bases is given. It can be set to 1, and then the full
//             representation of the bases is given.
//            The value of "rep" is 0 by default, and then the segments
//            are given in canonical P-representation. It can be set to 1
//            and then they are given in canonical C-representation.
//            If it is set to 2, then both representations are given.
// 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 E: The default is ("null",ideal(0)).
            "nonnull",ideal N: The default ("nonnull",ideal(1)).
                When options "null" and/or "nonnull" are given, then
                the parameter space is restricted to V(E)\V(N).
            "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 shorter.
            "ext",0-1: The default is ("ext",0). With the default
                ("ext",0), only the generic representation is computed
                (single polynomials, but not specializing to non-zero at
                each point of the segment. With option ("ext",1) the
                full representation of the bases is computed (possible
                shaves) and sometimes a simpler result is obtained.
            "rep",0-1-2: The default is ("rep",0) and then the segments
                are given in canonical P-representation. Option ("rep",1)
                represents the segments in canonical C-representation,
                and option ("rep",2) gives both representations.
            "comment",0-3: The default is ("comment",0). Setting
                "comment" higher will provide information about the
                development of the computation.
          One can give none or whatever of these options.
RETURN:   The list
          (
           (lpp_1,basis_1,segment_1,lpph_1),
           ...
           (lpp_s,basis_s,segment_s,lpph_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.
          The lpph corresponds to the lpp of the homogenized ideal
          and is different for each segment. It is given as a string.

          Basis: to each element of lpp corresponds an I-regular function given
          in full representation (by option ("ext",1)) or in
          generic representation (default option ("ext",0)). The
          I-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, specializes to the corresponding
          element of the reduced Groebner basis. In the full representation
          at least one of the polynomials representing the I-regular
          function specializes to non-zero.

          With the default option ("rep",0) the representation of the
          segment is the P-representation.
          With option ("rep",1) the representation of the segment is
          the C-representation.
          With option ("rep",2) both representations of the segment are
          given.

          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.

          The C-representation of a segment is of the form
          (E,N) representing V(E)\V(N), and the ideals E and N are
          radical and N contains E.

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.
EXAMPLE:  grobcov; shows an example"
{
  list S; int i; int ish=1; list GBR; list BR; int j; int k;
  ideal idp; ideal idq; int s; ideal ext; list SS;
  ideal E; ideal N; int canop;  int extop; int repop;
  int comment=0; int m;
  def RR=basering;
  setglobalrings();
  list L0=#;
  int out=0;
  L0[size(L0)+1]="res"; L0[size(L0)+1]=ideal(1);
  // default options
  int start=timer;
  E=ideal(0);
  N=ideal(1);
  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=0; // extop=0 if only generic representation of the bases are to be computed
           // extop=1 if the full representation of the bases are to be computed
  for(i=1;i<=size(L0) div 2;i++)
  {
    if(L0[2*i-1]=="can"){canop=L0[2*i];}
    else
    {
      if(L0[2*i-1]=="ext"){extop=L0[2*i];}
      else
      {
        if(L0[2*i-1]=="rep"){repop=L0[2*i];}
        else
        {
          if(L0[2*i-1]=="null"){E=L0[2*i];}
          else
          {
            if(L0[2*i-1]=="nonnull"){N=L0[2*i];}
            else
            {
              if (L0[2*i-1]=="comment"){comment=L0[2*i];}
            }
          }
        }
      }
    }
  }
  if(not((canop==0) or (canop==1)))
  {
    "Option can = ",canop," is not supported. It is changed to can = 1";
    canop=1;
  }
  for(i=1;i<=size(L0) div 2;i++)
  {
    if(L0[2*i-1]=="can"){L0[2*i]=canop;}
  }
  if ((printlevel) and (comment==0)){comment=printlevel;}
  list LL;
  LL[1]="can";     LL[2]=canop;
  LL[3]="comment"; LL[4]=comment;
  LL[5]="out";     LL[6]=0;
  LL[7]="null";    LL[8]=E;
  LL[9]="nonnull"; LL[10]=N;
  LL[11]="ext";    LL[12]=extop;
  LL[13]="rep";    LL[14]=repop;
  if (comment>=1)
  {
    "Begin grobcov with options: ",string(LL);
  }
  kill S;
  def S=gcover(F,LL);
  // NOW extend
  if(extop)
  {
    S=extend(S,LL);
  }
  else
  {
    // NOW representation of the segments by option repop
    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],S[i][5]);
        nS[size(nS)+1]=Si;
      }
      kill S;
      def S=nS;
    }
    else
    {
      if(repop==1)
      {
        for(i=1;i<=size(S);i++)
        {
          Si=list(S[i][1],S[i][2],S[i][4],S[i][5]);
          nS[size(nS)+1]=Si;
        }
        kill S;
        def S=nS;
      }
      else
      {
        for(i=1;i<=size(S);i++)
        {
          Si=list(S[i][1],S[i][2],S[i][3],S[i][4],S[i][5]);
          nS[size(nS)+1]=Si;
        }
        kill S;
        def S=nS;
      }
    }
  }
  if (comment>=1)
  {
    "Time in grobcov = ", timer-start;
    "Number of segments of grobcov = ", size(S);
  }
  if(defined(@P)==1){kill @R; kill @P; kill @RP;}
  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. GC the grobcov of an ideal in generic representation of the
//        bases computed with option option ("rep",2).
// output The grobcov in full representation.
// option ("comment",1) shows the time.
proc extend(list GC, list #);
"USAGE:   extend(GC); When the grobcov of an ideal has been computed
          with the default option ("ext",0) and the explicit option
          ("rep",2) (which is not the default), then one can call
          extend (GC) (and options) to obtain the full representation
          of the bases. With the default option ("ext",0) only the
          generic representation of the bases are computed, and one can
          obtain the full representation using extend.
            "rep",0-1-2: The default is ("rep",0) and then the segments
                are given in canonical P-representation. Option ("rep",1)
                represents the segments in canonical C-representation,
                and option ("rep",2) gives both representations.
            "comment",0-1: The default is ("comment",0). Setting
                "comment" higher will provide information about the
                time used in the computation.
          One can give none or whatever of these options.
RETURN:   The list
          (
           (lpp_1,basis_1,segment_1,lpph_1),
           ...
           (lpp_s,basis_s,segment_s,lpph_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.
          The lpph corresponds to the lpp of the homogenized ideal
          and is different for each segment. It is given as a string.

          Basis: to each element of lpp corresponds an I-regular function given
          in full representation. The
          I-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, specializes to the corresponding
          element of the reduced Groebner basis. In the full representation
          at least one of the polynomials representing the I-regular
          function specializes to non-zero.

          With the default option ("rep",0) the segments are given
          in P-representation.
          With option ("rep",1) the segments are given
          in C-representation.
          With option ("rep",2) both representations of the segments are
          given.

          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.

          The C-representation of a segment is of the form
          (E,N) representing V(E)\V(N), and the ideals E and N are
          radical and N contains E.

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, full representation.
EXAMPLE:  extend; shows an example"
{
  list L=#;
  list S=GC;
  ideal idp;
  ideal idq;
  int i; int j; int m; int s;
  m=0; i=1;
  while((i<=size(S)) and (m==0))
  {
    if(typeof(S[i][2])=="list"){m=1;}
    i++;
  }
  if(m==1){"Warning! grobcov has already extended bases"; return(S);}
  if(size(GC[1])!=5){"Warning! extend make sense only when grobcov has been called with options 'rep',2,'ext',0"; " "; return();}
  int repop=0;
  int start3=timer;
  int comment;
  for(i=1;i<=size(L) div 2;i++)
  {
    if(L[2*i-1]=="comment"){comment=L[2*i];}
    else
    {
      if(L[2*i-1]=="rep"){repop=L[2*i];}
    }
  }
  poly leadc;
  poly ext;
  int te=0;
  list SS;
  def R=basering;
  if (defined(@R)){te=1;}
  else{setglobalrings();}
  // Now extend
  for (i=1;i<=size(S);i++)
  {
    m=size(S[i][2]);
     for (j=1;j<=m;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=extend0(S[i][2][j],idp,idq);
        if (typeof(ext)=="poly")
        {
          S[i][2][j]=pnormalf(ext,idp,idq);
        }
        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]=pnormalf(ext[s],idp,idq);
            }
            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;
          }
        }
      }
    }
  }
  // NOW representation of the segments by option repop
  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],S[i][5]);
      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],S[i][5]);
        nS[size(nS)+1]=Si;
      }
      S=nS;
    }
    else
    {
      for(i=1;i<=size(S);i++)
      {
        Si=list(S[i][1],S[i][2],S[i][3],S[i][4],S[i][5]);
        nS[size(nS)+1]=Si;
      }

    }
  }
  if(comment>=1){"Time in extend = ",timer-start3;}
  if(te==0){kill @R; kill @RP; kill @P;}
  return(S);
}
example
{
  ring R=(0,a0,b0,c0,a1,b1,c1,a2,b2,c2),(x), dp;
  short=0;
  ideal S=a0*x^2+b0*x+c0,
          a1*x^2+b1*x+c1,
          a2*x^2+b2*x+c2;
  "System S="; S;

  def GCS=grobcov(S,"rep",2,"comment",1);
  "grobcov(S,'rep',2,'comment',1)="; GCS;
  def FGC=extend(GCS,"rep",0,"comment",1);
  "Full representation="; FGC;
}


// nonzerodivisor
// input:
//    poly g in K[a],
//    list P=(p_1,..p_r) representing a minimal prime decomposition
// output
//    poly f such that f notin p_i for all i and
//           g-f in p_i for all i such that g notin p_i
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);
}

// deltai
// 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.
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));
}

// combine
// input: a list of pairs ((p1,P1),..,(pr,Pr)) where
//    ideal pi is a prime component
//    poly Pi is the polynomial in Q[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)
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);
  f=primepartZ(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
proc elimconstfac(poly f)
{
  int cond; int i; int j; int t;
  if (f==0){return(f);}
  def RR=basering;
  setring(@R);
  def ff=imap(RR,f);
  def 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);
};

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

// monoms
proc monoms(poly f)
{
  list L;
  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[size(L)+1]=list(lc,lp);
    i++;
  }
  return(L);
}

// extend0
// 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:
proc extend0(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=pnormalf(fout,idp,W);
      if(ncols(NewMonoms[1][1])>1){idout[j]=fout;}
    }
    if(ncols(NewMonoms[1][1])==1){return(fout);} else{return(idout);}
  }
  else
  {
    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);}
  }
}

// findindexpolys
// 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.
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,
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);
}

// selectregularfun
// 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
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))){te=1; break;}
    }
    if(te){break;}
  }
  ci=T[1];
  def Cs=submat(C1,1..2,ci);
  setring(RR);
  return(imap(@P,Cs));
}

// searchinlist
// input:
//   intvec c:
//   list L=( (c1,T1),..(ck,Tk) )
//      where the c's are assumed to be intvects
// output:
//   object T with index c
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);
}

// comb: the list of combinations of elements (1,..n) of order p
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.
proc selectminsheaves(list L)
{
  list C=allsheaves(L);
  return(smsheaves(C[1],C[2]));
}

// smsheaves
// 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).
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
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.
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);
}

// pos
// Input:  intvec p of zeros and ones
// Output: intvec W of the positions where p has ones.
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.
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
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);
}

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

// intersp: computes the intersection of the ideals in S in @P
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));
}

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

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

// selectextendcoef
// 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);
proc selectextendcoef(matrix CC, ideal ida, ideal idb)
{
  def RR=basering;
  setring(@P);
  def ca=imap(RR,CC);
  def E0=imap(RR,ida);
  ideal E;
  def N=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))
  {
    com=comb(r,i);
    j=1;
    while((j<=size(com)) and (te))
    {
      E=E0;
      c=com[j];
      for (k=1;k<=i;k++)
      {
        E=E+ca[2,c[k]];
      }
      prep=Prep(E,N);
      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){"error: extendcoef does not extend to the whole S";}
  setring(RR);
  return(imap(@P,caout));
}

// input:
//   ideal E1: in some basering (depends only on the parameters)
//   ideal E2: in some basering (depends only on the parameters)
// output:
//   ideal Ep=E1+E2; computed in P
proc plusP(ideal E1,ideal E2)
{
  def RR=basering;
  setring(@P);
  def E1p=imap(RR,E1);
  def E2p=imap(RR,E2);
  def Ep=E1p+E2p;
  setring(RR);
  return(imap(@P,Ep));
}

// reform
// 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
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])
    {
      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);
}

// nonnullCrep
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);
}

// precombint
// 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
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) 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){"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).
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]))
    {
      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){selind=i1; break;}
  }
  vv=tes,selind;
  return(vv);
}

// minAssGTZ eliminating denominators
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);
}

//********************* Begin KapurSunWang *************************

// inconsistent
// Input:
//   ideal E: of null conditions
//   ideal N: of non-null conditions representing V(E)\V(N)
// Output:
//   1 if V(E) \V(N) = empty
//   0 if not
proc inconsistent(ideal E, ideal N)
{
  int j;
  int te=1;
  def R=basering;
  setring(@P);
  def EP=imap(R,E);
  def NP=imap(R,N);
  poly @t;
  ring H=0,@t,dp;
  def RH=@P+H;
  setring(RH);
  def EH=imap(@P,EP);
  def NH=imap(@P,NP);
  ideal G;
  j=1;
  while((te==1) and j<=size(NH))
  {
    G=EH+(1-@t*NH[j]);
    option(redSB);
    G=std(G);
    if (G[1]!=1){te=0;}
    j++;
  }
  setring(R);
  return(te);
}

// MDBasis: Minimal Dickson Basis
proc MDBasis(ideal G)
{
  int i; int j; int te=1;
  G=sortideal(G);
  ideal MD=G[1];
  poly lm;
  for (i=2;i<=size(G);i++)
  {
    te=1;
    lm=leadmonom(G[i]);
    j=1;
    while ((te==1) and (j<=size(MD)))
    {
      if (lm/leadmonom(MD[j])!=0){te=0;}
      j++;
    }
    if (te==1)
    {
      MD[size(MD)+1]=(G[i]);
    }
  }
  return(MD);
}

// primepartZ
proc primepartZ(poly f);
{
  def R=basering;
  def cp=content(f);
  def fp=f/cp;
  return(fp);
}

// LCMLC
proc LCMLC(ideal H)
{
  int i;
  def R=basering;
  setring(@RP);
  def HH=imap(R,H);
  poly h=1;
  for (i=1;i<=size(HH);i++)
  {
    h=lcm(h,HH[i]);
  }
  setring(R);
  def hh=imap(@RP,h);
  return(hh);
}

// KSW: Kapur-Sun-Wang algorithm for computing a CGS
// Input:
//   F:   parametric ideal to be discussed
//   Options:
//     "out",0 Transforms the description of the segments into
//     canonical P-representation form.
//     "out",1 Original KSW routine describing the segments as
//     difference of varieties
//   The ideal must be defined on C[parameters][variables]
// Output:
//   With option "out",0 :
//     ((lpp,
//       (1,B,((p_1,(p_11,..,p_1k_1)),..,(p_s,(p_s1,..,p_sk_s)))),
//       string(lpp)
//      )
//      ,..,
//      (lpp,
//       (k,B,((p_1,(p_11,..,p_1k_1)),..,(p_s,(p_s1,..,p_sk_s)))),
//       string(lpp))
//      )
//     )
//   With option "out",1 ((default, original KSW) (shorter to be computed,
//                    but without canonical description of the segments.
//     ((B,E,N),..,(B,E,N))
proc KSW(ideal F, list #)
{
  setglobalrings();
  int start=timer;
  ideal E=ideal(0);
  ideal N=ideal(1);
  int comment=0;
  int out=1;
  int i;
  def L=#;
  if (size(L)>0)
  {
    for (i=1;i<=size(L)/2;i++)
    {
      if (L[2*i-1]=="null"){E=L[2*i];}
      else
      {
        if (L[2*i-1]=="nonnull"){N=L[2*i];}
        else
        {
          if (L[2*i-1]=="comment"){comment=L[2*i];}
          else
          {
            if (L[2*i-1]=="out"){out=L[2*i];}
          }
        }
      }
    }
  }
  if (comment>0){"Begin KSW with null = ",string(E)," nonnull = ",string(N);}
  def CG=KSW0(F,E,N,comment);
  if (comment>0)
  {
    "Number of segments in KSW (total) = ",size(CG);
    "Time in KSW = ",timer-start;
  }
  if(out==0)
  {
    CG=KSWtocgsdr(CG);
    CG=groupKSWsegments(CG);
    if (comment>0)
    {
      "Number of lpp segments = ",size(CG);
      "Time in KSW + group + Prep = ",timer-start;
    }
  }
  if(defined(@P)){kill @P; kill @R; kill @RP;}
  return(CG);
}

// sqf
// This is for releases of Singular before 3-5-1
// proc sqf(poly f)
// {
//  def RR=basering;
//  setring(@P);
//  def ff=imap(RR,f);
//  def G=sqrfree(ff);
//  poly fff=1;
//  int i;
//  for (i=1;i<=size(G);i++)
//  {
//    fff=fff*G[i];
//  }
//  setring(RR);
//   def ffff=imap(@P,fff);
//   return(ffff);
// }

// sqf
proc sqf(poly f)
{
  def RR=basering;
  setring(@P);
  def ff=imap(RR,f);
  poly fff=sqrfree(ff,3);
  setring(RR);
  def ffff=imap(@P,fff);
  return(ffff);
}



// KSW0: Kapur-Sun-Wang algorithm for computing a CGS, called by KSW
// Input:
//   F:   parametric ideal to be discussed
//   Options:
//   The ideal must be defined on C[parameters][variables]
// Output:
proc KSW0(ideal F, ideal E, ideal N, int comment)
{
  def R=basering;
  int i; int j; list emp;
  list CGS;
  ideal N0;
  for (i=1;i<=size(N);i++)
  {
    N0[i]=sqf(N[i]);
  }
  ideal E0;
  for (i=1;i<=size(E);i++)
  {
    E0[i]=sqf(leadcoef(E[i]));
  }
  setring(@P);
  ideal E1=imap(R,E0);
  E1=std(E1);
  ideal N1=imap(R,N0);
  N1=std(N1);
  setring(R);
  E0=imap(@P,E1);
  N0=imap(@P,N1);
//   E0=elimrepeated(E0);
//   N0=elimrepeated(N0);
  if (inconsistent(E0,N0)==1)
  {
    return(emp);
  }
  setring(@RP);
  def FRP=imap(R,F);
  def ERP=imap(R,E);
  FRP=FRP+ERP;
  option(redSB);
  def GRP=std(FRP);
  setring(R);
  def G=imap(@RP,GRP);
  if (memberpos(1,G)[1]==1)
  {
    if(comment>1){"Basis 1 is found"; E; N;}
    return(E0,N0,ideal(1));
  }
  ideal Gr; ideal Gm; ideal GM;
  for (i=1;i<=size(G);i++)
  {
    if (variables(G[i])[1]==0){Gr[size(Gr)+1]=G[i];}
    else{Gm[size(Gm)+1]=G[i];}
  }
  ideal Gr0;
  for (i=1;i<=size(Gr);i++)
  {
    Gr0[i]=sqf(Gr[i]);
  }


  Gr=elimrepeated(Gr0);
  ideal GrN;
  for (i=1;i<=size(Gr);i++)
   {
    for (j=1;j<=size(N0);j++)
    {
      GrN[size(GrN)+1]=sqf(Gr[i]*N0[j]);
    }
  }
  if (inconsistent(E,GrN)){;}
  else
  {
    if(comment>1){"Basis 1 is found in a branch with arguments"; E; GrN;}
    CGS[size(CGS)+1]=list(E,GrN,ideal(1));
  }
  if (inconsistent(Gr,N0)){return(CGS);}
  GM=Gm;
  Gm=MDBasis(Gm);
  ideal H;
  for (i=1;i<=size(Gm);i++)
  {
    H[i]=sqf(leadcoef(Gm[i]));
  }
  H=facvar(H);
  poly h=sqf(LCMLC(H));
  if(comment>1){"H = "; H; "h = "; h;}
  ideal Nh=N0;
  if(size(N0)==0){Nh=h;}
  else
  {
    for (i=1;i<=size(N0);i++)
    {
      Nh[i]=sqf(N0[i]*h);
    }
  }
  if (inconsistent(Gr,Nh)){;}
  else
  {
    CGS[size(CGS)+1]=list(Gr,Nh,Gm);
  }
  poly hc=1;
  list KS;
  ideal GrHi;
  for (i=1;i<=size(H);i++)
  {
    kill GrHi;
    ideal GrHi;
    Nh=N0;
    if (i>1){hc=sqf(hc*H[i-1]);}
    for (j=1;j<=size(N0);j++){Nh[j]=sqf(N0[j]*hc);}
    if (equalideals(Gr,ideal(0))==1){GrHi=H[i];}
    else {GrHi=Gr,H[i];}
//     else {for (j=1;j<=size(Gr);j++){GrHi[size(GrHi)+1]=Gr[j]*H[i];}}
    if(comment>1){"Call to KSW with arguments "; GM; GrHi;  Nh;}
    KS=KSW0(GM,GrHi,Nh,comment);
    for (j=1;j<=size(KS);j++)
    {
      CGS[size(CGS)+1]=KS[j];
    }
    if(comment>1){"CGS after KSW = "; CGS;}
  }
  return(CGS);
}

// KSWtocgsdr
proc KSWtocgsdr(list L)
{
  int i; list CG; ideal B; ideal lpp; int j; list NKrep;
  for(i=1;i<=size(L);i++)
  {
    B=redgbn(L[i][3],L[i][1],L[i][2]);
    lpp=ideal(0);
    for(j=1;j<=size(B);j++)
    {
      lpp[j]=leadmonom(B[j]);
    }
    NKrep=KtoPrep(L[i][1],L[i][2]);
    CG[i]=list(lpp,B,NKrep);
  }
  return(CG);
}

// KtoPrep
// Computes the P-representaion of a K-representation (N,W) of a set
// input:
//    ideal E (null conditions)
//    ideal N (non-null conditions ideal)
// output:
//    the ((p_1,(p_11,..,p_1k_1)),..,(p_r,(p_r1,..,p_rk_r)));
//    the Prep of V(N) \ V(W)
proc KtoPrep(ideal N, ideal W)
{
  int i; int j;
  if (N[1]==1)
  {
    L0[1]=list(ideal(1),list(ideal(1)));
    return(L0);
  }
  def RR=basering;
  setring(@P);
  ideal B; int te; poly f;
  ideal Np=imap(RR,N);
  ideal Wp=imap(RR,W);
  list L;
  list L0; list T0;
  L0=minGTZ(Np);
  for(j=1;j<=size(L0);j++)
  {
    option(redSB);
    L0[j]=std(L0[j]);
  }
  for(i=1;i<=size(L0);i++)
  {
    if(inconsistent(L0[i],Wp)==0)
    {
      B=L0[i]+Wp;
      T0=minGTZ(B);
      option(redSB);
      for(j=1;j<=size(T0);j++)
      {
        T0[j]=std(T0[j]);
      }
      L[size(L)+1]=list(L0[i],T0);
    }
  }
  setring(RR);
  def LL=imap(@P,L);
  return(LL);
}

// groupKSWsegments
// input:  the list of vertices of KSW
// output: the same terminal vertices grouped by lpp
proc groupKSWsegments(list T)
{
  int i; int j;
  list L;
  list lpp; list lppor;
  list kk;
  lpp[1]=T[1][1]; j=1;
  lppor[1]=intvec(1);
  for(i=2;i<=size(T);i++)
  {
    kk=memberpos(T[i][1],lpp);
    if(kk[1]==0){j++; lpp[j]=T[i][1]; lppor[j]=intvec(i);}
    else{lppor[kk[2]][size(lppor[kk[2]])+1]=i;}
  }
  list ll;
  for (j=1;j<=size(lpp);j++)
  {
    kill ll; list ll;
    for(i=1;i<=size(lppor[j]);i++)
    {
      ll[size(ll)+1]=list(i,T[lppor[j][i]][2],T[lppor[j][i]][3]);
    }
    L[j]=list(lpp[j],ll,string(lpp[j]));
  }
  return(L);
}

//********************* End KapurSunWang *************************
;
//********************* Begin locus2d ****************************

// selfindimsols
// auxilliary routine called by locus2d
// input:  L the list of the Grobner Cover
// output: S the list of the union of segments where only a finite number
//         of solutions exists.
//         Supposed to be the set of points of the parameter space with
//         non degenerate solutions, for example in
//         automatic discovering of geometric theorems
proc selfindimsols(list L)
{
  int te=0;
  if (defined(@R)){te=1;}
  if(te==0){setglobalrings();}
  int i; int j;
  ideal v=variables(L[1][2]);
  ideal vv;
  for(i=2;i<=size(L);i++)
  {
    vv=variables(L[i][2]);
    for(j=1;j<=size(vv);j++)
    {
      if(memberpos(vv[j],v)[1]==0)
      {
        v[size(v)+1]=vv[j];
      }
    }
  }
  v=elimintfromideal(v);
  int nvartot=size(v);
  ideal lpp;
  int isovarlpp;
  ideal empty;
  list LL;
  ideal B;
  list SL;
  for (i=1;i<=size(L);i++)
  {
    lpp=L[i][1];
    isovarlpp=0;
    for (j=1;j<=size(lpp);j++)
    {
      if (size(variables(lpp[j]))==1)
      {
        isovarlpp=isovarlpp+1;
      }
    }
    if (isovarlpp==nvartot)
    {
      for(j=1;j<=size(L[i][3]);j++)
      {
        B=L[i][2],L[i][3][j][1];
        if(size(L[i][3][j][1])==1)
        {
          if(indepparameters(B))
          {
            SL=L[i][3][j];
            SL[3]="Special";
            LL[size(LL)+1]=SL;
          }
          else
          {
            LL[size(LL)+1]=L[i][3][j];
          }
        }
        else
        {
          LL[size(LL)+1]=L[i][3][j];
        }
      }
    }
  }
  if(te==0){kill @R; kill @P; kill @RP};
  return(LL);
}

// locus2d: Special routine for determining the locus of points
// of a two dimensional object. Given an ideal J with two
// parameters (a,b) and so many variables as needed, representing
// the system determining the locus of points (a,b) who verify
// certain geometrical properties, computing the grobcov of
// J and applying to it locus2d, determines the locus.
// input:
//    list GC, the output of grobcov
// output:
//    list, the locus of points of the parameter-space
//          for which the number of solutions in the variables
//          is finite.
//          If some component corresponds to a fixed single
//          solution in the variables but to a curve of the
//          parameter-sapace, then "Special" stands as
//          the third element of the component
//    ((p1,(p11,..p1s_1)),..,(pk,(pk1,..pks_k))
//    Possibly some component can be  (p1,(p11,..p1s_1),"Special")
//    These components of the locus correspond to locus curves
//    determined by a single or a finite number of points of
//    the geometrical construction.
proc locus2d(list GC)
"USAGE:   locus2d(G);
          The argument must be the grobcov of a two dimensional
          locus parametrical system with two parameters (a,b)
          and so many variables as needed, representing the locus
          points (a,b) who verify certain geometrical properties.
          Possibly some component can be  (p1,(p11,..p1s_1),'Special')
          These components of the locus correspond to locus curves
          determined by a single or a finite number of points of
          the geometrical construction.
RETURN:   The two dimensional locus.
NOTE:     It can only be called after computing the grobcov of the
          parametrical ideal in generic representation ('ext',0),
          which is the default.
          The basering R, must be of the form Q[a,b][x,y,..].
KEYWORDS: geometrical locus, locus, loci.
EXAMPLE:  locus2d; shows an example"
{
  def R=basering;
  setglobalrings();
  def LL=selfindimsols(GC);
  setring(@P);
  def L=imap(R,LL);
  int i; int j; int k; int n;
  list LL;
  intvec Lprals;
  intvec Ldep;
  list empty;
  poly f;
  list Ladd;
  intvec Lp;
  ideal N;
  intvec si;
  intvec sj;
  intvec elimin;
  for(i=1;i<=size(L);i++)
  {
    if(size(L[i][1])==1)
    {
      if(Lprals==intvec(0)){Lprals=i;}
      else{Lprals=Lprals,i;}
    }
    else
    {
      if(Ldep==intvec(0)){Ldep=i;}
      else{Ldep=Ldep,i;}
    }
  }
  for(i=1;i<=size(Lprals);i++)
  {
    Lp=Lprals[i];
    if(Ldep!=0)
    {
      for(j=1;j<=size(Ldep);j++)
      {
        N=L[Ldep[j]][1];
        attrib(N,"isSB",1);
        f=reduce(L[Lprals[i]][1][1],N);
        if(f==0)
        {
          Lp=Lp,Ldep[j];
        }
      }
    }
    Ladd[size(Ladd)+1]=Lp;
  }
  list Lfi;
  list La;
  list Lb;
  for (i=1;i<=size(Ladd);i++)
  {
    si=Ladd[i][1];
    n=size(L[si[1]][2]);
    kill elimin;
    intvec elimin;
    for (j=2;j<=size(Ladd[i]);j++)
    {
      sj=Ladd[i][j];
      for(k=1;k<=n;k++)
      {
        if (equalideals(L[sj][1],L[si[1]][2][k])==1)
        {
          if(elimin==intvec(0)){elimin=k;}
          else{elimin=elimin,k;}
        }
      }
    }
    kill Lb; list Lb;
    for (k=1;k<=n;k++)
    {
      if (not(memberpos(k,elimin)[1]))
      {
        Lb[size(Lb)+1]=L[si[1]][2][k];
      }
    }
    if (size(Lb)==0){Lb=ideal(1);}
    La=list(L[si[1]][1],Lb);
    if(size(L[si[1]])==3){La[3]=L[si[1]][3];}
    Lfi[size(Lfi)+1]=La;
  }
  setring(R);
  list Lout=imap(@P,Lfi);
  kill @R; kill @RP; kill @P;
  return(Lout);
}
example
{"EXAMPLE:"; echo = 2;
  ring R=(0,a,b),(x,y),dp;
  short=0;
  ideal H=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
  def G=grobcov(H);
  "grobcov(H)="; G; " ";
  def Gp=locus2d(G);
  "locus2d(G)="; Gp;
}

// locus2dto: Transforms the output of locus2d to a string that
//    can be reed from different computational systems.
// input:
//   list L: The output of locus2d
// output:
//   string s: The output of locus2d converted to a string readable
//             by other programs
proc locus2dto(list L)
"USAGE:   locus2dto(G);
          The argument must be the output of locus2d  of a two dimensional
          locus parametrical system with two parameters (a,b)
          and so many variables as needed, representing the locus
          points (a,b) who verify certain geometrical properties.
          It transforms the output to a string in standard form
          readable in many languages (Geogebra).

RETURN: The two dimensional locus in string standard form
NOTE:     It can only be called after computing the locus2d(grobcov(F)) of the
          parametrical ideal.
          The basering R, must be of the form Q[a,b][x,y,..].
KEYWORDS: geometrical locus, locus, loci.
EXAMPLE:  locus2dto; shows an example"
{
  int i; int j; int k;
  string s;
  s="[";
  ideal p;
  ideal q;
  for(i=1;i<=size(L);i++)
  {
    s=string(s,"[[");
    for (j=1;j<=size(L[i][1]);j++)
    {
      s=string(s,L[i][1][j],",");
    }
    s[size(s)]="]";
    s=string(s,",[");
    for(j=1;j<=size(L[i][2]);j++)
    {
      s=string(s,"[");
      for(k=1;k<=size(L[i][2][j]);k++)
      {
        s=string(s,L[i][2][j][k],",");
      }
      s[size(s)]="]";
      s=string(s,",");
    }
    s[size(s)]="]";
    s=string(s,"],");
    if(size(L[i])==3)
    {
      s[size(s)]=",";
      s=string(s,"[",L[i][3],"]],");
    }
  }
  s[size(s)]="]";
  return(s);
}
example
{"EXAMPLE:"; echo = 2;
  ring R=(0,a,b),(x,y),dp;
  short=0;
  ideal H=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
  def G=grobcov(H);
  "grobcov(H)="; G; " ";
  def Gp=locus2d(G);
  "locus2d(G)="; Gp;
  def L=locus2dto(Gp); " ";
  "locus2dto(Gp)="; L;
}

// indepparameters
// Auxiliary routine to detect "Special" components of the locus2d
// Input: ideal B
// Output:
//   1 if the solutions of the ideal do not depend on the parameters
//   0 if they depend
proc indepparameters(ideal B)
{
  def R=basering;
  ideal B0; ideal B00;
  int te;
  int i; int j;
  list s;
  poly t;
  ideal v=variables(B); // all the variables on B but not the parameters
  setring(@RP);
  ideal vv=imap(R,v);
  def BP=imap(R,B);
  option(redSB);
  BP=std(BP);
  setring(R);
  B0=imap(@RP,BP);
  for(i=1;i<=size(B0);i++)
  {
    if (equalideals(variables(B0[i]),ideal(0))){;}
    else {B00[size(B00)+1]=B0[i];}
  }
  for(i=1;i<=size(B00);i++)
  {
    s=factorize(B00[i]);
    for(j=1;j<=size(s[1]);j++)
    {
      if (equalideals(variables(s[1][j]),ideal(0))){;}
      else{B00[i]=s[1][j];}
    }
  }
  setring(@RP);
  BP=imap(R,B00);
  ideal vp=variables(BP);
  if(equalideals(vv,vp)){te=1;} else{te=0;}
  setring(R);
  return(te);
}

// lsolve
proc lsolve(ideal B)
{
  int i;
  list L;
  matrix c;
  def v=variables(B);
  ideal vi;
  poly v0;
  int te=1;
  i=1;
  while ((i<=size(B)) and te==1)
  {
    vi=variables(B[i]);
    if (size(vi)==1)
    {
      v0=vi[1];
      //"B[i]="; B[i];
      c=coeffs(B[i],v0);
      if (size(c)==2)
      {
        L[size(L)+1]=list(v0,-c[1,1]/c[2,1]);
      }
      else{te=0;}
    }
    else{te=0;}
    i++;
  }
  if(te==1){return(L);}
}