singular-data 1:4.1.0-p3+ds-2build1 / usr / share / singular / LIB / grobcov.lib

version="version grobcov.lib 4.0.3.4 Oct_2016 "; // $Id: 07e6e1e61cfc6becad39801eba0c7c567a96bf2e $
           // version M;  October_2016;
category="General purpose";
info="
LIBRARY:  grobcov.lib  Oktober 2016 Groebner Cover for parametric ideals.

          Groebner Cover for parametric ideals.
          Comprehensive Groebner Systems, Groebner Cover, Canonical Forms, Parametric Polynomial Systems,
          Dynamic Geometry, Loci, Envelop, Constructible sets. See:
          A. Montes A, M. Wibmer,
          \"Groebner Bases for Polynomial Systems with parameters\",
          Journal of Symbolic Computation 45 (2010) 1391-1425.
          (https://www.mat.upc.edu//en/people/antonio.montes/).

IMPORTANT: The book,  not yet published:
           A. Montes. \"Discussing Parametric Polynomial Systems: The Groebner Cover\"
          can be used as a user manual of all the routines included in this library.
          It defines and proves all the theoretic results used here, and shows examples of all the routines.
          It will be published soon.

AUTHORS:  Antonio Montes (Universitat Politecnica de Catalunya),
          Hans Schoenemann (Techische Universitaet Kaiserslautern).

OVERVIEW: In 2010, the library was designed to contain Montes-Wibmer's algorithms for computing the
          Canonical Groebner Cover of a parametric ideal.
          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
          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 actually 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.
          The library has evolved to include new applications of the Groebner Cover, and new theoretical
          developments have been done.

          The actual version also includes a routine (ConsLevels) for computing the canonical form of a
          constructible set, given as a union of locally closed sets. It is used in the new version for the
          computation of loci and envelops. It determines the canonical locally closed level sets of a
          constructible set. It is described in:
           J.M. Brunat, A. Montes, \"Computing the canonical representation of constructible sets\".
          Math.  Comput. Sci. (2016) 19: 165-178.

          A new set of routines (locus, locusdg, locusto) has been included to compute loci of points.
          The routines are used in the Dynamic Geometry software Geogebra. They are described in:
          M.A. Abanades, F. Botana, A. Montes, T. Recio:
          \''An Algebraic Taxonomy for Locus Computation in Dynamic Geometry\''.
          Computer-Aided Design 56 (2014) 22-33.

          Recently also routines for computing the generalized envelop of a family of hyper-surfaces (envelop),
          to be used in Dynamic Geometry, has been included and is described in the book (not yet published)
          A. Montes. \"Discussing Parametric Polynomial Systems: The Groebner Cover\"

          This version was finished on 31/10/2016

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.

PROCEDURES:
setglobalrings();  Generates the global rings @R, @P and @RP 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, locus, locusdg and
          killed before the output.
          In the actual version, after the call of setglobalrings on Q[a][x], the public names of the defined ideals
          generated by setglobalrings are Grobcov::@R, Grobcov::@P,  Grobcov::@RP.

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 can also be used independently if only a CGS is required. It is a more efficient routine
          than buildtree (the own routine of 2010 that is no more used).
          Now, Kapur-Sun-Wang (KSW) algorithm is used.

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.

Crep(N,M); Computes the canonical C-representation of V(N) - V(M). It can be called in Q[a] or in Q[a][x],
          but the ideals N,M can only contain parameters of Q[a].

Prep(N,M); Computes the canonical P-representation of V(N) - V(M). It can be called in Q[a] or in Q[a][x],
          but the ideals N,M can only contain parameters of Q[a].

PtoCrep(L)  Starting from the canonical Prep of a locally closed set computes its Crep.

extendpoly(f,p,q); Given the generic representation f of an I-regular function F defined by poly f on V(p) - V(q)
          it returns its full representation.

extendGC(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 extendGC(GC) (and options) to obtain the
          full representation of the bases. With the default option (\"ext\",0) only the generic representation of
          the bases is computed, and one can obtain the full representation using extendGC.

ConsLevels(L); Given a list L of locally closed sets, it returns the closures of the canonical levels of the constructible
          set and its complements of the union of them. It is described in
          J.M. Brunat, A. Montes, \"Computing the canonical representation of constructible sets\".
          Math.  Comput. Sci. (2016) 19: 165-178.

ConsLevelsToLevels(L);Transforms the output of ConsLevels into the proper Levels of the constructible set.

locus(G); Special routine for determining the geometrical locus of points verifying given conditions. Given a
          parametric ideal J with parameters (x,y) and variables (x_1,..,xn), representing the system determining
          the locus of points (x,y) who verify certain properties, one can apply locus to the output of  grobcov(J).
          locus determines the different classes of locus components, following the taxonomy described in
          M. Abanades, F. Botana, A. Montes, T. Recio,
          \"An Algebraic Taxonomy for Locus Computation in Dynamic Geometry\",
          Computer-Aided Design 56 (2014) 22-33.

          The components can be \"Normal\", \"Special\", \"Accumulation\", \"Degenerate\".
          The output are the components given in P-canonical form. It also detects automatically a possible point
          that is to be avoided by the mover, whose coordinates must be the last coordinates in the definition of
          the ring. If such a point is detected, then it eliminates the segments of the grobcov depending on the
          point that is to be avoided.

locusdg(G); Is a special routine that determines the  \"Relevant\" components of the locus in dynamic geometry.
          It is to be called to the output of locus and selects from it the useful components.

envelop(F,C); Special routine for determining the envelop of a family of hyper-surfaces F in Q[x1,..,xn][t1,..,tm]
          depending on a ideal of constraints C in Q[t1,..,tm]. It detemines the different components as well as its type:
          \"Normal\", \"Special\", \"Accumulation\", \"Degenerate\". And it also classifies the \"Special\" components,
          determining the zero dimensional antiimage of the component and verifying if the component is a special
          hyper-surface of the family or not. It calls internally first grobcov and then locus with special options
          to obtain the complete result. The taxonomy that it provides, as well as the algorithms involved are
          described in the book: (not yet published)
          A. Montes. \"Discussing Parametric Polynomial Systems: The Groebner Cover\"

locusto(L); Transforms the output of locus, locusdg, envelop into a string that can be reed from different
          computational systems.

AssocTanToEnv(F,C,E); Having computed an envelop component E of a family of hyper-surfaces F, with constraints C,
          it returns the parameter values of the associated tangent hyper-surface of the family passing
          at one point of the envelop component E.

FamElemsAtEnvCompPoints(F,C,E) Having computed an envelop component E of a family of hyper-surfaces F,
          with constraints C, it returns the parameter values of all the hyper-surfaes of the family
          passing at one point of the envelop component E.

discrim(f,x); Determines the factorized discriminant of a degree 2 polynomial in the variable x. The polynomial
          can be defined on any ring where x is a variable. The polynomial f can depend on parameters and variables.

WLemma(F,A); Given an ideal F in K[a][x] and an ideal A in K[a], it returns the list (lpp,B,S)  were B is the reduced
          Groebner basis of the specialized F over the segment computed in P-representation (or optionally in
          C-representation). The basis is given by I-regular functions.

SEE ALSO: compregb_lib
";

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

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

// Development of the library:
// Library grobcov.lib
// (Groebner Cover):
// Release 0: (public)
// Initial data: 21-1-2008
// Uses buildtree for cgsdr
// Final data: 3-7-2008
// Release 2: (prived)
// Initial data: 6-9-2009
// Last version using buildtree for cgsdr
// Final data: 25-10-2011
// Release B: (prived)
// Initial data: 1-7-2012
// Uses both buildtree and KSW for cgsdr
// Final data: 4-9-2012
// Release G: (public)
// Initial data: 4-9-2012
// Uses KSW algorithm for cgsdr
// Final data: 21-11-2013
// Release L: (public)
// New routine ConsLevels: 25-1-2016
// Updated locus: 10-7-2015 (uses ConsLevels)
// Release M: (public)
// New routines for computing the envelop of a family of
//    hyper-surfaces and associated questions: 22-4-2016: 20-9-2016
// New routine WLemma for computing the result of
//    Wibmer's Lemma:  19-9-2016
// Final version October 2016

//*************Auxiliary routines**************

// 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 K[x1,..,xm]
static proc elimintfromideal(ideal J)
{
  int i;
  int j=0;
  ideal K;
  if (size(J)==0){return(ideal(0));}
  for (i=1;i<=ncols(J);i++){if (size(variables(J[i])) !=0){j++; K[j]=J[i];}}
  return(K);
}

// delfromideal: deletes the i-th polynomial from the ideal F
//    Works in any kind of ideal
static proc delfromideal(ideal F, int i)
{
  int j;
  ideal G;
  if (size(F)<i){ERROR("delfromideal was called with incorrect arguments");}
  if (size(F)<=1){return(ideal(0));}
  if (i==0){return(F)};
  for (j=1;j<=ncols(F);j++)
  {
    if (j!=i){G[ncols(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
//   Works in any kind of ideal
static proc delidfromid(ideal I, ideal J)
{
  int i; list r;
  ideal JJ=J;
  for (i=1;i<=size(I);i++)
  {
    r=memberpos(I[i],JJ);
    if (r[1])
    {
      JJ=delfromideal(JJ,r[2]);
    }
  }
  return(JJ);
}

// eliminates the ith element from a list or an intvec
static proc elimfromlist(l, int i)
{
  if(typeof(l)=="list"){list L;}
  if (typeof(l)=="intvec"){intvec L;}
  if (typeof(l)=="ideal"){ideal L;}
  int j;
  if((size(l)==0) or (size(l)==1 and i!=1)){return(l);}
  if (size(l)==1 and i==1){return(L);}
  // L=l[1];
  if(i>1)
  {
    for(j=1;j<=i-1;j++)
    {
      L[size(L)+1]=l[j];
    }
  }
  for(j=i+1;j<=size(l);j++)
  {
    L[size(L)+1]=l[j];
  }
  return(L);
}

// eliminates repeated elements form an ideal or matrix or module or intmat or bigintmat
static proc elimrepeated(F)
{
  int i;
  int nt;
  if (typeof(F)=="ideal"){nt=ncols(F);}
  else{nt=size(F);}

  def FF=F;
  FF=F[1];
  for (i=2;i<=nt;i++)
  {
    if (not(memberpos(F[i],FF)[1]))
    {
      FF[size(FF)+1]=F[i];
    }
  }
  return(FF);
}

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

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

// idcontains
// Input: ideal p, ideal q
// Output: 1 if p contains q,  0 otherwise
// If the routine is to be called from the top, a previous call to
// setglobalrings() is needed.
static proc idcontains(ideal p, ideal q)
{
  int t; int i;
  t=1; i=1;
  def P=p; def Q=q;
  attrib(P,"isSB",1);
  poly r;
  while ((t) and (i<=size(Q)))
  {
    r=reduce(Q[i],P);
    if (r!=0){t=0;}
    i++;
  }
  return(t);
}

// selectminideals
//   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
//   Works in Q[u_1,..,u_m]
static proc selectminideals(list L)
{
  list P; int i; int j; int t;
  if(size(L)==0){return(L)};
  if(size(L)==1){P[1]=1; return(P);}
  for (i=1;i<=size(L);i++)
  {
    t=1;
    j=1;
    while ((t) and (j<=size(L)))
    {
      if (i!=j)
      {
        if(idcontains(L[i],L[j])==1)
        {
          t=0;
        }
      }
      j++;
    }
    if (t){P[size(P)+1]=i;}
  }
  return(P);
}

// Auxiliary routine
// elimconstfac: eliminate the factors in the polynom f that are in Q[a]
// Input:
//   poly f:
//   list L: of components of the segment
// Output:
//   poly f2  where the factors of f in Q[a] that are non-null on any component
//   have been dropped from f
static proc elimconstfac(poly f)
{
  int cond; int i; int j; int t;
  if (f==0){return(f);}
  def RR=basering;
  setring(@R);
  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);
};

static proc memberpos(f,J)
//"USAGE:  memberpos(f,J);
//         (f,J) expected (polynomial,ideal)
//               or       (int,list(int))
//               or       (int,intvec)
//               or       (intvec,list(intvec))
//               or       (list(int),list(list(int)))
//               or       (ideal,list(ideal))
//               or       (list(intvec),  list(list(intvec))).
//         Works in any kind of ideals
//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);
//}

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

// Input:
//  A,B: lists of ideals
// Output:
//   1 if both lists of ideals are equal, or 0 if not
static proc equallistsofideals(list A, list B)
{
 int i;
 int tes=0;
 if (size(A)!=size(B)){return(tes);}
 tes=1; i=1;
 while(tes==1 and i<=size(A))
 {
   if (equalideals(A[i],B[i])==0){tes=0; return(tes);}
   i++;
 }
 return(tes);
}

// Input:
//  A,B:  lists of P-rep, i.e. of the form [p_i,[p_{i1},..,p_{ij_i}]]
// Output:
//   1 if both lists of P-reps are equal, or 0 if not
static proc equallistsA(list A, list B)
{
  int tes=0;
  if(equalideals(A[1],B[1])==0){return(tes);}
  tes=equallistsofideals(A[2],B[2]);
  return(tes);
}

// Input:
//  A,B:  lists lists of of P-rep, i.e. of the form [[p_1,[p_{11},..,p_{1j_1}]] .. [p_i,[p_{i1},..,p_{ij_i}]]
// Output:
//   1 if both lists of lists of P-rep are equal, or 0 if not
static proc equallistsAall(list A,list B)
{
 int i; int tes;
 if(size(A)!=size(B)){return(tes);}
 tes=1; i=1;
 while(tes and i<=size(A))
 {
   tes=equallistsA(A[i],B[i]);
   i++;
 }
 return(tes);
}

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

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

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

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

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

proc setglobalrings()
"USAGE:   setglobalrings();
          No arguments.
          Can be called when a parametric ideal Q[a][x] is in use. (a=parameters, x=variables).
RETURN: After its call the rings Grobcov::@R=Q[a][x], Grobcov::@P=Q[a],  Grobcov::@RP=Q[x,a] are defined as
          global variables. (a=parameters, x=variables).
NOTE: It is called internally by many basic routines of the library grobcov, cgsdr, extendGC, pdivi, pnormalf, locus,
          locusdg, envelop, WLemma, and killed before the output. The user does not need to call it.
          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 Q[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;
  export(@R);      // global ring Q[a][x]
  export(@P);      // global ring Q[a]
  export(@RP);     // global ring K[x,a] with product order
  setring(RR);
};
example
{
  "EXAMPLE:"; echo = 2;
  if(defined(R)){kill R;}
  ring R=(0,a,b),(x,y,z),dp;
  setglobalrings();

  R;

  Grobcov::@R;

  Grobcov::@P;

  Grobcov::@RP;

  ringlist(Grobcov::@R);

  ringlist(Grobcov::@P);

 ringlist(Grobcov::@RP);
}

// cld : clears denominators of an ideal and normalizes to content 1
//        can be used in @R or @P or @RP
// input:
//        ideal J (J can be also poly), but the output is an ideal;
// output:
//        ideal Jc (the new form of ideal J without denominators and
//        normalized to content 1)
static proc cld(ideal J)
{
  if (size(J)==0){return(ideal(0));}
  int te=0;
  def RR=basering;
  if(not(defined(@RP)))
  {
    te=1;
    setglobalrings();
  }
  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);
  if(te){kill @R; kill @RP; kill @P;}
  return(Jc);
};

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

// pdivi : pseudodivision of a parametric polynomial f by a parametric ideal F in Q[a][x].
// input:
//   poly  f
//   ideal F
// output:
//   list (poly r, ideal q, poly mu)
proc pdivi(poly f,ideal F)
"USAGE: pdivi(poly f,ideal F);
          poly f: the polynomial in Q[a][x] to be divided
          ideal F: the divisor ideal in Q[a][x].
          (a=parameters, x=variables).
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 Q[a][x]. Returns a list (r,q,m) such that
          m*f=r+sum(q.F),
          and no lpp of a divisor divides a pp of r.
KEYWORDS: division; reduce
EXAMPLE:  pdivi; shows an example"
{
  F=simplify(F,2);
  int i;
  int j;
  poly v=1;
  for(i=1;i<=nvars(basering);i++){v=v*var(i);}
  poly r=0;
  poly mu=1;
  def p=f;
  ideal q;
  for (i=1; i<=ncols(F); i++){q[i]=0;};
  ideal lpf;
  ideal lcf;
  for (i=1;i<=ncols(F);i++){lpf[i]=leadmonom(F[i]);}
  for (i=1;i<=ncols(F);i++){lcf[i]=leadcoef(F[i]);}
  poly lpp;
  poly lcp;
  poly qlm;
  poly nu;
  poly rho;
  int divoc=0;
  ideal qlc;
  while (p!=0)
  {
    i=1;
    divoc=0;
    lpp=leadmonom(p);
    lcp=leadcoef(p);
    while (divoc==0 and i<=size(F))
    {
      qlm=lpp/lpf[i];
      if (qlm!=0)
      {
        qlc=simpqcoeffs(lcp,lcf[i]);
        nu=qlc[2];
        mu=mu*nu;
        rho=qlc[1]*qlm;
        p=nu*p-rho*F[i];
        r=nu*r;
        for (j=1;j<=size(F);j++){q[j]=nu*q[j];}
        q[i]=q[i]+rho;
        divoc=1;
      }
      else {i++;}
    }
    if (divoc==0)
    {
      r=r+lcp*lpp;
      p=p-lcp*lpp;
    }
  }
  list res=r,q,mu;
  return(res);
}
example
{
  "EXAMPLE:"; echo = 2;
  if(defined(R)){kill R;}
  ring R=(0,a,b,c),(x,y),dp;
  short=0;
  // 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]);
}

// 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)
static proc pspol(poly f,poly g)
{
  def lcf=leadcoef(f);
  def lcg=leadcoef(g);
  def lpf=leadmonom(f);
  def lpg=leadmonom(g);
  def v=gcd(lpf,lpg);
  def s=simpqcoeffs(lcf,lcg);
  def vf=lpf/v;
  def vg=lpg/v;
  poly S=s[2]*vg*f-s[1]*vf*g;
  return(S);
}

// facvar: Returns all the free-square factors of the elements
//         of ideal J (non repeated). Integer factors are ignored,
//         even 0 is ignored. It can be called from ideal @R, but
//         the given ideal J must only contain poynomials in the
//         parameters.
//         Operates in the ring @P, but can be called from ring @R,
//         and the ideal @P must be defined calling first setglobalrings();
// input:  ideal J
// output: ideal Jc: Returns all the free-square factors of the elements
//         of ideal J (non repeated). Integer factors are ignored,
//         even 0 is ignored. It can be called from ideal @R.
static proc facvar(ideal J)
//"USAGE:   facvar(J);
//          J: an ideal in the parameters
//RETURN:   all the free-square factors of the elements
//          of ideal J (non repeated). Integer factors are ignored,
//          even 0 is ignored. It can be called from ideal @R, but
//          the given ideal J must only contain poynomials in the
//          parameters.
//NOTE:     Operates in the ring @P, and the ideal J must contain only
//          polynomials in the parameters, but can be called from ring @R.
//KEYWORDS: factor
//EXAMPLE:  facvar; shows an example"
{
  int i;
  def RR=basering;
  setring(@P);
  def Ja=imap(RR,J);
  if(size(Ja)==0){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
static 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(poly f,ideal E,ideal N);
          f: the polynomial in Q[a][x]  (a=parameters, x=variables) to be reduced modulo V(E) - V(N) of a segment in Q[a].
          E: the null conditions ideal in Q[a]
          N: the non-null conditions in Q[a]
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 Q[a].
KEYWORDS: division; pdivi; reduce
EXAMPLE: pnormalf; shows an example"
{
    def RR=basering;
    int te=0;
    if (defined(@P)){te=1;}
    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;
  if(defined(R)){kill R;}
  ring R=(0,a,b,c),(x,y),dp;
  short=0;
  poly f=(b^2-1)*x^3*y+(c^2-1)*x*y^2+(c^2*b-b)*x+(a-bc)*y;
  ideal p=c-1;
  ideal q=a-b;
  pnormalf(f,p,q);
}

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

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

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

// redgbn: given a minimal basis (gb reducing) it
// reduces each polynomial wrt to V(E) \ V(N)
static 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);
}

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

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

// postredgb: given a minimal basis (gb reducing) it
// reduces each polynomial wrt to the others
static 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 reduced Groebner basis called in @R and computed in @P
static proc gbR(ideal N)
{
  def RR=basering;
  setring(@P);
  def Np=imap(RR,N);
  option(redSB);
  Np=std(Np);
  setring(RR);
  return(imap(@P,Np));
}

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

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

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

// Auxiliary routine to define an order for ideals
// Returns 1 if the ideal a is shoud precede ideal b by sorting them in idbefid order
//             2 if the the contrary happen
//             0 if the order is not relevant
static proc idbefid(ideal a, ideal b)
{
  poly fa; poly fb; poly la; poly lb;
  int te=1; int i; int j;
  int na=size(a);
  int nb=size(b);
  int nm;
  if (na<=nb){nm=na;} else{nm=nb;}
  for (i=1;i<=nm; i++)
  {
    fa=a[i]; fb=b[i];
    while((fa!=0) or (fb!=0))
    {
      la=lead(fa);
      lb=lead(fb);
      fa=fa-la;
      fb=fb-lb;
      la=leadmonom(la);
      lb=leadmonom(lb);
      if(leadmonom(la+lb)!=la){return(1);}
      else{if(leadmonom(la+lb)!=lb){return(2);}}
    }
  }
  if(na<nb){return(1);}
  else
  {
    if(na>nb){return(2);}
    else{return(0);}
  }
}

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

// Crep
// Computes the C-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 list (a,b) of the canonical ideals
//    the Crep of V(N) \ V(M)
// Assumed to be called in the ring @R or the ring @P or a ring ring Q[a]
proc Crep(ideal N, ideal M)
"USAGE:  Crep(ideal N,ideal M);
           ideal N (null ideal) (not necessarily radical nor maximal) in Q[a]. (a=parameters, x=variables).
           ideal M (hole ideal) (not necessarily containing N) in Q[a].
           To be called in a ring Q[a][x] or a ring Q[a]. But the ideals can contain only the parameters in Q[a].
RETURN: The canonical C-representation [P,Q] of the locally closed set, formed by a pair of radical ideals with P
           included in Q, representing the set V(P) - V(Q) = V(N) - V(M)
 KEYWORDS: locally closed set; canoncial form
 EXAMPLE:  Crep; shows an example"
{
  int te;
  def RR=basering;
  if(defined(@P)){te=1;  setring(@P); ideal Np=imap(RR,N); ideal Mp=imap(RR,M);}
  else {te=0; def Np=N; def Mp=M;}
  def La=Crep0(Np,Mp);
  if(size(La)==0)
  {
    if(te==1) {setring(RR); list LL;}
    if(te==0){list LL;}
    return(LL);
  }
  else
  {
    if(te==1) {setring(RR); def LL=imap(@P,La);}
    if(te==0){def LL=La;}
  return(LL);
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  short=0;
  if(defined(R)){kill R;}
  ring R=0,(a,b,c),lp;
  ideal p=a*b;
  ideal q=a,b-2;;
  Crep(p,q);
}

// Crep0
// Computes the C-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 list (a,b) of the canonical ideals
//    the Crep0 of V(N) \ V(M)
// Assumed to be called in a ring Q[x] (example @P)
static proc Crep0(ideal N, ideal M)
{
  list l;
  ideal Np=std(N);
  if (equalideals(Np,ideal(1)))
  {
    l=ideal(1),ideal(1);
    return(l);
  }
  int i;
  list L;
  ideal Q=Np+M;
  ideal P=ideal(1);
  L=minGTZ(Np);
  //"T_Np="; Np;
  //"T_minGTZ(Np)="; L;
  for(i=1;i<=size(L);i++)
  {
    L[i]=std(L[i]);
    if(idcontains(L[i],Q)==0)
    {
      P=intersect(P,L[i]);
    }
  }
  P=std(P);
  Q=std(radical(Q+P));
  if(equalideals(P,Q)){return(l);}
  list T=P,Q;
  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 be called in the ring @R or the ring @P or a ring ring Q[a]
proc Prep(ideal N, ideal M)
 "USAGE: Prep(ideal N,ideal M);
           ideal N (null ideal) (not necessarily radical nor maximal) in Q[a]. (a=parameters, x=variables).
           ideal M (hole ideal) (not necessarily containing N) in Q[a].
           To be called in a ring Q[a][x] or a ring Q[a]. But the ideals can contain only the parameters in Q[a].
 RETURN: The canonical P-representation of the locally closed set V(N) - V(M)
           Output: [Comp_1, .. , Comp_s ] where
           Comp_i=[p_i,[p_i1,..,p_is_i]]
 KEYWORDS: locally closed set; canoncial form
 EXAMPLE:  Prep; shows an example"
{
  int te;
  def RR=basering;
  if(defined(@P))
  {
    te=1; setring(@P); ideal Np=imap(RR,N); ideal Mp=imap(RR,M);
  }
  else {te=0; def Np=N; def Mp=M;}
  def La=Prep0(Np,Mp);
  if(te==1) {setring(RR); def LL=imap(@P,La); }
  if(te==0){def LL=La;}
  return(LL);
}
example
{
  "EXAMPLE:"; echo = 2;
  short=0;
  if(defined(R)){kill R;}
  ring R=0,(a,b,c),lp;
  ideal p=a*b;;
  ideal q=a,b-1;
  Prep(p,q);
}

// Prep0
// 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 be called in a ring Q[x] (example @P)
static proc Prep0(ideal N, ideal M)
{
  int te;
  if (N[1]==1)
  {
    return(list(list(ideal(1),list(ideal(1)))));
  }
  int i; int j; list L0;
  list Ni=minGTZ(N);
  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]+M);
    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);
    }
  }
  //"T_before="; prep;
  if (size(prep)==0){prep=list(list(ideal(1),list(ideal(1))));}
  //"T_Prep="; prep;
  //def Lout=CompleteA(prep,prep);
  //"T_Lout="; Lout;
  return(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 be called in the ring @R of the ring @P or a ring Q[a]
proc PtoCrep(list L)
"USAGE: PtoCrep(list L)
          list L=  [ Comp_1, .. , Comp_s ] where
          list Comp_i=[p_i,[p_i1,..,p_is_i] ], is the P-representation of a locally closed set V(N) - V(M).
          To be called in a ring Q[a][x] or a ring Q[a]. But the ideals can contain only the parameters in Q[a].
 RETURN:The canonical C-representation [P,Q] of the locally closed set. A pair of radical ideals with P included in Q,
          representing the set V(P) - V(Q)
 KEYWORDS: locally closed set; canoncial form
 EXAMPLE:  PtoCrep; shows an example"
{
  int te;
  def RR=basering;
  if(defined(@P)){te=1; setring(@P); list Lp=imap(RR,L);}
  else {te=0; def Lp=L;}
  def La=PtoCrep0(Lp);
  if(te==1) {setring(RR); def LL=imap(@P,La);}
  if(te==0){def LL=La;}
  return(LL);
}
example
{
 "EXAMPLE:"; echo = 2;
  if(defined(R)){kill R;}
  ring R=0,(a,b,c),lp;
  short=0;
  ideal p=a*(a^2+b^2+c^2-25);
  ideal q=a*(a-3),b-4;
  def Cr=Crep(p,q);
  Cr;
  def L=Prep(p,q);
  L;
  def Cr1=PtoCrep(L);
  Cr1;
}

// PtoCrep0
// 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-representation of V(N) \ V(M) = V(ida) \ V(idb)
// Assumed to be called in a ring Q[x] (example @P)
static proc PtoCrep0(list L)
{
  int te=0;
  def Lp=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);
    //"T_Lp[i]="; Lp[i];
    N=Lp[i][1];
    Lb=Lp[i][2];
    //"T_Lb="; Lb;
    ida=intersect(ida,N);
    for(j=1;j<=size(Lb);j++)
    {
      idb=intersect(idb,Lb[j]);
    }
  }
  //idb=radical(idb);
  def La=list(ida,idb);
  return(La);
}

// input: F a parametric ideal in Q[a][x]
// output: a disjoint and reduced Groebner System.
//      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(ideal F);
          F: ideal in Q[a][x] (a=parameters, x=variables) to be discussed.
          Computes a disjoint, reduced Comprehensive Groebner System (CGS).
          cgsdr is the starting point of the fundamental routine grobcov.
          The basering R, must be of the form Q[a][x], (a=parameters, x=variables), and should be defined previously.
RETURN: Returns a list T describing a reduced and disjoint Comprehensive Groebner System (CGS). The output is
          a list of  (full,hole,basis), where the ideals full and hole represent the segment V(full) - V(hole).
          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],lpph],  ... ,
             [lpp, [num,basis,segment],...,[num,basis,segment],lpph] ].
          The bases are the reduced Groebner bases (after normalization) for each point of the corresponding segment.
          The third element lpph of each lpp segment is the lpp of the homogenized ideal used ideal in the CGS as a string,
          that is shown only when option  (\"can\",1) is used.
          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
OPTIONS: An option is a pair of arguments: string, integer. To modify the default options, pairs of arguments
         -option name, value- of valid options must be added to the call.
          Inside grobcov the default option is \"can\",1. It can be used also with option \"can\",0 but then the output
          is not the canonical Groebner Cover. There, it cannot be used with option \"can\",2.
          Ehen cgsdr is called directly, the options are \"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: (default is 1) the output segments are given as as difference of varieties.
          With option \"out\",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 efficient but is only the starting point for the computation
          of grobcov. When grobcov is computed, the call to cgsdr inside uses specific options that are more
          expensive (\"can\",0-1,\"out\",0).

KEYWORDS: CGS; disjoint; reduced; Comprehensive Groebner System
EXAMPLE:  cgsdr; shows an example"
{
  int te;
  def RR=basering;
  if(defined(@P)){te=1;}
  else{te=0; setglobalrings();}
  // INITIALIZING OPTIONS
  int i; int j;
  def E=ideal(0);
  def N=ideal(1);
  int comment=0;
  int can=2;
  int out=1;
  poly f;
  ideal B;
  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)
  {
    " "; string("Begin cgsdr with options: ",LL);
  }
  int ish;
  for (i=1;i<=size(B);i++){ish=ishomog(B[i]); if(ish==0){break;};}
  if (ish)
  {
    if(comment>0){" "; string("The given system is homogneous");}
    def GS=KSW(B,LL);
    //can=0;
  }
  else
  {
  // ACTING DEPENDING ON OPTIONS
  if(can==2)
  {
    // WITHOUT HOMOHGENIZING
    if(comment>0){" "; string("Option of cgsdr: do not homogenize");}
    def GS=KSW(B,LL);
    setglobalrings();
  }
  else
  {
    if(can==1)
    {
      // COMPUTING THE HOMOGOENIZED IDEAL
      if(comment>0){" "; string("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);
      if(comment>0){" ";string("Basis before computing its std basis="); B1;}
      B1=std(B1);
      if(comment>0){" ";string("Basis after computing its std basis="); B1;}
      setring(RR);
      def B2=imap(RP,B1);
    }
    else
    { // (can=0)
       if(comment>0){" "; string( "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>0){" "; string("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]));
          if (typeof(GS[i][7])=="ideal")
          { GS[i][7]=postredgb(mingb(GS[i][7]));}
        }
      }
    }
  }
  if(te==0){kill @P; kill @R; kill @RP;}
  return(GS);
}
example
{
  "EXAMPLE:"; echo = 2;
  // Casas conjecture for degree 4:
  if(defined(R)){kill R;}
  ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp;
  short=0;
  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
static proc grsegments(list T)
{
  int i;
  list L;
  list lpp;
  list lp;
  list ls;
  int n=size(T);
  lpp[1]=T[n][1];
  L[1]=list(lpp[1],list(list(T[n][2],T[n][3],T[n][4])));
  if (n>1)
  {
    for (i=1;i<=size(T)-1;i++)
    {
      lp=memberpos(T[n-i][1],lpp);
      if(lp[1]==1)
      {
        ls=L[lp[2]][2];
        ls[size(ls)+1]=list(T[n-i][2],T[n-i][3],T[n-i][4]);
        L[lp[2]][2]=ls;
      }
      else
      {
        lpp[size(lpp)+1]=T[n-i][1];
        L[size(L)+1]=list(T[n-i][1],list(list(T[n-i][2],T[n-i][3],T[n-i][4])));
      }
    }
  }
  return(L);
}

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

// Auxiliary routine
// 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),pos1),..,(q_k,(q_k1,..,q_kl_k),posk))
// posi=(i,j) position of the component
//        the list of segments to be added to the holes
static proc addpart(list H, list C)
{
  list Q; int i; int j; int k; int l; int t; int t1;
  Q=H; intvec notQ; list QQ; list addq;
  // @Q2=list of (i,j) positions of the components that have been aded to some hole of the maximal ideals
  //          plus those of the components added to the holes.
  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]))
        {
          // \\ @Q2[size(@Q2)+1]=C[j][3];
          t=0;
          for (k=1;k<=size(C[j][2]);k++)
          {
            t1=1;
            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];
              // \\ @Q2[size(@Q2)+1]=C[j][3];
            }
          }
          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;
      }
    }
    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));
}

// Auxiliary routine called by addpart to finish the modification of the holes of a primepart
// of the union by the addition of the segments that do not correspond to
// that part.
// Works on @P ring.
static proc addpartfine(list H, list C0)
{
  //"T_H="; H;
  int i; int j; int k; int te; intvec notQ; int l; list sel;
  intvec jtesC;
  if ((size(H)==1) and (equalideals(H[1],ideal(1)))){return(H);}
  if (size(C0)==0){return(H);}
  list newQ; list nQ; list Q; list nQ1; list Q0;
  def Q1=H;
  //Q1=sortlistideals(Q1,idbefid);
  def C=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)
          {
            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;
  }
  if(size(Q1)==0){Q1=ideal(1),ideal(1);}
  return(Q1);
}

// Auxiliary rutine for gcover
// Deciding if combine is needed
// input: list LCU=( (basis1, p_1, (p11,..p1s1)), .. (basisr, p_r, (pr1,..prsr))
// output: (tes); if tes==1 then combine is needed, else not.
static proc needcombine(list LCU,ideal N)
{
  //"Deciding if combine is needed";;
  ideal BB;
  int tes=0; int m=1; int j; int k; poly sp;
  while((tes==0) and (m<=size(LCU[1][1])))
  {
    j=1;
    while((tes==0) and (j<=size(LCU)))
    {
      k=1;
      while((tes==0) 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=1;}
        }
        k++;
      }
      j++;
    }
    if(tes){break;}
    m++;
  }
  return(tes);
}

// Auxiliary routine
// precombine
// input:  L: list of ideals (works in @P)
// output: F0: ideal of polys. F0[i] is a poly in the intersection of
//             all ideals in L except in the ith one, where it is not.
//             L=(p1,..,ps);  F0=(f1,..,fs);
//             F0[i] \in intersect_{j#i} p_i
static proc precombine(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);
}

// Auxiliary routine
// 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)
static proc combine(list L, ideal F)
{
  // ATTENTION REVISE AND USE Pci and F
  int i; poly f;
  f=0;
  for(i=1;i<=size(L);i++)
  {
    f=f+F[i]*L[i][2];
  }
//   f=elimconstfac(f);
  f=primepartZ(f);
  return(f);
}

// Central routine for grobcov: ideal F is assumed to be homogeneous
// 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 )
static 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;
  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;
    kill crep; list crep;
    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";
    crep=PtoCrep(prep);
    if(size(LCU)>1){tes=1;}
    else
    {
      tes=0;
      for(k=1;k<=size(B);k++){B[k]=pnormalf(B[k],crep[1],crep[2]);}
    }
    // tes=needcombine(LCU,N);
    // if(tes==1){" ";string("combine is needed for segment ",i);" ";}
    //crep=PtoCrep(prep);
    if(tes)
    {
      // combine is needed
      kill B; ideal B;
      for (j=1;j<=size(LCU);j++)
      {
        LL[j]=LCU[j][2];
      }
      FF=precombine(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]);
        }
        B[k]=combine(L,FF);
      }
    }
    for(j=1;j<=size(B);j++)
    {
      B[j]=pnormalf(B[j],crep[1],crep[2]);
    }
    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)
  {
    string("Time in LCUnion + combine = ",timer-start);
    if(comment>=2){string("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], (a=parameters, x=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(ideal F[,options]);
          F: ideal in Q[a][x] (a=parameters, x=variables) to be discussed.
          This is the fundamental routine of the library. It computes the Groebner Cover of a
          parametric ideal F in Q[a][x]. See
          A. Montes , M. Wibmer, \"Groebner Bases for Polynomial Systems with parameters\".
          JSC 45 (2010) 1391-1425.)
          or the not yet published book
          A. Montes. \"Discussing Parametric Polynomial Systems: The Groebner Cover\"
          The Groebner Cover of a parametric ideal F 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 reducedGroebner bases of the ideal on every point
          of S_i. The ideal F must be defined on a parametric ring Q[a][x] (a=parameters, x=variables).
RETURN:   The list  [ [ lpp_1,basis_1,segment_1],  ...,  [lpp_s,basis_s,segment_s] ]
          optionally  [ [ 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.
          With option (\"showhom\",1) the lpph will be shown: The lpph corresponds to the lpp of the homogenized ideal
          and is different for each segment. It is given as a string, and shown only for information. With the default option
          \"can\",1, the segments have different lpph.
          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 Union_i ( V(p_i) - ( Union_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.
OPTIONS: An option is a pair of arguments: string, integer. To modify the default options, pairs of arguments
          -option name, value- of valid options must be added to the call.
          \"null\",ideal E: The default is (\"null\",ideal(0)).
          \"nonnull\",ideal N: The default is (\"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 of the bases is
          computed (single polynomials, but not specializing to non-zero for every point of the segment.
          With option (\"ext\",1) the full representation of the bases is computed (possible sheaves) and sometimes a
          simpler result is obtained, but the computation is more time consuming.
          \"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.
          \"showhom\",0-1: The default is (\"showhom\",0). Setting \"showhom\",1 will output the set of lpp of the
          homogenized ideal of each segment as last element.
          One can give none or whatever of these options.
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=#;
  list Se;
  int out=0;
  int showhom=0;
  int hom;
  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];}
              else
              {
                if (L0[2*i-1]=="showhom"){showhom=L0[2*i];}
              }
            }
          }
        }
      }
    }
  }
  if(not((canop==0) or (canop==1)))
  {
    string("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;
  LL[15]="showhom";    LL[16]=showhom;
  if (comment>=1)
  {
    string("Begin grobcov with options: ",LL);
  }
  kill S;
  def S=gcover(F,LL);
  // NOW extendGC
  if(extop)
  {
    S=extendGC(S,LL);
  }
  // NOW repop and showhom
  list Si; list nS;
  for(i=1;i<=size(S);i++)
  {
    if(repop==0){Si=list(S[i][1],S[i][2],S[i][3]);}
    if(repop==1){Si=list(S[i][1],S[i][2],S[i][4]);}
    if(repop==2){Si=list(S[i][1],S[i][2],S[i][3],S[i][4]);}
    if(showhom==1){Si[size(Si)+1]=S[i][5];}
    nS[size(nS)+1]=Si;
  }
  S=nS;
  if (comment>=1)
  {
    string("Time in grobcov = ", timer-start);
    string("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:
  if(defined(R)){kill R;}
  ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp;
  short=0;
  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);
  // EXAMPLE
  // M. Rychlik robot;
  // Complexity and Applications of Parametric Algorithms of Computational Algebraic Geometry.;
  // In: Dynamics of Algorithms, R. de la Llave, L. Petzold and J. Lorenz eds.;
  // IMA Volumes in Mathematics and its Applications, Springer-Verlag 118: 1-29 (2000).;
  // (18. Mathematical robotics: Problem 4, two-arm robot)."
  if (defined(R)){kill R;}
  ring R=(0,a,b,l2,l3),(c3,s3,c1,s1), dp;
  short=0;
  ideal S12=a-l3*c3-l2*c1,b-l3*s3-l2*s1,c1^2+s1^2-1,c3^2+s3^2-1; S12;
  grobcov(S12);
}

// Auxiliary routine called by extendGC
// extendpoly
// 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 extendpoly(poly f, ideal idp, ideal idq)
"USAGE: extendGC(poly f,ideal p,ideal q);
          f is a polynomial in Q[a][x] in generic representation of an I-regular function F defined on the locally closed
          segment S=V(p) - V(q).
          p,q are ideals in Q[a], representing the Crep of segment S.
RETURN: the extended representation of F in S.
          It can consist of a single polynomial or a set of polynomials when needed.
NOTE: The basering R, must be of the form Q[a][x], (a=parameters,x=variables), and should be defined previously.
          The ideals must be defined on R.
KEYWORDS: Groebner cover; parametric ideal; locally closed set; parametric ideal; generic representation;
           full representation
EXAMPLE:  extendpoly; shows an example"
{
  int te=0;
  if(defined(@P)){te=0;}
  else{te=1; setglobalrings();}
  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)
    {
      if(te==1){kill @R; kill @P; kill @RP;}
      return(fout);
    }
    else
    {
      if(te==1){kill @R; kill @P; kill @RP;}
      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)
    {
      if(te==1){kill @R; kill @P; kill @RP;}
      return(fout);
    }
    else
    {
      if(te==1){kill @R; kill @P; kill @RP;}
      return(idout);
    }
  }
}
example
{
"EXAMPLE:"; echo = 2;
  if(defined(R)){kill R;}
  ring R=(0,a1,a2),(x),lp;
  short=0;
  poly f=(a1^2-4*a1+a2^2-a2)*x+(a1^4-16*a1+a2^3-4*a2);
  ideal p=a1*a2;
  ideal q=a2^2-a2,a1*a2,a1^2-4*a1;
  extendpoly(f,p,q);
  // EXAMPLE;
  if (defined(R)){kill R;}
  ring R=(0,a0,b0,c0,a1,b1,c1,a2,b2,c2),(x), dp;
  short=0;
  poly f=(b1*a2*c2-c1*a2*b2)*x+(-a1*c2^2+b1*b2*c2+c1*a2*c2-c1*b2^2);
  ideal p=
  (-a0*b1*c2+a0*c1*b2+b0*a1*c2-b0*c1*a2-c0*a1*b2+c0*b1*a2),
  (a1^2*c2^2-a1*b1*b2*c2-2*a1*c1*a2*c2+a1*c1*b2^2+b1^2*a2*c2-b1*c1*a2*b2+c1^2*a2^2),
  (a0*a1*c2^2-a0*b1*b2*c2-a0*c1*a2*c2+a0*c1*b2^2+b0*b1*a2*c2-b0*c1*a2*b2-c0*a1*a2*c2+c0*c1*a2^2),
  (a0^2*c2^2-a0*b0*b2*c2-2*a0*c0*a2*c2+a0*c0*b2^2+b0^2*a2*c2-b0*c0*a2*b2+c0^2*a2^2),
  (a0*a1*c1*c2-a0*b1^2*c2+a0*b1*c1*b2-a0*c1^2*a2+b0*a1*b1*c2-b0*a1*c1*b2-c0*a1^2*c2+c0*a1*c1*a2),
  (a0^2*c1*c2-a0*b0*b1*c2-a0*c0*a1*c2+a0*c0*b1*b2-a0*c0*c1*a2+b0^2*a1*c2-b0*c0*a1*b2+c0^2*a1*a2),
  (a0^2*c1^2-a0*b0*b1*c1-2*a0*c0*a1*c1+a0*c0*b1^2+b0^2*a1*c1-b0*c0*a1*b1+c0^2*a1^2),
  (2*a0*a1*b1*c1*c2-a0*a1*c1^2*b2-a0*b1^3*c2+a0*b1^2*c1*b2-a0*b1*c1^2*a2
  -b0*a1^2*c1*c2+b0*a1*b1^2*c2-b0*a1*b1*c1*b2+b0*a1*c1^2*a2-c0*a1^2*b1*c2+c0*a1^2*c1*b2);
  ideal q=
  (-a1*c2+c1*a2),
  (-a1*b2+b1*a2),
  (-a0*c2+c0*a2),
  (-a0*b2+b0*a2),
  (-a0*c1+c0*a1),
  (-a0*b1+b0*a1),
  (-a1*b1*c2+a1*c1*b2),
  (-a0*b1*c2+a0*c1*b2),
  (-a0*b0*c2+a0*c0*b2),
  (-a0*b0*c1+a0*c0*b1);
  extendpoly(f,p,q);
}

// if L is a list(ideal,ideal)  return 1 else returns 0;
static proc typeofCrep(L)
{
  if(typeof(L)!="list"){return(0);}
  if(size(L)!=2){return(0);}
  if((typeof(L[1])!="ideal") or (typeof(L[2])!="ideal")){return(0);}
  return(1);
}

// 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.
// Can be called from the top
proc extendGC(list GC)
"USAGE: extendGC(list GC);
          list GC must the grobcov of a parametric ideal computed with option \"rep\",2.
          It determines the full representation.
          The default option of grobcov provides the bases in generic representation (the I-regular functions forming
          the bases are then given by a single polynomial. They can specialize to zero for some points of
          the segments, but in general, it is sufficient for many pouposes. Nevertheless the I-regular functions allow a
          full representation given by a set of polynomials specializing to the value of the function (after normalization)
          or to zero, but at least one of the polynomials specializes to non-zero.
          The full representation can be obtained by computing the grobcov with option \"ext\",1. (The default option
          there is \"ext\",0). With option \"ext\",1 the computation can be much more time consuming, but the result
          can be simpler.
          Alternatively, one can compute the full representation of the bases after computing grobcov with the default
          option for \"ext\" and the option \"rep\",2, that outputs both the Prep and the Crep of the segments, and then
          call \"extendGC\" to its output.
RETURN: When calling extendGC(grobcov(S,\"rep\",2)) the result is of the form
          [ [ [lpp_1,basis_1,segment_1,lpph_1], ... ,[lpp_s,basis_s,segment_s,lpph_s]] ],
          where each function of the basis can be given by an ideal of representants.
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:  extendGC; shows an example"
{
  int te;
  if(defined(@P)){te=1;}
  else{setglobalrings();}
  list S=GC;
  ideal idp;
  ideal idq;
  int i; int j; int m; int s; int k;
  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";
    if(te==0){kill @R; kill @RP; kill @P;}
    return(S);
  }
  if(typeofCrep(S[1][3])){k=3;}
  else{if(typeofCrep(S[1][4])){k=4;};}
  if(k==0)
  {
    "Warning! extendGC make sense only when grobcov has been called with option 'rep',1 or 'rep',2";
    if(te==0){kill @R; kill @RP; kill @P;}
    return(S);
  }
  poly leadc;
  poly ext;
  list SS;
  // Now extendGC
  for (i=1;i<=size(S);i++)
  {
    m=size(S[i][2]);
     for (j=1;j<=m;j++)
    {
      idp=S[i][k][1];
      idq=S[i][k][2];
      if (size(idp)>0)
      {
        leadc=leadcoef(S[i][2][j]);
        kill ext;
        def ext=extendpoly(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;
          }
        }
      }
    }
  }
  if(te==0){kill @R; kill @RP; kill @P;}
  return(S);
}
example
{
  "EXAMPLE:"; echo = 2;
  if(defined(R)){kill R;}
  ring R=(0,a0,b0,c0,a1,b1,c1),(x), dp;
  short=0;
  ideal S=a0*x^2+b0*x+c0,
          a1*x^2+b1*x+c1;
  def GCS=grobcov(S,"rep",2);
  GCS;
  def FGC=extendGC(GCS,"rep",0);
  // Full representation=
  FGC;
}

// Auxiliary routine
// nonzerodivisor
// input:
//    poly g in Q[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
static proc nonzerodivisor(poly gr, list Pr)
{
  def RR=basering;
  setring(@P);
  def g=imap(RR,gr);
  def P=imap(RR,Pr);
  int i; int k;  list J; ideal F;
  def f=g;
  ideal Pi;
  for (i=1;i<=size(P);i++)
  {
    option(redSB);
    Pi=std(P[i]);
    //attrib(Pi,"isSB",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);
}

//Auxiliary routine
// 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.
static 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);
}

// Auxiliary routine
// monoms
// Input: A polynomial f
// Output: The list of leading terms
static 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);
}


// Auxiliary routine
// 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 procedures.
static proc findindexpolys(list coefs)
{
  int i; int j; intvec numdens;
  for(i=1;i<=size(coefs);i++)
  {
    numdens[i]=size(coefs[i]);
  }
  def RR=basering;
  setring(@P);
  def coefsp=imap(RR,coefs);
  ideal cof; list combpolys; intvec v; int te; list mp;
  for(i=1;i<=size(coefsp);i++)
  {
    cof=ideal(0);
    for(j=1;j<=size(coefsp[i]);j++)
    {
      cof[j]=factorize(coefsp[i][j],3);
    }
    coefsp[i]=cof;
  }
  for(j=1;j<=size(coefsp[1]);j++)
  {
    v[1]=j;
    te=1;
    for (i=2;i<=size(coefsp);i++)
    {
      mp=memberpos(coefsp[1][j],coefsp[i]);
      if(mp[1])
      {
        v[i]=mp[2];
      }
      else{v[i]=0;}
    }
    combpolys[j]=v;
  }
  combpolys=reform(combpolys,numdens);
  setring(RR);
  return(combpolys);
}

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

// Auxiliary routine
// 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
static proc selectregularfun(matrix CC, ideal NN, ideal MM)
{
  int numcombused;
  def RR=basering;
  setring(@P);
  def C=imap(RR,CC);
  def N=imap(RR,NN);
  def M=imap(RR,MM);
  if (ncols(C)==1){return(C);}

  int i; int j; int k; list c; intvec ci; intvec c0; intvec c1;
  list T; list T0; list T1; list LL; ideal N1;ideal M1; int te=0;
  for(i=1;i<=ncols(C);i++)
  {
    if((C[1,i]!=0) and (C[2,i]!=0))
    {
      if(c0==intvec(0)){c0[1]=i;}
      else{c0[size(c0)+1]=i;}
    }
  }
  def C1=submat(C,1..2,c0);
  for (i=1;i<=ncols(C1);i++)
  {
    c=comb(ncols(C1),i);
    for(j=1;j<=size(c);j++)
    {
      ci=c[j];
      numcombused++;
      if(i==1){N1=N+C1[2,j]; M1=M;}
      if(i>1)
      {
        kill c0; intvec c0 ; kill c1; intvec c1;
        c1=ci[size(ci)];
        for(k=1;k<size(ci);k++){c0[k]=ci[k];}
        T0=searchinlist(c0,LL);
        T1=searchinlist(c1,LL);
        N1=T0[1]+T1[1];
        M1=intersect(T0[2],T1[2]);
      }
      T=list(ci,PtoCrep0(Prep0(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));
}

// Auxiliary routine
// 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
static proc searchinlist(intvec c,list L)
{
  int i; list T;
  for(i=1;i<=size(L);i++)
  {
    if (L[i][1]==c)
    {
      T=L[i][2];
      break;
    }
  }
  return(T);
}

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

// Auxiliary routine
// 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).
static proc smsheaves(list C, list L)
{
  int i; int i0; intvec W;
  int nor; int norn;
  intvec p;
  int sp=size(L[1]); int j0=1;
  for (i=1;i<=sp;i++){p[i]=1;}
  while (p!=0)
  {
    i0=0; nor=0;
    for (i=1; i<=size(L); i++)
    {
      norn=numones(L[i],pos(p));
      if (nor<norn){nor=norn; i0=i;}
    }
    W[j0]=i0;
    j0++;
    p=actualize(p,L[i0]);
  }
  list LL;
  for (i=1;i<=size(W);i++)
  {
    LL[size(LL)+1]=C[W[i]];
  }
  return(LL);
}

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

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

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


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

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

// Auxiliary routine
// radicalmember
static 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);
}

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

// Auxiliary routine
// plusP
// 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
static 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));
}

// Auxiliary routine
// 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
static proc reform(list combpolys, intvec numdens)
{
  list combp0; list combp1; int i; int j; int k; int l; list rest; intvec notfree;
  list free; intvec free1; int te; intvec v;  intvec w;
  int nummonoms=size(combpolys[1]);
  for(i=1;i<=size(combpolys);i++)
  {
    if(memberpos(0,combpolys[i])[1])
    {
      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);
}

// Auxiliary routine
// 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
static proc precombint(list L)
{
  int i; int j; int tes;
  def RR=basering;
  setring(@P);
  list L0; list L1; list L2; list L3; ideal F;
  L0=imap(RR,L);
  L1[1]=L0[1]; L2[1]=L0[size(L0)];
  for (i=2;i<=size(L0)-1;i++)
  {
    L1[i]=intersect(L1[i-1],L0[i]);
    L2[i]=intersect(L2[i-1],L0[size(L0)-i+1]);
  }
  L3[1]=L2[size(L2)];
  for (i=2;i<=size(L0)-1;i++)
  {
    L3[i]=intersect(L1[i-1],L2[size(L0)-i]);
  }
  L3[size(L0)]=L1[size(L1)];
  for (i=1;i<=size(L3);i++)
  {
    option(redSB); L3[i]=std(L3[i]);
  }
  for (i=1;i<=size(L3);i++)
  {
    tes=1; j=0;
    while((tes) 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);
}

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

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

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

// Auxiliary routine
// MDBasis: Minimal Dickson Basis
static 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);
}

// Auxiliary routine
// primepartZ
static proc primepartZ(poly f);
{
  def cp=content(f);
  def fp=f/cp;
  return(fp);
}

// LCMLC
static 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))
static 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) 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]=="out"){out=L[2*i];}
          }
        }
      }
    }
  }
  if (comment>0){string("Begin KSW with null = ",E," nonnull = ",N);}
  def CG=KSW0(F,E,N,comment);
  if (comment>0)
  {
    string("Number of segments in KSW (total) = ",size(CG));
    string("Time in KSW = ",timer-start);
  }
  if(out==0)
  {
    CG=KSWtocgsdr(CG);
    CG=groupKSWsegments(CG);
    if (comment>0)
    {
      string("Number of lpp segments = ",size(CG));
      string("Time in KSW + group + Prep = ",timer-start);
    }
  }
  if(defined(@P)){kill @P; kill @R; kill @RP;}
  return(CG);
}

// Auxiliary routine
// sqf
static 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);
}

// Auxiliary routine
// 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:
static 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);
  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;}
    list KK; KK[1]=list(E0,N0,ideal(1));
    return(KK);
   }
  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];}
    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);
}

// Auxiliary routine
// KSWtocgsdr
static 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);
}

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

// Auxiliary routine
// groupKSWsegments
// input:  the list of vertices of KSW
// output: the same terminal vertices grouped by lpp
static 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 ConsLevels ***************************

static proc zeroone(int n)
{
  list L; list L2;
  intvec e; intvec e2; intvec e3;
  int j;
  if(n==1)
  {
    e[1]=0;
    L[1]=e;
    e[1]=1;
    L[2]=e;
    return(L);
  }
  if(n>1)
  {
    L=zeroone(n-1);
    for(j=1;j<=size(L);j++)
    {
      e2=L[j];
      e3=e2;
      e3[size(e3)+1]=0;
      L2[size(L2)+1]=e3;
      e3=e2;
      e3[size(e3)+1]=1;
      L2[size(L2)+1]=e3;
    }
  }
  return(L2);
}

// Auxiliary routine
// subsets: the list of subsets of (1,..n)
static proc subsets(int n)
{
  list L; list L1;
  int i; int j;
  L=zeroone(n);
  intvec e; intvec e1;
  for(i=1;i<=size(L);i++)
  {
    e=L[i];
    kill e1; intvec e1;
    for(j=1;j<=n;j++)
    {
      if(e[n+1-j]==1)
      {
        if(e1==intvec(0)){e1[1]=j;}
        else{e1[size(e1)+1]=j};
      }
    }
    L1[i]=e1;
  }
  return(L1);
}

// Input a list A of locally closed sets in C-rep
// Output a list B of a simplified list of A
static proc SimplifyUnion(list A)
{
  int i; int j;
  list L=A;
  int n=size(L);
  if(n<2){return(A);}
  intvec w;
  for(i=1;i<=size(L);i++)
  {
    for(j=1;j<=size(L);j++)
    {
      if(i != j)
      {
        if(equalideals(L[i][2],L[j][1])==1)
        {
          L[i][2]=L[j][2];
          w[size(w)+1]=j;
        }
      }
    }
  }
  if(size(w)>0)
  {
    for(i=1; i<=size(w);i++)
    {
      j=w[size(w)+1-i];
      L=elimfromlist(L, j);
    }
  }
  ideal T=ideal(1);
  intvec v;
  for(i=1;i<=size(L);i++)
  {
    if(equalideals(L[i][2],ideal(1)))
    {
      v[size(v)+1]=i;
      T=intersect(T,L[i][1]);
    }
  }
  if(size(v)>0)
  {
    for(i=1; i<=size(v);i++)
    {
      j=v[size(v)+1-i];
      L=elimfromlist(L, j);
    }
  }
  if(equalideals(T,ideal(1))==0){L[size(L)+1]=list(std(T),ideal(1))};
  return(L);
}

// input list A=[[p1,q1],...,[pn,qn]] :
//                    the list of segments of a constructible set S, where each [pi,qi] is given in C-representation
// output list [topA,C]
//       where topA is the closure of A
//                 C is the list of segments of the complement of A given in C-representation
static proc FirstLevel(list A)
{
  int n=size(A);
  list T=zeroone(n);
  ideal P; ideal Q;
  list Cb;  ideal Cc=1;
  int i; int j;
  intvec t;
  ideal topA=1;
  list C;
  for(i=1;i<=n;i++)
  {
    topA=intersect(topA,A[i][1]);
  }
  //topA=std(topA);
  for(i=2; i<=size(T);i++)
  {
    t=T[i];
    //"T_t="; t;
    P=0; Q=1;
    for(j=1;j<=n;j++)
    {
      if(t[n+1-j]==1)
      {
        P=P+A[j][2];
      }
      else
      {
        Q=intersect(Q,A[j][1]);
      }
    }
    Cb=Crep0(P,Q);
    //"T_Cb="; Cb;
    if(size(Cb)!=0)
    {
      if( Cb[1][1]<>1)
      {
        C[size(C)+1]=Cb;
      }
    }
  }
  if(size(C)>1){C=SimplifyUnion(C);}
  return(list(topA,C));
}

// Input:
// Output:
static proc ConstoPrep(list L)
{
  list L1;
  int i; int j;
  list aux;
  for(i=1;i<=size(L);i++)
  {
    aux=Prep(L[i][2][1],L[i][2][2]);
    L1[size(L1)+1]=list(L[i][1],aux);
  }
  return(L1);
}

// Input:
//     list A =  [[P1,Q1], .. [Pn,Qn]]
//                  A constructible set as union of locally closed sets represented by pairs of ideals
// Output:
//     list L =[p1,p2,p3,...,pk]
//        where pi is the ideal of the closure of level i alternatively of  A or its complement
//        Note: the levels of A are [p1,p2], [p3,p4], [p5,p6],...
//                 the levels of C are [p2,p3],[p4,p5], ...
//                 expressed in C-representation
//    Assumed to be called in the ring @R or the ring @P or a ring ring Q[a]
proc ConsLevels(list A0)
"USAGE: ConsLevels(list L);
         L=[[P1,Q1],...,[Ps,Qs]] is a list of lists of of pairs of ideals represening the constructible set
         S=V(P1)\ V(Q1) u ... u V(Ps)\ V(Qs).
         To be called in a ring Q[a][x] or a ring Q[a]. But the ideals can contain only the parameters in Q[a].
RETURN: The list of ideals [a1,a2,...,at] representing the closures of the canonical levels of S and its complement C
         wrt to the closure of S. The canonical levels of S are represented by theirs Crep. So we have:
         Levels of S:  [a1,a2],[a3,a4],...
         Levels of C:  [a2,a3],[a4,a5],...
         S=V(a1)\ V(a2) u V(a3)\ V(a4) u ...
         C=V(a2\ V(a3) u V(a4)\ V(a5) u ...
         It is used internally by locus procedure.
         The expression of S can be obtained from the output of ConsLevels by the call to ConsLevelsToLevels.
NOTE: Th ring can be Q[a][x] or Q[a], but the ideals can only contain parmeters in Q[a].
         The algorithm is described in
         J.M. Brunat, A. Montes. \"Computing the canonical representation of constructible sets.\"
         Math.  Comput. Sci. (2016) 19: 165-178.
KEYWORDS: constructible set; locally closed set; canonical form
EXAMPLE:  ConsLevels; shows an example"
{
  def RR=basering;
  int te;
  if(defined(@P)){te=1; setring(@P); list A=imap(RR,A0);}
  else {te=0; def A=A0;}

  list L; list C;
  list B; list T; int i;
  for(i=1; i<=size(A);i++)
  {
    T=Crep0(A[i][1],A[i][2]);
    B[size(B)+1]=T;
  }
  list K;
  while(size(B)>0)
  {
    K=FirstLevel(B);
    //"T_K="; K;
    L[size(L)+1]=K[1];
    B=K[2];
   }
  L[size(L)+1]=ideal(1);
  if(te==1) {setring(RR); def LL=imap(@P,L);}
  if(te==0){def LL=L;}
  return(LL);
}
example
{
"EXAMPLE:"; echo = 2;
  if(defined(R)){kill R;}
  ring R=0,(x,y,z),lp;
  short=0;
  ideal P1=x*(x^2+y^2+z^2-1);
  ideal Q1=z,x^2+y^2-1;
  ideal P2=y,x^2+z^2-1;
  ideal Q2=z*(z+1),y,x*(x+1);

  list Cr1=Crep(P1,Q1);
  list Cr2=Crep(P2,Q2);

  list L=list(Cr1,Cr2);
  L;
  ConsLevels(L);
}

// Converts the output of ConsLevels, given by the set of closures of the Levels of the constructible S
//     to an expression where the Levels are apparent.
// Input: The ouput of ConsLevels of the form
//    [A1,A2,..,Ak], where the Ai's are the closures of the levels.
// Output: An expression of the form
//      L1=[[1,[A1,A2]],[3,[A3,A4]],..,[2l-1,[A_{2l-1},A_{2l}]]] the list of Levels of S
proc ConsLevelsToLevels(list L)
"USAGE: ConsLevelsToLevels(list L);
          The input list L must be the output of the call to the routine ConsLevels of a constructible set:
          L=[a1,a2,..,ak], where the a's are the closures of the levels, determined by ConsLevels.
          ConsLevelsToLevels selects the levels of the constructible set.
          To be called in a ring Q[a][x] or a ring Q[a]. But the ideals can contain only the parameters in Q[a].
RETURN: The levels of the constructible set:
          Lc=[ [1,[a1,a2]],[3,[a3,a4]],..,[2l-1,[a_{2l-1},a_{2l}]] ] the list of  Levels of S
KEYWORDS: constructible sets; canonical form
EXAMPLE:  ConsLevelsToLevels shows an example"
{
  int n=size(L) div 2;
  int i;
  list L1; list L2;
  for(i=1; i<=n;i++)
  {
    L1[size(L1)+1]=list(2*i-1,list(L[2*i-1],L[2*i]));
  }
  return(L1);
}
example
{
"EXAMPLE:"; echo = 2;
  if(defined(R)){kill R;}
  ring R=0,(x,y,z),lp;
  short=0;
  ideal P1=(x^2+y^2+z^2-1);
  ideal Q1=z,x^2+y^2-1;
  ideal P2=y,x^2+z^2-1;
  ideal Q2=z*(z+1),y,x*(x+1);
  ideal P3=x;
  ideal Q3=5*z-4,5*y-3,x;

  list Cr1=Crep(P1,Q1);
  list Cr2=Crep(P2,Q2);
  list Cr3=Crep(P3,Q3);

  list L=list(Cr1,Cr2,Cr3);
  L;

  def LL=ConsLevels(L);
  LL;
  ConsLevelsToLevels(LL);
}

//**************************** End ConsLevels ******************

//******************** Begin locus ******************************

// indepparameters
// Auxiliary routine to detect 'Special' components of the locus
// Input: ideal B
// Output:
//   1 if the solutions of the ideal do not depend on the parameters
//   0 if they depend
static proc indepparameters(ideal B)
{
  def R=basering;
  ideal v=variables(B);
  setring @RP;
  def BP=imap(R,B);
  def vp=imap(R,v);
  ideal varpar=variables(BP);
  int te;
  te=equalideals(vp,varpar);
  setring(R);
  if(te){return(1);}
  else{return(0);}
}

// dimP0: Auxiliary routine
// if the dimension in @P of an ideal in the parameters has dimension 0 then it returns 0
// else it retuns 1
static proc dimP0(ideal N)
{
  def R=basering;
  setring(@P);
  int te=1;
  def NP=imap(R,N);
  attrib(NP,"IsSB",1);
  int d=dim(std(NP));
  if(d==0){te=0;}
  setring(R);
  return(te);
}

// Takes a list of intvec and sorts it and eliminates repeated elements.
// Auxiliary routine
static proc sortpairs(L)
{
  def L1=sort(L);
  def L2=elimrepeated(L1[1]);
  return(L2);
}

// Eliminates the pairs of L1 that are also in L2.
// Auxiliary routine
static proc minuselements(list L1,list L2)
{
  int i;
  list L3;
  for (i=1;i<=size(L1);i++)
  {
    if(not(memberpos(L1[i],L2)[1])){L3[size(L3)+1]=L1[i];}
  }
  return(L3);
}

// NorSing
// Input:
//   ideal B: the basis of a component of the grobcov
//   ideal E: the top of the component (assumed to be of dimension > 0 (single equation)
//   ideal N: the holes of the component
// Output:
//   int d: the dimension of B on the variables (antiimage).
//     if d>0    then the component is 'Normal'
//     if d==0 then the component is 'Singular'
static proc NorSing(ideal B, ideal E, ideal N, list #)
{
  int i; int j; int Fenv=0; int env; int dd;
  list DD=#;
  def RR=basering;
  int moverdim=2;
  int version=0;
  int nv=nvars(RR);
  if(nv<4){version=1;}
  int d;
  poly F;
  for(i=1;i<=(size(DD) div 2);i++)
  {
    if(DD[2*i-1]=="movdim"){moverdim=DD[2*i];}
    if(DD[2*i-1]=="version"){version=DD[2*i];}
    if(DD[2*i-1]=="family"){F=DD[2*i];}
  }
  if(F!=0){Fenv=1;}
  list L0;
  if(dimP0(E)==0){L0=2,"Normal";} // 2 es fals pero ha de ser >0 encara que sigui 0
  else
  {
    if(version==0)
    {
      // Computing std(B+E,plex(x,y,x1,..xn)) one can detect if there is a first part
      // independent of parameters giving the variables with dimension 0
     dd=indepparameters(B);
      if (dd==1){d=0; L0=d,string(B);} // ,determineF(B,F,E)
      else{d=1; L0=2,"Normal";}
    }
    else
    {
      def RH=ringlist(RR);
      //"T_RH="; RH;
      def H=RH;
      H[1]=0;
      H[2]=RH[1][2]+RH[2];
      int n=size(H[2]);
      intvec ll;
      for(i=1;i<=n;i++)
      {
        ll[i]=1;
      }
      H[3][1][1]="lp";
      H[3][1][2]=ll;
      def RRH=ring(H);
      setring(RRH);
      ideal BH=imap(RR,B);
      ideal EH=imap(RR,E);
      ideal NH=imap(RR,N);
      if(Fenv==1){poly FH=imap(RR,F);}
      for(i=1;i<=size(EH);i++){BH[size(BH)+1]=EH[i];}
      BH=std(BH);  // MOLT COSTOS!!!
      ideal G;
      ideal r; poly q;
      for(i=1;i<=size(BH);i++)
      {
        r=factorize(BH[i],1);
        q=1;
        for(j=1;j<=size(r);j++)
        {
          if((pdivi(r[j],NH)[1] != 0) or (equalideals(ideal(NH),ideal(1))))
          {
            q=q*r[j];
          }
        }
        if(q!=1){G[size(G)+1]=q;}
      }
      setring RR;
      def GG=imap(RRH,G);
      ideal GGG;
      if(defined(L0)){kill L0; list L0;}
      for(i=1;i<=size(GG);i++)
      {
        if(indepparameters(GG[i])){GGG[size(GGG)+1]=GG[i];}
      }
      GGG=std(GGG);
      ideal GLM;
      for(i=1;i<=size(GGG);i++)
      {
        GLM[i]=leadmonom(GGG[i]);
      }
      attrib(GLM,"IsSB",1);
      d=dim(std(GLM));
      string antiim=string(GGG);
      L0=d,antiim;
      if(d==0)
      {
        //" ";string("Antiimage of Special component = ", GGG);
      }
      else
      {
        L0[2]="Normal";
      }
    }
  }
  //"T_L0="; L0;
  return(L0);
}

static proc determineF(ideal A,poly F,ideal E)
{
  int env; int i;
  def RR=basering;
  def RH=ringlist(RR);
  def H=RH;
  H[1]=0;
  H[2]=RH[1][2]+RH[2];
  int n=size(H[2]);
  intvec ll;
  for(i=1;i<=n;i++)
  {
    ll[i]=1;
  }
  H[3][1][1]="lp";
  H[3][1][2]=ll;
  def RRH=ring(H);

        //" ";string("Antiimage of Special component = ", GGG);

   setring(RRH);
   list LL;
   def AA=imap(RR,A);
   def FH=imap(RR,F);
   def EH=imap(RR,E);
   ideal M=std(AA+FH);
   def rh=reduce(EH,M);
   //"T_AA="; AA; "T_FH="; FH; "T_EH="; EH; "T_rh="; rh;
   if(rh==0){env=1;} else{env=0;}
   setring RR;
          //L0[3]=env;
    //"T_env="; env;
    return(env);
}

//  DimPar
//  Auxilliary routine to NorSing determining the dimension of a parametric ideal
//  Does not use @P and define it directly because is assumes that
//                              variables and parameters have been inverted
 static proc DimPar(ideal E)
 {
   def RRH=basering;
   def RHx=ringlist(RRH);
   def Prin=ring(RHx[1]);
   setring(Prin);
   def E2=std(imap(RRH,E));
   def d=dim(E2);
   setring RRH;
   return(d);
 }

// locus0(G): Private routine used by locus (the public routine), that
//                builds the diferent components.
// input:      The output G of the grobcov (in generic representation, which is the default option)
//       Options: The algorithm allows the following options as pair of arguments:
//                "movdim", d  : by default movdim is 2 but it can be set to other values
//                    when locus is called by envelop then as default is uses d=dim @P
//                "version", v   :  There are two versions of the algorithm. ('version',1) is
//                 a full algorithm that always distinguishes correctly between 'Normal'
//                 and 'Special' components, whereas ('version',0) can decalre a component
//                 as 'Normal' being really 'Special', but is more effective. By default ('version',1)
//                 is used when the number of variables is less than 4 and 0 if not.
//                 The user can force to use one or other version, but it is not recommended.
//                 "system", ideal F: if the initial systrem is passed as an argument. This is actually not used.
//                 "comments", c: by default it is 0, but it can be set to 1.
//                 Usually locus problems have mover coordinates, variables and tracer coordinates.
//                 The mover coordinates are to be placed as the last variables, and by default,
//                 its number is dim @P, but as option it can be set to other values.
// output:
//         list, the canonical P-representation of the Normal and Non-Normal locus:
//              The Normal locus has two kind of components: Normal and Special.
//              The Non-normal locus has two kind of components: Accumulation and Degenerate.
//              This routine is compemented by locus that calls it in order to eliminate problems
//              with degenerate points of the mover.
//         The output components are given as
//              ((p1,(p11,..p1s_1),type_1,level_1),..,(pk,(pk1,..pks_k),type_k,level_k)
//         The components are given in canonical P-representation of the subset.
//              If all levels of a class of locus are 1, then the set is locally closed. Otherwise the level
//              gives the depth of the component.
static proc locus0(list GG, list #)
{
  int te=0;
  int t1=1; int t2=1; int i;
  def R=basering;
  def RL=ringlist(R);
  if(defined(@P)==1){te=1; kill @P; kill @R; kill @RP; }
  setglobalrings();
  //Options
  list DD=#;
  ideal vmov;
  int moverdim=size(ringlist(R)[1][2]);
  int dimpar=size(RL[1][2]);
  int dimvar=size(RL[2]);
  int nv=nvars(R);
  if(moverdim>nv){moverdim=nv;}
  for(i=1;i<=moverdim;i++)
  {
    //string("T_moverdim=",moverdim,"T_nv=",nv,"T_var(",i,")=",var(i));
    vmov[size(vmov)+1]=var(i+nv-moverdim);
  }
  int version=0;
  if(nv<4){version=1;}
  int comment=0;
  ideal Fm;
  poly F;
  for(i=1;i<=(size(DD) div 2);i++)
  {
    if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    if(DD[2*i-1]=="version"){version=DD[2*i];}
    if(DD[2*i-1]=="system"){Fm=DD[2*i];}
    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
    if(DD[2*i-1]=="family"){F=DD[2*i];}
  }
  list HHH;
  if (GG[1][1][1]==1 and GG[1][2][1]==1 and GG[1][3][1][1][1]==0 and GG[1][3][1][2][1]==1)
    {return(HHH);}
  list G1; list G2;
  def G=GG;
  list Q1; list Q2;
  int d; int j; int k;
  t1=1;
  for(i=1;i<=size(G);i++)
  {
    attrib(G[i][1],"IsSB",1);
    d=dim(std(G[i][1]));
    if(d==0){G1[size(G1)+1]=G[i];}
    else
    {
      if(d>0){G2[size(G2)+1]=G[i];}
    }
  }
  if(size(G1)==0){t1=0;}
  if(size(G2)==0){t2=0;}
  setring(@RP);
  if(t1)
  {
    list G1RP=imap(R,G1);
  }
  else {list G1RP;}
  list P1RP;
  ideal B;
  for(i=1;i<=size(G1RP);i++)
  {
    kill B;
     ideal B;
     for(k=1;k<=size(G1RP[i][3]);k++)
     {
       attrib(G1RP[i][3][k][1],"IsSB",1);
       G1RP[i][3][k][1]=std(G1RP[i][3][k][1]);
       for(j=1;j<=size(G1RP[i][2]);j++)
       {
         B[j]=reduce(G1RP[i][2][j],G1RP[i][3][k][1]);
       }
       P1RP[size(P1RP)+1]=list(G1RP[i][3][k][1],G1RP[i][3][k][2],B);
     }
  }  //In P1RP the basis has been reduced wrt the top
  setring(R);
  ideal h;
  list P1;
  if(t1)
  {
    P1=imap(@RP,P1RP);
    for(i=1;i<=size(P1);i++)
    {
      for(j=1;j<=size(P1[i][3]);j++)
      {
        h=factorize(P1[i][3][j],1);
        P1[i][3][j]=h[1];
        for(k=2;k<=size(h);k++)
        {
          P1[i][3][j]=P1[i][3][j]*h[k];
        }
      }
    }
  } //In P1 factors in the basis independent of the parameters have been obtained
  ideal BB; int dd; list NS;
  for(i=1;i<=size(P1);i++)
  {
     NS=NorSing(P1[i][3],P1[i][1],P1[i][2][1],DD);
     dd=NS[1];
     if(dd==0){P1[i][3]=NS;}  //"Special";
     else{P1[i][3]="Normal";}
  }
  list P2;
  for(i=1;i<=size(G2);i++)
  {
    for(k=1;k<=size(G2[i][3]);k++)
    {
      P2[size(P2)+1]=list(G2[i][3][k][1],G2[i][3][k][2]);
    }
  }
  list l;
  for(i=1;i<=size(P1);i++){Q1[i]=l; Q1[i][1]=P1[i];} P1=Q1;
  for(i=1;i<=size(P2);i++){Q2[i]=l; Q2[i][1]=P2[i];} P2=Q2;

  setring(@P);
  ideal J;
  if(t1==1)
  {
    def C1=imap(R,P1);
    def L1=AddLocus(C1);
   }
  else{list C1; list L1; kill P1; list P1;}
  if(t2==1)
  {
    def C2=imap(R,P2);
    def L2=AddLocus(C2);
  }
  else{list L2; list C2; kill P2; list P2;}
  for(i=1;i<=size(L2);i++)
  {
    J=std(L2[i][2]);
    d=dim(J);
    if(d+1==dimpar)
    {
      L2[i][4]=string("Degenerate",L2[i][4]);
    }
    else{L2[i][4]=string("Accumulation",L2[i][4]);}
  }
  list LN;
  if(t1==1)
  {
    for(i=1;i<=size(L1);i++){LN[size(LN)+1]=L1[i];}
  }
  if(t2==1)
  {
    for(i=1;i<=size(L2);i++){LN[size(LN)+1]=L2[i];}
  }
  int tLN=1;
  if(size(LN)==0){tLN=0;}
  setring(R);
  if(tLN)
  {
   def L=imap(@P,LN);
    for(i=1;i<=size(L);i++){if(size(L[i][2])==0){L[i][2]=ideal(1);}}
    list LL;
    for(i=1;i<=size(L);i++)
    {
      if(typeof(L[i][4])=="list") {L[i][4][1]="Special";}
      l[1]=L[i][2];
      l[2]=L[i][3];
      l[3]=L[i][4];
      l[4]=L[i][5];
      L[i]=l;
    }
  }
  else{list L;}
  return(L);
}

//  locus(G):  Special routine for determining the locus of points
//                 of  geometrical constructions.
//  input:      The output G of the grobcov (in generic representation, which is the default option for grobcov)
//  output:
//          list, the canonical P-representation of the Normal and Non-Normal locus:
//               The Normal locus has two kind of components: Normal and Special.
//               The Non-normal locus has two kind of components: Accumulation and Degenerate.
//          Normal component:
//           - the component has non-zero antiimage
//           - each point in the component has 0-dimensional antiimage.
//          Special component:
//           - Non-zero dimensional component whose antiimage is 0-dimensional.
//          Accumulation points:
//           - Zero-dimesnional  component whose antiimage is non-zero-dimensional.
//          Degenerate components:
//           - Non-zero-dimensional component
//           - each point in the component has non-zero-dimensional antiimage.
//          The output components are given as
//               ((p1,(p11,..p1s_1),type_1,level_1),..,(pk,(pk1,..pks_k),type_k,level_k)
//          The components are given in canonical P-representation of the subset.
//               If all levels of a class of locus are 1, then the set is locally closed. Otherwise the level
//               gives the depth of the component of the constructible set.
proc locus(list GG, list #)
"USAGE: locus(list G[,options])
          The list G must be the output of grobcov. Calling sequence: locus(grobcov(S))[,options]);
          The input must be the grobcov  of a parametrical ideal in Q[a][x], (a=parameters, x=variables). In practice
          a must be the tracer coordinates and x the mover coordinates and remaining auxiliary variables.  (Invert the
          concept of parameters and variables of the ring).
          Special routine for determining the locus of points of  geometrical constructions. Given a parametric ideal J
          representing the system determining the locus of points (a) which verify certain properties, the call to locus
          on the output of grobcov(J) determines the different classes of locus components, following the taxonomy
          defined in
          Abanades, Botana, Montes, Recio:
          \"An Algebraic Taxonomy for Locus Computation in Dynamic Geometry\".
          Computer-Aided Design 56 (2014) 22-33.
          The components can be \"Normal\", \"Special\", \"Accumulation\" or \"Degenerate\".
OPTIONS:An option is a pair of arguments: string, integer. To modify the default options, pairs of arguments
          -option name, value- of valid options must be added to the call.
          The algorithm allows the following options as pair of arguments:@*
          \"vmov\", ideal(mover variables)  : by default vmov are the last n variables, where n is the number of parameters
          of the ring (tracer variables plus extra parameters). Thus, if the mover coordinates are not indicated, locus
          algorithm will assume that they are the last n ring variables.
          When locus is called internally by envelop, by default, the mover variables are assumed to be all the ring variables.@*
          \"version\", v   :  There are two versions of the algorithm. (\"version\",1) is a full algorithm that always
          distinguishes correctly between \"Normal\" and \"Special\" components, whereas (\"version\",0) can declare a
          component to be \"Normal\" being in fact \"Special\", but it is more effective. By default,  (\"version\",1) is used
          when the number of variables is less than 4 and 0 if not. The user can force to use one or other version, but it is
          not recommended.@*
          \"comments\", c: by default it is 0, but it can be set to 1.
          Usually locus problems have mover coordinates, variables and tracer coordinates.
          Example of option call:
          locus(S,\"version\",1,\"vmov\",ideal(x1,y1))
RETURN:The output is a list of the components [C_1, .. , C_n] of the locus. Each component is given by
           Ci=[[pi,[pi1,..pi_s_i],type_i,level_i]]
          where the first element is the canonical P-representation of the subset. The type is one of
          \"Normal\", \"Special\", \"Accumulation\" or \"Degenerate\", and level is the depth of the segment in the
          constructible set of the locus. Generally it is 1, because the locus components are locally closed.
          The locus is divided into two class of subsets: the normal and the non-normal locus.
          The Normal locus has two kind of components: \"Normal\" and \"Special\".
          The Non-normal locus has two kind of components: \"Accumulation\" and \"Degenerate\".
          Normal component is n-1-dimensional component, where each point in the component has 0-dimensional
          antiimage, and the anti-image depends on the point in the component.
          Special component is n-1-dimensional component, where each point in the component has 0-dimensional
          antiimage, and the anti-image does not depend on the point in it. The Special components return more
          information, namely the antiimage of the component, that is 0-dimensional, and is independent of the point
          in the locus component.
          Accumulation component is of dimension less that n-1 (less than an hyper-surface)  whose anti-image is
          non-zero dimensional.
          Degenerate components is n-1-dimensional component, and each point in the component has non-zero-
          dimensional anti-image.
          The level is the depth of the segment of the constructible locus subset (normal and non-normal subsets).
          If all levels of a locus are 1, then all subsets are locally closed.
NOTE: The input must be the grobcov of the locus system in generic representation (\"ext\",0), which is the default.
KEYWORDS: geometrical locus; locus; dynamic geometry
EXAMPLE: locus; shows an example"
{
  int tes=0; int i;
  def R=basering;
  if(defined(@P)==1){tes=1; kill @P; kill @R; kill @RP;}
  setglobalrings();
  //Options
  list DD=#;
  ideal vmov;
  int moverdim=size(ringlist(R)[1][2]);
  int nv=nvars(R);
  if(moverdim>nv){moverdim=nv;}
  for(i=1;i<=moverdim;i++){vmov[size(vmov)+1]=var(i+nv-moverdim);}
  int version=0;
  if(nv<4){version=1;}
  int comment=0;
  ideal Fm;
  for(i=1;i<=(size(DD) div 2);i++)
  {
    if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    if(DD[2*i-1]=="version"){version=DD[2*i];}
    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
  }
  DD=list("vmov",vmov,"version",version,"comment",comment);
  int j; int k; int te;
  def B0=GG[1][2];
  def H0=GG[1][3][1][1];
  if (equalideals(B0,ideal(1)) )
  {
    return(locus0(GG,DD));
  }
  else
  {
    int n=nvars(R);
    ideal vB;
    ideal N;
    for(i=1;i<=size(B0);i++)
    {
      if(subset(variables(B0[i]),vmov)){N[size(N)+1]=B0[i];}
    }
    te=indepparameters(N);
    if(te)
    {
      string("locus detected that the mover must avoid points (",N,") in order to obtain the correct locus");" ";
      //eliminates segments of GG where N is contained in the basis
      list nGP;
      def GP=GG;
      ideal BP;
      for(j=1;j<=size(GP);j++)
      {
        te=1; k=1;
        BP=GP[j][2];
        while((te==1) and (k<=size(N)))
        {
          if(pdivi(N[k],BP)[1]!=0){te=0;}
          k++;
        }
        if(te==0){nGP[size(nGP)+1]=GP[j];}
      }
     }
    else
    {
      " ";string("Warning! Problem with more than one mover.");
      string("Try option 'vmov',ideal(of mover variables) to avoid some point of the mover");
      " ";"Elements of the basis of the generic segment in mover variables=";  N;" ";
      list L; return(L);
    }
  }
  def LL=locus0(nGP,DD);
  kill @RP; kill @P; kill @R;
  return(LL);
}
example
{
  "EXAMPLE:"; echo = 2;
  if(defined(R)){kill R;}
  ring R=(0,a,b),(x,y),dp;
  short=0;

  // Concoid
  ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
  S96;

  locus(grobcov(S96));
}

//  locusdg(G):  Special routine for determining the locus of points
//                 of  geometrical constructions in Dynamic Geometry.
//                 It is to be applied to the output of locus and selects
//                 as 'Relevant' the 'Normal' and the 'Accumulation'
//                 components.
//  input:      The output of locus(G);
//  output:
//          list, the canonical P-representation of the 'Relevant' components of the locus.
//          The output components are given as
//               ((p1,(p11,..p1s_1),type_1,level_1),..,(pk,(pk1,..pks_k),type_k,level_k)
//          The components are given in canonical P-representation of the subset.
//               If all levels of a class of locus are 1, then the set is locally closed. Otherwise the level
//               gives the depth of the component of the constructible set.
proc locusdg(list L)
"USAGE: locusdg(list L)
           Calling sequence:
           locusdg(locus(grobcov(S))).
RETURN: The output is the list of the \"Relevant\" components of the locus in Dynamic Geometry [C1,..,C:m], where
           C_i= [p_i,[p_i1,..p_is_i], \"Relevant\", level_i]
           The \"Relevant\" components are the \"Normal\" and \"Accumulation\" components of the locus. (See help for
           locus).
KEYWORDS: geometrical locus; locus; dynamic geometry
EXAMPLE: locusdg; shows an example"
{
  list LL;
  int i;
  for(i=1;i<=size(L);i++)
  {
    if(typeof(L[i][3])=="string")
    {
      if((L[i][3]=="Normal") or (L[i][3]=="Accumulation")){L[i][3]="Relevant"; LL[size(LL)+1]=L[i];}
    }
  }
  return(LL);
}
example
{
  "EXAMPLE:"; echo = 2;
  if(defined(R)){kill R;};
  ring R=(0,a,b),(x,y),dp;
  short=0;
  // Concoid
  ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
  S96;
  locus(grobcov(S96));
  locusdg(locus(grobcov(S96)));
}

// locusto: Transforms the output of locus, locusdg, envelop
//             into a string that can be reed from different computational systems.
// input:
//     list L: The output of locus or locusdg or envelop.
// output:
//     string s: Converts the input into a string readable by other programs
proc locusto(list L)
"USAGE: locusto(list L);
          The argument must be the output of locus or locusdg or envelop.
          It transforms the output into a string in standard form readable in other languages, not only Singular (Geogebra).
RETURN: The locus in string standard form
NOTE: It can only be called after computing either
           - locus(grobcov(F))                -> locusto( locus(grobcov(F)) )
           - locusdg(locus(grobcov(F)))  -> locusto( locusdg(locus(grobcov(F))) )
           - envelop(F,C)                       -> locusto( envelop(F,C) )
KEYWORDS: geometrical locus; locus; envelop; string
EXAMPLE:  locusto; shows an example"
{
  int i; int j; int k;
  string s="["; string sf="]"; string st=s+sf;
  if(size(L)==0){return(st);}
  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=string(s,",[");
      if(typeof(L[i][3])=="string")
      {
        s=string(s,string(L[i][3]),"]]");
      }
      else
      {
        for(k=1;k<=size(L[i][3]);k++)
        {
          s=string(s,"[",L[i][3][k],"],");
        }
        s[size(s)]="]";
        s=string(s,"]");
      }
    }
    if(size(L[i])>=4)
    {
      s[size(s)]=",";
      s=string(s,string(L[i][4]),"],");
    }
    s[size(s)]="]";
    s=string(s,",");
  }
  s[size(s)]="]";
  return(s);
}
example
{
  "EXAMPLE:"; echo = 2;
  if(defined(R)){kill R;}
  ring R=(0,x,y),(x1,y1),dp;
  short=0;
  ideal S=x1^2+y1^2-4,(y-2)*x1-x*y1+2*x,(x-x1)^2+(y-y1)^2-1;
  locusto(locus(grobcov(S)));
  locusto(locusdg(locus(grobcov(S))));
}

// envelop
// Input:
//   poly F: the polynomial defining the family of hypersurfaces in ring R=0,(x_1,..,x_n),(u_1,..,u_m),lp;
//   ideal C=g1,..,g_{n-1}:  the set of constraints;
//   options.
// Output: the components of the envolvent;
proc envelop(poly F, ideal C, list #)
"USAGE: envelop(poly F,ideal C[,options]);
          poly F must represent the family of hyper-surfaces for which on want to compute its envelop.
          ideal C must be the ideal of restrictions on the variables defining the family, and should contain less polynomials
          than the number of variables.
          (x_1,..,x_n) are the variables of the hyper-surfaces of F, that are considered as parameters of the parametric ring.
          (u_1,..,u_m) are the parameteres of the hyper-surfaces, that are considered as variables of the parametric ring.
          Calling sequence:
          ring R=(0,x_1,..,x_n),(u_1,..,u_m),lp;
          poly F=F(x_1,..,x_n,u_1,..,u_m);
          ideal C=g_1(u_1,..u_m),..,g_s(u_1,..u_m);
          envelop(F,C[,options]);   where s<m.
RETURN: The output is a list of the components [C_1, .. , C_n] of the locus. Each component is given by
           Ci=[pi,[pi1,..pi_s_i],type_i,level_i]
          where the first element is the canonical P-representation of the subset. The type is one of
          \"Normal\", \"Special\", \"Accumulation\" or \"Degenerate\", and level is the depth of the segment in the
          constructible set of the locus. Generally it is 1, because the locus components are locally closed.
          The locus is divided into two class of subsets: the normal and the non-normal locus.
          The Normal locus has two kind of components: \"Normal\" and \"Special\".
          The Non-normal locus has two kind of components: \"Accumulation\" and \"Degenerate\".
          Normal component is n-1-dimensional component, where each point in the component has 0-dimensional
          antiimage, and the anti-image depends on the point in the component.
          Special component is n-1-dimensional component, where each point in the component has 0-dimensional
          antiimage, and the anti-image does not depend on the point in it. The Special components return more
          information, namely the antiimage of the component, that is 0-dimensional, and is independent of the point
          in the locus component.
          Accumulation component is of dimension less that n-1 (less than an hyper-surface)  whose anti-image is
          non-zero dimensional.
          Degenerate components is n-1-dimensional component, and each point in the component has non-zero-
          dimensional anti-image.
          The level is the depth of the segment of the constructible locus subset (normal and non-normal subsets).
          If all levels of a locus are 1, then all subsets are locally closed.
OPTIONS: An option is a pair of arguments: string, integer. To modify the default options, pairs of arguments
         -option name, value- of valid options must be added to the call.
          The algorithm allows the following options as pair of arguments:
          \"vmov\", ideal(mover variables)  : by default vmov are  u_1,..,u_m. But it can be restricted by the user to the
          more convenient ones.
          \"version\", v   :  There are two versions of the algorithm. (\"version\",1) is a full algorithm that always
          distinguishes correctly between \"Normal\" and \"Special\" components, whereas (\"version\",0) can declare
          a component as \"Normal\" being really \"Special\", but is more effective. By default (\"version\",1) is used
          when the number of variables is less than 4 and 0 if not. The user can force to use one or other version, but
          it is not recommended.
          \"comments\", c: by default it is 0, but it can be set to 1.
NOTE: grobcov and locus are called internally.
          The basering R, must be of the form Q[a][x] (a=parameters, x=variables).
          This routine uses the generalized definition of envelop introduced in the book
           A. Montes. \"Discussing Parametric Polynomial Systems: The Groebner Cover\" not yet published.
KEYWORDS: geometrical locus; locus; envelop
EXAMPLE:  envelop; shows an example"
{
  def R=basering;
  int tes=0; int i;   int j;  int k; int m;
  int d;
  int dp;
  ideal BBB;
  if(defined(@P)==1){tes=1; kill @P; kill @R; kill @RP;}
  setglobalrings();
  //Options
  list DD=#;
  ideal vmov;
  int nv=nvars(R);
  for(i=1;i<=nv;i++){vmov[size(vmov)+1]=var(i);}
  int numpars=size(ringlist(R)[1][2]);
  int version=0;
  if(nv<4){version=1;}
  int comment=0;
  int familyinfo;
  ideal Fm;
  for(i=1;i<=(size(DD) div 2);i++)
  {
    if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    if(DD[2*i-1]=="version"){version=DD[2*i];}
    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
    if(DD[2*i-1]=="familyinfo"){familyinfo=DD[2*i];}
  };
  DD=list("vmov",vmov,"version",version,"comment",comment);
  int ng=size(C);
  ideal S=F;
  for(i=1;i<=size(C);i++){S[size(S)+1]=C[i];}
  int s=nv-ng;
  if(s>0)
  {
    matrix M[ng+1][ng+1];
    def cc=comb(nv,ng+1);
    poly J;
    for(k=1;k<=size(cc);k++)
    {
      for(j=1;j<=ng+1;j++)
      {
        M[1,j]=diff(F,var(cc[k][j]));
      }
      for(i=1;i<=ng;i++)
      {
        for(j=1;j<=ng+1;j++)
        {
          M[i+1,j]=diff(C[i],var(cc[k][j]));
        }
      }
      J=det(M);
      S[size(S)+1]=J;
    }
 }
 if(comment>0){"System S before grobcov ="; S;}
  def G=grobcov(S,DD);
  list HHH;
  if (G[1][1][1]==1 and G[1][2][1]==1 and G[1][3][1][1][1]==0 and G[1][3][1][2][1]==1)
  {return(HHH);}
   //DD[size(DD)+1]="vmov";
   //DD[size(DD)+1]=4;
  def L=locus(G,DD);
  list GL;
  ideal fam; ideal env;

   def Rx=ringlist(R);
   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=P+D;
  list LL;
  list NormalComp;
  ideal Gi;
  ideal BBBB;
  poly B0;
  if(familyinfo==1)
  {
    for(i=1;i<=size(L);i++)
    {
      if(typeof(L[i][3])=="string")
      {
        if(L[i][3]=="Normal")
        {
          NormalComp[size(NormalComp)+1]=L[i][1];
          Gi=S;
          Gi[size(Gi)+1]=L[i][1][1];
          //"T_grobcov(Gi)="; grobcov(Gi);
          kill HHH; list HHH;
          //"T_L[i]="; L[i];
          //d=DimPar(L[i][1]);
          if(defined(SL)){kill SL;}
          def SL=C;
          SL[size(SL)+1]=F;
          for(j=1;j<=size(L[i][1]);j++)
          {
            SL[size(SL)+1]=L[i][1][j];
          }
          setring RP;
          if(defined(BBBB)){kill BBBB;}
          ideal BBBB;
          if(defined(BB)){kill BB;}
          def BB=imap(R,SL);
          if(defined(B0)){kill B0;}
          poly B0;
          if(defined(LLL)){kill LLL;}
          def LLL=imap(R,L);
          //"T_BB="; BB;
          BB=std(BB);
           for(j=1;j<=size(BB);j++)
           {
             B0=reduce(BB[j],LLL[i][1]);
             if(not(B0==0)){BBBB[size(BBBB)+1]=B0;}
           }
          setring R;
          BBB=imap(RP,BBBB);
           L[i][5]=BBB;
        }
      }
    }
    LL[1]=L;
    LL[2]=NormalComp;
    list LLL; list LLLL;
    int t;
    for(k=1;k<=size(LL[2]);k++)
    {
      for(i=1;i<=size(G);i++)
      {
        j=1; t=0;
        while(t==0 and j<=size(G[i][3]))
        {
          //"T_LL[2][k]="; LL[2][k];
          //"T_G[i][3][j][1]="; G[i][3][j][1];
          if(equalideals(LL[2][k],G[i][3][j][1]))
          {
            LLL[size(LLL)+1]=list(k,i,j);
            t=1;
          }
          j++;
        }
      }
    }
    LL[3]=LLL;

    for(k=1;k<=size(LLL);k++)
    {
      for(m=k+1;m<=size(LLL);m++)
      {
        for(i=1;i<=size(G[LLL[k][2]][3][LLL[k][3]][2]);i++)
        {
          for(j=1;j<=size(G[LLL[m][2]][3][LLL[m][3]][2]);j++)
          {
            //string("T_(G[",LLL[k][2],"][3][",LLL[k][3],"][2][",i,"])=");  G[LLL[k][2]][3][LLL[k][3]][2][i];
            //string("T_(G[",LLL[m][2],"][3][",LLL[m][3],"][2][",j,"])=");  G[LLL[m][2]][3][LLL[m][3]][2][j];
            if(equalideals( G[LLL[k][2]][3][LLL[k][3]][2][i] ,G[LLL[m][2]][3][LLL[m][3]][2][j]))
            {
              //"T_GGG="; LL[2][LLL[k][1]]; LL[2][LLL[m][1]]; G[LLL[k][2]][3][LLL[k][3]][2][i];
              LLLL[size(LLLL)+1]=list(LL[2][LLL[k][1]],LL[2][LLL[m][1]], G[LLL[k][2]][3][LLL[k][3]][2][i]);
            }
          }
        }
      }
    }
    LL[3]=LLLL;
  }
  else{LL=L;}
  return(LL);
}
example
{
  "EXAMPLE:"; echo = 2;
  // Steiner Deltoid
  // 1. Consider the circle x1^2+y1^2-1=0, and a mover point M(x1,y1) on it.
  // 2. Consider the triangle A(0,1), B(-1,0), C(1,0).
  // 3. Consider lines passing through M perpendicular to two sides of ABC triangle.
  // 4. Obtain the envelop of the lines above.
  if(defined(R)){kill R;}
  ring R=(0,x,y),(x1,y1,x2,y2),lp;
  short=0;
  ideal C=(x1)^2+(y1)^2-1,
               x2+y2-1,
               x2-y2-x1+y1;
  matrix M[3][3]=x,y,1,x2,y2,1,x1,0,1;
  poly F=det(M);
  // Curves Family F
  F;
  // Conditions C=
  C;
  envelop(F,C);
}

proc AssocTanToEnv(poly F,ideal C, ideal E,list #)
"USAGE: AssocTanToEnv(poly F,ideal C,ideal E);
          poly F must be the family of hyper-surfaces whose envelop is analyzed.
          ideal C must be the ideal of restrictions on the variables for defining the family.
          Must contain less  polynomials than the dimension n of the space.
          ideal E must be a component of the envelop of (F,C), previously computed by envelop.
          (x_1,..,x_n) are the variables of the hypersurfaces of F, that are considered as parameters of the parametric ring.
          (u_1,..,u_m) are the parameteres of the hyper-surfaces, that are considered as variables of the parametric ring.
          Having computed an envelop component E of a family of hyper-surfaces F, with constraints C,
          it returns the parameter values of the associated tangent hyper-surface of the family passing
          at one point of the envelop component E.
          Calling sequence:
          ring R=(0,x_1,..,x_n),(u_1,..,u_m),lp;
          poly F=F(x_1,..,x_n,u_1,..,u_m);
          ideal C=g_1(u_1,..u_m),..,g_s(u_1,..u_m);
          poly E(x_1,..,x_n);
          AssocTanToEnv(F,C,E,[,options]);
RETURN: list [lpp,basis,segment]. The basis determines the associated tangent hyper-surface at a point of the envelop
          component E. The segment is given in Prep. See the not yet published book
          A. Montes. \"Discussing Parametric Polynomial Systems: The Groebner Cover\".
NOTE: grobcov is called internally. The basering R, must be of the form Q[a][x] (a=parameters, x=variables).
KEYWORDS: geometrical locus; locus; envelop; associated tangent
EXAMPLE:  AssocTanToEnv; shows an example"
{
  def R=basering;
  int tes=0; int i;   int j;  int k; int m;
  int d;
  int dp;
  ideal EE=E;
  int moreinfo=1;
  ideal BBB;
  if(defined(@P)==1){tes=1; kill @P; kill @R; kill @RP;}
  setglobalrings();
  //Options
  list DD=#;
  ideal vmov;
  int nv=nvars(R);
  for(i=1;i<=nv;i++){vmov[size(vmov)+1]=var(i);}
  int numpars=size(ringlist(R)[1][2]);
  int version=0;
  if(nv<4){version=1;}
  int comment=0;
  int familyinfo;
  ideal Fm;
  for(i=1;i<=(size(DD) div 2);i++)
  {
    if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    if(DD[2*i-1]=="version"){version=DD[2*i];}
    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
    if(DD[2*i-1]=="familyinfo"){familyinfo=DD[2*i];}
    if(DD[2*i-1]=="moreinfo"){moreinfo=DD[2*i];}
  };
  DD=list("vmov",vmov,"version",version,"comment",comment);
  int ng=size(C);
  ideal S=F;
  for(i=1;i<=size(C);i++){S[size(S)+1]=C[i];}
  int s=nv-ng;
  if(s>0)
  {
    matrix M[ng+1][ng+1];
    def cc=comb(nv,ng+1);
    poly J;
    for(k=1;k<=size(cc);k++)
    {
      for(j=1;j<=ng+1;j++)
      {
        M[1,j]=diff(F,var(cc[k][j]));
      }
      for(i=1;i<=ng;i++)
      {
        for(j=1;j<=ng+1;j++)
        {
          M[i+1,j]=diff(C[i],var(cc[k][j]));
        }
      }
      J=det(M);
      S[size(S)+1]=J;
    }
 }
 for(i=1;i<=size(EE);i++)
 {
   S[size(S)+1]=EE[i];
 }
 if(comment>0){"System S before grobcov ="; S;}
  def G=grobcov(S,DD);
  //"T_G=";G;
  list GG;
  for(i=2;i<=size(G);i++)
  {
    GG[size(GG)+1]=G[i];
  }
  G=GG;
  //"T_G=";G;
  if(moreinfo>0){return(G);}
  else
  {
    int t=0;
    ideal H;
    i=1;
    while(t==0 and i<=size(G))
    {
      //string("T_G[",i,"][3][1][1][1]="); G[i][3][1][1][1];
      //string("T_EE="); EE;
      if(G[i][3][1][1][1]==EE)
      {
         t=1;
         H=G[i][2];
      }
      i++;
    }
    return(H);
  }
  return(G);
}
example
{
  "EXAMPLE:"; echo = 2;
  if(defined(R)){kill R;}
  ring R=(0,y0,x,y),(t),dp;
  short=0;
  poly F=(x-5*t)^2+y^2-3^2*t^2;
  F;
  ideal C;
  C;

  def Env=envelop(F,C);
  Env;
  // E is a component of the envelop:
  ideal E=Env[1][1];
  E;
  def A=AssocTanToEnv(F,C,E);
  A;
  // The basis of the parameter values of the associated tangent component is
  A[1][2][1];
  // Thus t=(5/12)*y0 and the associated tangent component at (x0,y0) is
  subst(F,t,(5/12)*y0);
  // EXAMPLE
  if(defined(R)){kill R;}
  ring R=(0,x,y,z),(x1,y1,z1),dp;
  short=0;
  poly F=(x-x1)^2+(y-y1)^2+(z-z1)^2-1;
  ideal C=(x1)^2+(y1)^2-1;
  short=0;
  def Env=envelop(F,C); Env;
  def E=Env[1][1];
  AssocTanToEnv(F,C,E);
  def E1=Env[2][1];
  AssocTanToEnv(F,C,E1);
}

proc FamElemsAtEnvCompPoints(poly F,ideal C, ideal E)
"USAGE: FamElemsAtEnvCompPoints(poly F,ideal C,ideal E);
          poly F must be the family of hyper-surfaces whose envelop is analyzed.
          ideal C must be the ideal of restrictions on the variables for defining the family. Must contain less  polynomials
          than the dimension n of the space.
          ideal E must be a component of the envelop of (F,C), previously computed by envelop.
          (x_1,..,x_n) are the variables of the hypersurfaces of F, that are considered as parameters of the parametric ring.
          (u_1,..,u_m) are the parameteres of the hyper-surfaces, that are considered as variables of the parametric ring.
          Having computed an envelop component E of a family of hyper-surfaces F, with constraints C,
          it returns the parameter values of the set of all hyper-surfaces of the family passing
          at one point of the envelop component E.
          Calling sequence:
          ring R=(0,x_1,..,x_n),(u_1,..,u_m),lp;
          poly F=F(x_1,..,x_n,u_1,..,u_m);
          ideal C=g_1(u_1,..u_m),..,g_s(u_1,..u_m);
          poly E(x_1,..,x_n);
          FamElemsAtEnvCompPoints(F,C,E[,options]);
RETURN: list [lpp,basis,segment]. The basis determines the parameter values of the of hyper-surfaces that pass at a
          fixed point of the envelop component E. The lpp determines the dimension of the set. The segment is the
          component and is given in Prep.
          Fixing the values of (x_1,..,x_n) inside E, the basis allows to detemine the values of the parameters (u_1,..u_m),
          of the hyper-surfaces passing at a point of E. See the not yet published book
          A. Montes. \"Discussing Parametric Polynomial Systems: The Groebner Cover\"
NOTE: grobcov is called internally.
          The basering R, must be of the form Q[a][x] (a=parameters, x=variables).
KEYWORDS: geometrical locus; locus; envelop; associated tangent
EXAMPLE:  FamElemsAtEnvCompPoints; shows an example"
{
  int i;
  ideal S=C;
  S[size(S)+1]=F;
  for(i=1;i<=size(E);i++){S[size(S)+1]=E[i];}
  def G=grobcov(S);
  list GG;
  for(i=2; i<=size(G); i++)
  {
    GG[size(GG)+1]=G[i];
  }
  return(GG);
}
example
{
  "EXAMPLE:"; echo = 2;
  if(defined(R)){kill R;}
  ring R=(0,y0,x,y),(t),dp;
  short=0;
  poly F=(x-5*t)^2+y^2-3^2*t^2;
  F;
  ideal C;
  C;

  def Env=envelop(F,C);
  Env;

  // E is a component of the envelop:
  ideal E=Env[1][1];
  E;
  def A=AssocTanToEnv(F,C,E);
  A;
  // The basis of the parameter values of the associated tangent component is
  A[1][2][1];
  // Thus t=(5/12)*y0 the assocoated tangent family element at (x0,y0) is
  subst(F,t,(5/12)*y0);

  FamElemsAtEnvCompPoints(F,C,E);
  // Thus (12*t^2-5*y0)^2=0 and the unique circle of the family passing at (x0,y0) in E
  // is the associated   tangent circle:
  subst(F,t,(5/12)*y0);
}

// discrim
proc discrim(poly F0, poly x0)
"USAGE: discrim(f,x);
          poly f: the polynomial in Q[a][x] or Q[x] of degree 2 in x
          poly x: can be a variable or a parameter of the ring.
RETURN: the factorized discriminant of f wrt x for discussing its sign
KEYWORDS: second degree; solve
EXAMPLE:  discrim; shows an example"
{
  def RR=basering;
  int i;
  int te;
  int d;  int dd;
  if(size(ringlist(RR)[1])>0)
  {
    te=1;
    setglobalrings();
    setring @RP;
    poly F=imap(RR,F0);
    poly X=imap(RR,x0);
  }
  else
  {poly F=F0; poly X=x0;}
  matrix M=coef(F,X);
  d=deg(M[1,1]);
  if(d>2){"Degree is higher than 2. No discriminant"; setring RR; return();}
    poly dis=(M[2,2])^2-4*M[2,1]*M[2,3];
    def disp=factorize(dis,0);
    if(te==0){return(disp);}
    else
    {
      setring RR;
      def disp0=imap(@RP,disp);
      return(disp0);
    }
}
example
{
  "EXAMPLE:"; echo = 2;
  if(defined(R)){kill R;}
  ring R=(0,a,b,c),(x,y),dp;
  short=0;
  poly f=a*x^2*y+b*x*y+c*y;
  discrim(f,x);
}

// AddLocus: auxilliary routine for locus0 that computes the components of the constructible:
// Input:  the list of locally closed sets to be added, each with its type as third argument
//     L=[ [LC[11],..,LC[1k_1],.., [LC[r1],..,LC[rk_r] ] where
//            LC[1]=[p1,[p11,..,p1k],typ]
// Output:  the list of components of the constructible union of L, with the type of the corresponding top
//               and the level of the constructible
//     L4= [[v1,p1,[p11,..,p1l],typ_1,level]_1 ,.. [vs,ps,[ps1,..,psl],typ_s,level_s]
static proc AddLocus(list L)
{
  list L1; int i; int j;  list L2; list L3;
  list l1; list l2;
  intvec v;
  for(i=1; i<=size(L); i++)
  {
    for(j=1;j<=size(L[i]);j++)
    {
      l1[1]=L[i][j][1];
      l1[2]=L[i][j][2];
      l2[1]=l1[1];
      if(size(L[i][j])>2){l2[3]=L[i][j][3];}
      v[1]=i; v[2]=j;
      l2[2]=v;
      L1[size(L1)+1]=l1;
      L2[size(L2)+1]=l2;
    }
  }
  L3=LocusConsLevels(L1);
  list L4; int level;
  ideal p1; ideal pp1; int t; int k; int k0; string typ; list l4;
  for(i=1;i<=size(L3);i++)
  {
    level=L3[i][1];
    for(j=1;j<=size(L3[i][2]);j++)
    {
      p1=L3[i][2][j][1];
      t=1; k=1;
      while((t==1) and (k<=size(L2)))
      {
        pp1=L2[k][1];
        if(equalideals(p1,pp1)){t=0; k0=k;}
        k++;
      }
      if(t==0)
      {
        v=L2[k0][2];
        l4[1]=v; l4[2]=p1; l4[3]=L3[i][2][j][2];  l4[5]=level;
        if(size(L2[k0])>2){l4[4]=L2[k0][3];}
        L4[size(L4)+1]=l4;
      }
      else{"ERROR p1 NOT FOUND";}
    }
  }
  return(L4);
}

// Input L: list of components in P-rep to be added
//         [  [[p_1,[p_11,..,p_1,r1]],..[p_k,[p_k1,..,p_kr_k]]  ]
// Output:
//          list of lists of levels of the different locally closed sets of
//          the canonical P-rep of the constructible.
//          [  [level_1,[ [Comp_11,..Comp_1r_1] ] ], .. ,
//             [level_s,[ [Comp_s1,..Comp_sr_1] ]
//          ]
//          where level_i=i,   Comp_ij=[ p_i,[p_i1,..,p_it_i] ] is a prime component.
// LocusConsLevels: given a set of components of locally closed sets in P-representation, it builds the
//       canonical P-representation of the corresponding constructible set of its union,
//       including levels it they are.
static proc LocusConsLevels(list L)
{
  list Lc; list Sc;
  int i;
  for(i=1;i<=size(L);i++)
  {
    Sc=PtoCrep0(list(L[i]));
    Lc[size(Lc)+1]=Sc;
  }
  list S=ConsLevels(Lc);
  S=ConsLevelsToLevels(S);
  list Sout;
  list Lev;
  for(i=1;i<=size(S);i++)
  {
    Lev=list(S[i][1],Prep(S[i][2][1],S[i][2][2]));
    Sout[size(Sout)+1]=Lev;
  }
  return(Sout);
}

//********************* End locus ****************************

//********************* Begin WLemma **********************

// input ideal F in @R
//          ideal a in @R but only depending on parameters
//          F is a generating ideal in V(a);
// output:  ideal b in @R but depending only on parameters
//              ideal G=GBasis(F) in V(a) \ V(b)
proc WLemma(ideal F,ideal a, list #)
"USAGE: WLemma(F,A[,options]);
          The first argument ideal F in K[a][x];
          The second argument ideal A in K[a].
          Calling sequence:
          ring R=(0,x_1,..,x_n),(u_1,..,u_m),lp;
          ideal  F=f1(x_1,..,x_n,u_1,..,u_m),..,fs(x_1,..,x_n,u_1,..,u_m);
          ideal A=g_1(u_1,..u_m),..,g_s(u_1,..u_m);
          list # : Options
          Calling sequence:
          WLemma(F,A[,options]);
OPTIONS: either (\"rep\", 0) or (\"rep\",1) the representation of
          the resulting segment, by default is
          0 =P-representation, (default) but can be set to
          1=C-representation.
RETURN: list of [lpp,B,S] = [leading power product, basis, segment],
          being B the reduced Groebner Basis given by I-regular functions in full representation, of the specialized
          ideal F on the segment S given in P- or C-representation. It is the result of Wibmer's Lemma. See
          A. Montes , M. Wibmer, \"Groebner Bases for Polynomial Systems with parameters\".
          JSC 45 (2010) 1391-1425.)
          or the not yet published book
          A. Montes. \"Discussing Parametric Polynomial Systems: The Groebner Cover\"

NOTE: The basering R, must be of the form Q[a][x] (a=parameters, x=variables).
KEYWORDS: Wibmer's Lemma
EXAMPLE:  WLemma; shows an example"
{
  list L=#;
  int i; int j;
  def RR=basering;
  setglobalrings();
  setring(@RP);
  ideal FF=imap(RR,F);
  FF=std(FF);
  ideal AA=imap(RR,a);
  AA=std(AA);
  ideal FFa;
  poly r;
  for(i=1; i<=size(FF);i++)
  {
    r=reduce(FF[i],AA);
    if(r!=0){FFa[size(FFa)+1]=r;}
  }
  setring RR;
  ideal Fa=imap(@RP,FFa);
  ideal AAA=imap(@RP,AA);
  ideal lppFa;
  ideal lcFa;
  for(i=1;i<=size(Fa);i++)
  {
    lppFa[size(lppFa)+1]=leadmonom(Fa[i]);
    lcFa[size(lcFa)+1]=leadcoef(Fa[i]);
  }
  setring @RP;
  ideal lccr=imap(RR,lppFa);
  lccr=std(lccr);
  setring RR;
  ideal lcc=imap(@RP,lccr);
  list J; list Jx;
  ideal Jci;
  ideal Jxi;
  list B;
  for(i=1;i<=size(lcc);i++)
  {
    kill Jci; ideal Jci; kill Jxi; ideal Jxi;
    for(j=1;j<=size(Fa);j++)
    {
      if(lppFa[j]==lcc[i])
      {
        Jci[size(Jci)+1]=lcFa[j];
        Jxi[size(Jxi)+1]=Fa[j];
      }
    }
    J[size(J)+1]=Jci;
    B[size(B)+1]=Jxi;
  }
  setring @P;
  list Jp=imap(RR,J);
  ideal JL=product(Jp);
  setring(RR);
  def JLA=imap(@P,JL);
  list PR;
  if (size(L)>0)
  {
    if((L[1]=="rep") and (L[2]==1))
    {
      PR=Crep(AAA, JLA);
    }
    else
    {PR=Prep(AAA, JLA);}
  }
  else{PR=Prep(AAA, JLA);}
//  setring(RR);
 // list PRR=imap(@P,PR);
  return(list(lcc,B,PR));
}
example
{
  "EXAMPLE:"; echo = 2;
  if(defined(R)){kill R;}
  if(defined(R)){kill R;}
  ring R=(0,a,b,c,d,e,f),(x,y),lp;
  ideal F=a*x^2+b*x*y+c*y^2,d*x^2+e*x*y+f*y^2;
  ideal A=a*e-b*d;
  WLemma(F,A);
  WLemma(F,A,"rep",1);
}

//********************* End WLemma ************************