source: git/Singular/LIB/grobcov.lib @ ede2ad8

//
version="version grobcov.lib 4.0.2.0 Jul_2015 "; // $Id$
           // version L;  July_2015;
category="General purpose";
info="
LIBRARY:  grobcov.lib   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.
          (http://www-ma2.upc.edu/~montes/).

AUTHORS:  Antonio Montes , Hans Schoenemann.

OVERVIEW:
            In 2010, the library was designed to contain Montes-Wibmer's
            algorithms for compute the canonical Groebner Cover of a
            parametric ideal as described in the paper:

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

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

            grobcov uses a first algorithm cgsdr that outputs a disjoint
            reduced Comprehensive Groebner System with constant lpp.
            For this purpose, in this library, the implemented algorithm is
            Kapur-Sun-Wang algorithm, because it is the most efficient
            algorithm known for this purpose.

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

            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 constructuble. They will be described in a forthcoming paper:

             J.M. Brunat, A. Montes,
            \"Computing the canonical representation of constructible sets\".
            Submited to Mathematics in Computer Science. July 2015.

            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:

            Abanades, Botana, Montes, Recio:
            \''An Algebraic Taxonomy for Locus Computation in Dynamic Geometry\''.
            Computer-Aided Design 56 (2014) 22-33.

            Recently also routines for computing the envelop of a family
            of curves (enverlop, envelopdg), to be used in Dynamic Geometry,
            has been included and will be described in a forthcoming paper:

             Abanades, Botana, Montes, Recio:
             \''Envelops in Dynamic Geometry using the Groebner cover\''.


            This version was finished on 31/07/2015

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

Prep(N,M);  Computes the canonical P-representation of V(N) \ V(M).

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

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

ConsLevels(L); Given a list L of locally closed sets, it returns the canonical levels
             of the constructible set of the union of them, as well as the levels
             of the complement. It is described in the paper

             J.M. Brunat, A. Montes,
            \"Computing the canonical representation of constructible sets\".
            Submited to Mathematics in Computer Science. July 2015.

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.
              described in the paper:

             \"An Algebraic Taxonomy for Locus Computation in Dynamic Geometry\",
             M. Abanades, F. Botana, A. Montes, T. Recio,
             Computer-Aided Design 56 (2014) 22-33.

             The components can be 'Normal', 'Special', 'Accumulation', 'Degenerate'.
             The output are the components is 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 two
             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 curves
             F in Q[x,y][x_1,..xn] depending on a ideal of constraints C in Q[x_1,..,xn].
             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 curve 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
             will be described in a forthcoming paper:

             Abanades, Botana, Montes, Recio:
             \''Envelops in Dynamic Geometry using the Gr\"obner cover\''.

envelopdg(ev); Is a special routine to determine the 'Relevant' components
             of the envelop of a family of curves to be used in Dynamic Geometry.
             It must be called to the output of envelop(F,C).

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


SEE ALSO: compregb_lib
";

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

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

// Library grobcov.lib
// (Groebner Cover):
// Release 0: (public)
// Initial data: 21-1-2008
// Uses buildtree for cgsdr
// Final data: 3-7-2008
// Release 2: (private)
// Initial data: 6-9-2009
// Last version using buildtree for cgsdr
// Final data: 25-10-2011
// Release B: (private)
// 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: 10-7-2015
// Updated locus: 10-7-2015 (uses Conslevels)
// New routines for computing the envelop of a family of curves: 22-1-2015
// Final data: 22-7-2015

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

static proc setglobalrings()
// "USAGE:   setglobalrings();
//           No arguments
// 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, extend, pdivi, pnormalf, locus, locusdg,
//           envelop, envelopdg, and killed before the output.
//           The user does not need to call it, except when it is interested
//           in using some internal routine of the library that
//           uses these rings.
//           The basering R, must be of the form Q[a][x], (a=parameters,
//           x=variables), and should be defined previously.
// KEYWORDS: ring, rings
// EXAMPLE:  setglobalrings; shows an example"
{
  if (defined(@P))
  {
    kill @P; kill @R; kill @RP;
  }
  def RR=basering;
  def @R=basering;  // must be of the form 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;
//   ring R=(0,a,b),(x,y,z),dp;
//   setglobalrings();
//   " ";"R=";R;
//   " ";"Grobcov::@R=";Grobcov::@R;
//   " ";"Grobcov::@P=";Grobcov::@P;
//   " ";"Grobcov::@RP=";Grobcov::@RP;
//  " "; "ringlist(Grobcov::@R)=";  ringlist(Grobcov::@R);
//  " "; "ringlist(Grobcov::@P)=";  ringlist(Grobcov::@P);
//  " "; "ringlist(Grobcov::@RP)=";  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(f,F);
          poly f: the polynomialin Q[a][x] to be divided
          ideal F: the divisor ideal in Q[a][x].
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.G), 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;
  ring R=(0,a,b,c),(x,y),dp;
  poly f=(ab-ac)*xy+(ab)*x+(5c);
  // Divisor=";
  f;
  ideal F=ax+b,cy+a;
  // Dividends=";
  F;
  def r=pdivi(f,F);
  // (Remainder, quotients, factor)=";
  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(f,E,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;}
    else{setglobalrings();}
    setring(@RP);
    def fa=imap(RR,f);
    def Ea=imap(RR,E);
    def Na=imap(RR,N);
    option(redSB);
    Ea=std(Ea);
    def r=cld(reduce(fa,Ea));
    poly f1=r[1];
    f1=Ered(r[1],Ea,Na);
    setring(RR);
    def f2=imap(@RP,f1);
    if(te==0){kill @R; kill @RP; kill @P;}
    return(f2)
};
example
{
  "EXAMPLE:"; echo = 2;
  ring R=(0,a,b,c),(x,y),dp;
  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 E=(c-1);
  ideal N=a-b;

  pnormalf(f,E,N);
}

// 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
// Works on the ring @P
proc Crep(ideal N, ideal M)
 "USAGE:  Crep(N,M);
 Input: ideal N (null ideal) (not necessarily radical nor maximal)
           ideal M (hole ideal) (not necessarily containing N)
 RETURN:  The canonical C-representation of the locally closed set.
           [ P,Q ], a pair of radical ideals with P included in Q,
           representing the set V(P) \ V(Q) = V(N) \ V(M)
 NOTE:   Operates in a ring R=Q[a] (a=parameters)
 KEYWORDS: locally closed set, canoncial form
 EXAMPLE:  Crep; shows an example"
{
  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);
  for(i=1;i<=size(L);i++)
  {
    if(idcontains(L[i],Q)==0)
    {
      P=intersect(P,L[i]);
    }
  }
  P=std(P);
  Q=std(radical(Q+P));
  list T=P,Q;
  return(T);
}
example
{
  "EXAMPLE:"; echo = 2;
  if(defined(Grobcov::@P)){kill Grobcov::@R; kill Grobcov::@P; kill Grobcov::@RP;}
  ring R=0,(x,y,z),lp;
  short=0;
  ideal E=x*(x^2+y^2+z^2-25);
  ideal N=x*(x-3),y-4;
  def Cr=Crep(E,N);
  Cr;
  def L=Prep(E,N);
  L;
  def Cr1=PtoCrep(L);
  Cr1;
}

// 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
// Works on the ring @P
proc Prep(ideal N, ideal M)
 "USAGE:  Prep(N,M);
 Input: ideal N (null ideal) (not necessarily radical nor maximal)
           ideal M (hole ideal) (not necessarily containing N)
 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]]
 NOTE:   Operates in a ring R=Q[a]  (a=parameters)
 KEYWORDS: Locally closed set, Canoncial form
 EXAMPLE:  Prep; shows an example"
{
  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_abans="; 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);
}
example
{
  "EXAMPLE:"; echo = 2;
  if(defined(Grobcov::@P)){kill Grobcov::@R; kill Grobcov::@P; kill Grobcov::@RP;}
  short=0;
  ring R=0,(x,y,z),lp;
  ideal E=x*(x^2+y^2+z^2-25);
  ideal N=x*(x-3),y-4;
  Prep(E,N);
}

// 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
// Works on the ring @P
proc PtoCrep(list L)
"USAGE: PtoCrep(L)
 Input L: list  [ Comp_1, .. , Comp_s ] where
          Comp_i=[p_i,[p_i1,..,p_is_i] ], is
          the P-representation of a locally closed set V(N) \ V(M)
 RETURN:  The canonical C-representation of the locally closed set
          [ P,Q ], a pair of radical ideals with P included in Q,
          representing the set V(P) \ V(Q)
 NOTE: Operates in a ring R=Q[a]  (a=parameters)
 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(Grobcov::@P)){kill Grobcov::@R; kill Grobcov::@P; kill Grobcov::@RP;}
  short=0;
  ring R=0,(x,y,z),lp;

  // (P,Q) represents a locally closed set
  ideal P=x^3+x*y^2+x*z^2-25*x;
  ideal Q=y-4,x*z,x^2-3*x;

  // Now compute the P-representation=
  def L=Prep(P,Q);
  L;
  // Now compute the C-representation=
  def J=PtoCrep(L);
  J;
  // Now come back recomputing the P-represetation of the C-representation=
  Prep(J[1],J[2]);
}

// 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)
// Works in a ring Q[u_1,..,u_m] and is called on it.
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);
    N=Lp[i][1];
    ida=intersect(ida,N);
    Lb=Lp[i][2];
    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 rComprehensive Groebner System disjoint and reduced.
//      It uses Kapur-Sun-Wang algorithm, and with the options
//      can compute the homogenization before  (('can',0) or ( 'can',1))
//      and dehomogenize the result.
proc cgsdr(ideal F, list #)
"USAGE: cgsdr(F);
          F: ideal in Q[a][x] (a=parameters, x=variables) to be discussed.
          To compute a disjoint, reduced Comprehensive Groebner System (CGS).
          cgsdr is the starting point of the fundamental routine grobcov.
          Inside grobcov it is used with options 'can' set to 0,1 and
          not with options ('can',2).
          It is to be used if only a disjoint reduced CGS is required.

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

          Options:
            \"can\",0-1-2: The default value is \"can\",2. In this case no
                homogenization is done. With option (\"can\",0) the given
                basis is homogenized, and with option (\"can\",1) the
                whole given ideal is homogenized before computing the
                cgs and dehomogenized after.
                with option (\"can\",0) the homogenized basis is used
                with option (\"can\",1) the homogenized ideal is used
                with option (\"can\",2) the given basis is used
            \"null\",ideal E: The default is (\"null\",ideal(0)).
            \"nonnull\",ideal N: The default (\"nonnull\",ideal(1)).
                When options \"null\" and/or \"nonnull\" are given, then
                the parameter space is restricted to V(E) \ V(N).
            \"comment\",0-1: The default is (\"comment\",0). Setting (\"comment\",1)
                will provide information about the development of the
                computation.
            \"out\",0-1: 1 (default) the output segments are given as
                as difference of varieties.
                0: the output segments are given in P-representation
                and the segments grouped by lpp
                With options (\"can\",0) and (\"can\",1) the option (\"out\",1)
                is set to (\"out\",0) because it is not compatible.
          One can give none or whatever of these options.
          With the default options (\"can\",2,\"out\",1), only the
          Kapur-Sun-Wang algorithm is computed. This is very 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).
RETURN: Returns a list T describing a reduced and disjoint
          Comprehensive Groebner System (CGS),
          With option (\"out\",0)
          the segments are grouped by
          leading power products (lpp) of the reduced Groebner
          basis and given in P-representation.
          The returned list is of the form:
           (
             (lpp, (num,basis,segment),...,(num,basis,segment),lpp),
             ..,,
             (lpp, (num,basis,segment),...,(num,basis,segment),lpp)
           )
          The bases are the reduced Groebner bases (after normalization)
          for each point of the corresponding segment.

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

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

NOTE:  The basering R, must be of the form Q[a][x], (a=parameters,
          x=variables), and should be defined previously, and the ideal
          defined on R.
KEYWORDS: CGS, disjoint, reduced, Comprehensive Groebner System
EXAMPLE:  cgsdr; shows an example"
{
  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>=1){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);
      B1=std(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>=2){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:
  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; int used;
  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)
          {
            used++;
            if ((size(notQ)==1) and (notQ[1]==0)){notQ[1]=i;}
            else{notQ[size(notQ)+1]=i;}
            kill newQ; list newQ;
            for(k=1;k<=size(C[j][2]);k++)
            {
              nQ=minGTZ(Q[i]+C[j][2][k]);
              for(l=1;l<=size(nQ);l++)
              {
                option(redSB);
                nQ[l]=std(nQ[l]);
                newQ[size(newQ)+1]=nQ[l];
              }
            }
            sel=selectminideals(newQ);
            kill nQ1; list nQ1;
            for(l=1;l<=size(sel);l++)
            {
              nQ1[l]=newQ[sel[l]];
            }
            newQ=nQ1;
            for(l=1;l<=size(newQ);l++)
            {
              Q[size(Q)+1]=newQ[l];
            }
            break;
          }
        }
      }
    }
    kill Q1; list Q1;
    for(i=1;i<=size(Q);i++)
    {
      if(memberpos(i,notQ)[1]==0)
      {
        Q1[size(Q1)+1]=Q[i];
      }
    }
    //"T_Q1_0="; Q1;
    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);}
  //"T_Q1_1="; Q1;
  //if(used>0){string("addpartfine was ", used, " times used");}
  return(Q1);
}


// Auxiliary 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;
  ideal BB;
  for (i=1;i<=size(GP);i++)
  {
    kill LL;
    list LL;
    lpp=GP[i][1];
    GPi2=GP[i][2];
    lpph=GP[i][3];
    kill pairspP; list pairspP;
    for(j=1;j<=size(GPi2);j++)
    {
      pairspP[size(pairspP)+1]=GPi2[j][3];
    }
    LCU=LCUnion(pairspP);
    kill prep; list prep;
    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";
    kill BB;
    ideal BB;
    tes=1; m=1;
    while((tes) and (m<=size(LCU[1][1])))
    {
      j=1;
      while((tes) and (j<=size(LCU)))
      {
        k=1;
        while((tes) and (k<=size(LCU)))
        {
          if(j!=k)
          {
            sp=pnormalf(pspol(LCU[j][1][m],LCU[k][1][m]),LCU[k][2],N);
            if(sp!=0){tes=0;}
          }
          k++;
        }        //setglobalrings();
        if(tes)
        {
          BB[m]=LCU[j][1][m];
        }
        j++;
      }
      if(tes==0){break;}
      m++;
    }
    crep=PtoCrep(prep);
    if(tes==0)
    {
      // combine is needed
      kill B; ideal B;
      for (j=1;j<=size(LCU);j++)
      {
        LL[j]=LCU[j][2];
      }
      if (size(LCU)>1)
      {
        FF=precombint(LL);
      }
      for (k=1;k<=size(lpp);k++)
      {
        kill L; list L;
        for (j=1;j<=size(LCU);j++)
        {
          L[j]=list(LCU[j][2],LCU[j][1][k]);
        }
        if (size(LCU)>1)
        {
          B[k]=combine(L,FF);
        }
        else{B[k]=L[1][2];}
      }
    }
    else{B=BB;}
    for(j=1;j<=size(B);j++)
    {
      B[j]=pnormalf(B[j],crep[1],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\"i s 1 by default and can be set to 0 if we do not
//            need to obtain the canonical GC, but only a GC.
//            The value of \"ext\" is 0 by default and so the generic representation
//             of the bases is given. It can be set to 1, and then the full
//             representation of the bases is given.
//            The value of \"rep\" is 0 by default, and then the segments
//            are given in canonical P-representation. It can be set to 1
//            and then they are given in canonical C-representation.
//            If it is set to 2, then both representations are given.
// output:
//    list S: ((lpp,basis,(idp_1,(idp_11,..,idp_1s_1))), ..
//             (lpp,basis,(idp_r,(idp_r1,..,idp_rs_r))) ) where
//            each element of S corresponds to a lpp-segment
//            given by the lpp, the basis, and the P-representation of the segment
proc grobcov(ideal F,list #)
"USAGE:   grobcov(F); This is the fundamental routine of the
          library. It computes the Groebner cover of a parametric ideal. See
                Montes A., Wibmer M.,
               \"Groebner Bases for Polynomial Systems with parameters\".
               JSC 45 (2010) 1391-1425.)
          The Groebner Cover of a parametric ideal consist of a set of
          pairs(S_i,B_i), where the S_i are disjoint locally closed
          segments of the parameter space, and the B_i are the reduced
          Groebner bases of the ideal on every point of S_i.

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

          Options:
            \"null\",ideal E: The default is (\"null\",ideal(0)).
            \"nonnull\",ideal N: The default 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.
          One can give none or whatever of these options.
RETURN:   The list
          (
           (lpp_1,basis_1,segment_1,lpph_1),
           ...
           (lpp_s,basis_s,segment_s,lpph_s)
          )

          The lpp are constant over a segment and correspond to the
          set of lpp of the reduced Groebner basis for each point
          of the segment.
          The lpph corresponds to the lpp of the homogenized ideal
          and is different for each segment. It is given as a string,
          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.

NOTE: The basering R, must be of the form Q[a][x], (a=parameters,
          x=variables), and should be defined previously. The ideal must
          be defined on R.
KEYWORDS: Groebner cover, parametric ideal, canonical, discussion of
          parametric ideal.
EXAMPLE:  grobcov; shows an example"
{
  list S; int i; int ish=1; list GBR; list BR; int j; int k;
  ideal idp; ideal idq; int s; ideal ext; list SS;
  ideal E; ideal N; int canop;  int extop; int repop;
  int comment=0; int m;
  def RR=basering;
  setglobalrings();
  list L0=#;
  int out=0;
  L0[size(L0)+1]="res"; L0[size(L0)+1]=ideal(1);
  // default options
  int start=timer;
  E=ideal(0);
  N=ideal(1);
  canop=1; // canop=0 for homogenizing the basis but not the ideal (not canonical)
           // canop=1 for working with the homogenized ideal
  repop=0; // repop=0 for representing the segments in Prep
           // repop=1 for representing the segments in Crep
           // repop=2 for representing the segments in Prep and Crep
  extop=0; // extop=0 if only generic representation of the bases are to be computed
           // extop=1 if the full representation of the bases are to be computed
  for(i=1;i<=size(L0) div 2;i++)
  {
    if(L0[2*i-1]=="can"){canop=L0[2*i];}
    else
    {
      if(L0[2*i-1]=="ext"){extop=L0[2*i];}
      else
      {
        if(L0[2*i-1]=="rep"){repop=L0[2*i];}
        else
        {
          if(L0[2*i-1]=="null"){E=L0[2*i];}
          else
          {
            if(L0[2*i-1]=="nonnull"){N=L0[2*i];}
            else
            {
              if (L0[2*i-1]=="comment"){comment=L0[2*i];}
            }
          }
        }
      }
    }
  }
  if(not((canop==0) or (canop==1)))
  {
    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;
  if (comment>=1)
  {
    string("Begin grobcov with options: ",LL);
  }
  kill S;
  def S=gcover(F,LL);
  // NOW extend
  if(extop)
  {
    S=extend(S,LL);
  }
  else
  {
    // NOW representation of the segments by option repop
    list Si; list nS;
    if(repop==0)
    {
      for(i=1;i<=size(S);i++)
      {
        Si=list(S[i][1],S[i][2],S[i][3],S[i][5]);
        nS[size(nS)+1]=Si;
      }
      kill S;
      def S=nS;
    }
    else
    {
      if(repop==1)
      {
        for(i=1;i<=size(S);i++)
        {
          Si=list(S[i][1],S[i][2],S[i][4],S[i][5]);
          nS[size(nS)+1]=Si;
        }
        kill S;
        def S=nS;
      }
      else
      {
        for(i=1;i<=size(S);i++)
        {
          Si=list(S[i][1],S[i][2],S[i][3],S[i][4],S[i][5]);
          nS[size(nS)+1]=Si;
        }
        kill S;
        def S=nS;
      }
    }
  }
  if (comment>=1)
  {
    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:
  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);
}

// 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 extend(list GC, list #);
"USAGE: extend(GC); The default option of grobcov provides
          the bases in generic representation (the I-regular functions
          of the bases ara given by a single polynomial. It 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 bey 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 is \"ext\",0.
          With option \"ext\",1 the computation can be much more
          time consuming, even if the result can be simpler.
          Alternatively, one can compute the full representation of the
          bases after computing grobcov with the defaoult option \"ext\",0
          and the option \"rep\",2, that outputs both the Prep and the Crep
          of the segments and then call \"extend\" to the output.

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

    }
  }
  if(comment>=1){string("Time in extend = ",timer-start3);}
  if(te==0){kill @R; kill @RP; kill @P;}
  return(S);
}
example
{
  "EXAMPLE:"; echo = 2;
  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=extend(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
// deltai
// input:
//   int i:
//   list LPr: (p1,..,pr) of prime components of an ideal in Q[a]
// output:
//   list (fr,fnr) of two polynomials that are equal on V(pi)
//       and fr=0 on V(P) \ V(pi), and fnr is nonzero on V(pj) for all j.
static proc deltai(int i, list LPr)
{
  def RR=basering;
  setring(@P);
  def LP=imap(RR,LPr);
  int j; poly p;
  def F=ideal(1);
  poly f;
  poly fn;
  ideal LPi;
  for (j=1;j<=size(LP);j++)
  {
    if (j!=i)
    {
      F=idint(F,LP[j]);
    }
  }
  p=0; j=1;
  while ((p==0) and (j<=size(F)))
  {
    LPi=LP[i];
    attrib(LPi,"isSB",1);
    p=reduce(F[j],LPi);
    j++;
  }
  f=F[j-1];
  fn=nonzerodivisor(f,LP);
  setring(RR);
  def fr=imap(@P,f);
  def fnr=imap(@P,fn);
  return(list(fr,fnr));
}

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


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

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

// 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,PtoCrep(Prep(N1,M1)));
      LL[size(LL)+1]=T;
      if(equalideals(T[2][1],ideal(1))){te=1; break;}
    }
    if(te){break;}
  }
  ci=T[1];
  def Cs=submat(C1,1..2,ci);
  setring(RR);
  return(imap(@P,Cs));
}

// 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
// // NonNull: returns 1 if the poly f is nonnull on V(E) \ V(N), 0 otherwise.
// // Input:
// //   f: polynomial
// //   E: null ideal
// //   N: nonnull ideal
// // Output:
// //   1 if f is nonnul in V(P) \ V(Q),
// //   0 if it has zeroes inside.
// //  Works in @P
// proc NonNull(poly f, ideal E, ideal N)
// {
//   int te=1; int i;
//   def RR=basering;
//   setring(@P);
//   def fp=imap(RR,f);
//   def Ep=imap(RR,E);
//   def Np=imap(RR,N);
//   ideal H;
//   ideal Ef=Ep+fp;
//   for (i=1;i<=size(Np);i++)
//   {
//     te=radicalmember(Np[i],Ef);
//     if (te==0){break;}
//   }
//   setring(RR);
//   return(te);
// }

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

// 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];}
//     else {for (j=1;j<=size(Gr);j++){GrHi[size(GrHi)+1]=Gr[j]*H[i];}}
    if(comment>1){"Call to KSW with arguments "; GM; GrHi;  Nh;}
    KS=KSW0(GM,GrHi,Nh,comment);
    for (j=1;j<=size(KS);j++)
    {
      CGS[size(CGS)+1]=KS[j];
    }
    if(comment>1){"CGS after KSW = "; CGS;}
  }
  return(CGS);
}

// 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);}
  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];
        }
      }
    }
  }
  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))};
  //string("T_n=",n," new n0",size(L));
  return(L);
}

// Input: list(A)
//          A is a list of locally closed sets in Crep. A=[[P1,Q1],[P2,Q2],..,[Pr,Qr]]
//          whose union is a constructible set from
// Output list [Lev,C]:
//          where Lev is the Crep of the canonical first level of A, and
//          C is the complement of the first level Lev wrt to the closure of A. The elements are given in Crep,
static proc FirstLevel(list A)
{
  int n=size(A);
  list T=zeroone(n);
  ideal P; ideal Q;
  list Cb;  ideal Cc=ideal(1);
  int i; int j;
  intvec t;
  ideal Z=ideal(1);
  list C;
  for(i=1;i<=size(A);i++)
  {
    Z=intersect(Z,A[i][1]);
  }
  for(i=2; i<=size(T);i++)
  {
    t=T[i];
    P=ideal(1); Q=ideal(0);
    for(j=1;j<=n;j++)
    {
      if(t[n+1-j]==1)
      {
        Q=Q+A[j][2];
      }
      else
      {
        P=intersect(P,A[j][1]);
      }
    }
    //"T_Q="; Q; "T_P="; P;
    Cb=Crep(Q,P);
    //"T_Cb="; Cb;
    if(Cb[1][1]<>1)
    {
      C[size(C)+1]=Cb;
      Cc=intersect(Cc,Cb[1]);
    }
  }
  list Lev=list(Z,Cc);                //Crep(Z,Cc);
  if(size(C)>1){C=SimplifyUnion(C);}
  return(list(Lev,C));
}

// Input: list(A)
//          A is a list of locally closed sets in Crep. A=[[P1,Q1],[P2,Q2],..,[Pr,Qr]]
//          whose union is a constructible set from which the algorithm computes its canonical form.
// Output:
//     list [L,C] where
//          where L is the list of canonical levels of A,
//          and C is the list of canonical levels of the complement of A wrt to the closure of A.
//          All locally closed sets are given in Crep.
//          L=[[1,[p1,p2],[3,[p3,p4],..,[2r-1,[p_{2r-1},p_2r]]]] is the list of levels of A in Crep.
//          C=[[2,p2,p3],[4,[p4,p5],..,[2s,[p_{2s},p_{2s+1}]]]  is the list of levels of C,
//                                              the complement of S wrt the closure of A, in Crep.
proc ConsLevels(list A)
"USAGE:   ConsLevels(A);
          A is a list of locally closed sets in Crep. A=[[P1,Q1],[P2,Q2],..,[Pr,Qr]]
          whose union is a constructible set from which the algorithm computes its
          canonical form.
RETURN: A list with two elements: [L,C]
          where L is the list of canonical levels of A, and C is the list of canonical
          levels of the
          complement of A wrt to the closure of A.
          All locally closed sets are given in Crep.
          L=[[1,[p1,p2],[3,[p3,p4],..,[2r-1,[p_{2r-1},p_2r]]]]
          C=[[2,p2,p3],[4,[p4,p5],..,[2s,[p_{2s},p_{2s+1}]]]
NOTE: The basering R, must be of the form Q[a]
KEYWORDS: locally closed set, constructible set
EXAMPLE:  ConsLevels; shows an example"
{
  list L; list C; int i;
  list B; list T;
  for(i=1; i<=size(A);i++)
  {
    T=Crep(A[i][1],A[i][2]);
    B[size(B)+1]=T;
  }
  int level=0;
  list K;
  while(size(B)>0)
  {
    level++;
    K=FirstLevel(B);
    if(level mod 2==1){L[size(L)+1]=list(level,K[1]);}
    else{C[size(C)+1]=list(level,K[1]);}
    //"T_L="; L;
    //"T_C="; C;
    B=K[2];
    //"T_size(B)="; size(B);
  }
  return(list(L,C));
}
example
{ "EXAMPLE:"; echo = 2;
// Example of ConsLevels with nice geometrical interpretetion and 2 levels

if (defined(R)){kill R;}
if (defined(@P)){kill @P; kill @R; kill @RP;}

  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)=
  ConsLevels(L);

//----------------------------------------------------------------------
// Begin Problem S93
// Automatic theorem proving
// Generalized Steiner-Lehmus Theorem
// A.Montes, T.Recio

// Bisector of A(-1,0) = Bisector of B(1,0) varying C(a,b)

if(defined(R1)){kill R1;}
ring R1=(0,a,b),(x1,y1,x2,y2,p,r),dp;
ideal S93=(a+1)^2+b^2-(p+1)^2,
                    (a-1)^2+b^2-(1-r)^2,
                    a*y1-b*x1-y1+b,
                    a*y2-b*x2+y2-b,
                    -2*y1+b*x1-(a+p)*y1+b,
                    2*y2+b*x2-(a+r)*y2-b,
                    (x1+1)^2+y1^2-(x2-1)^2-y2^2;

short=0;
def GC93=grobcov(S93,"nonnull",ideal(b),"rep",1);
//"grobcov(S93,'nonnull',ideal(b),"rep",1)="; GC93;

list L0;
for(int i=1;i<=size(GC93);i++)
{
   L0[size(L0)+1]=GC93[i][3];
}

def GC93ab=grobcov(S93,"nonnull",ideal(a*b),"rep",1);

ring RR=0,(a,b),lp;

list L1;
L1=imap(R1,L0);
// L1=Total elements of the grobcov(S93,'nonnull',ideal(b),'rep',1) to be added=
L1;

// Adding all the elements of grobcov(S93,'nonnull',ideal(b),'rep',1)=
ConsLevels(L1);
}

//**************************** 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 reduced by the holes.
//     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;}
  //"T_B="; B; "E="; E; "N="; N;
  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)
    {
      //"T_B="; B;  // 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);
         if(Fenv==1)
        {
          L0[3]=determineF(GGG,F,E);
        }
      }
      else
      {
        L0[2]="Normal";
      }
    }
  }
  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);
   if(rh==0){env=1;} else{env=0;}
   setring RR;
          //L0[3]=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 ar pair of arguments:
//                "movdim", d  : by default movdim is 2 but it can be set to other values
//                "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 2. If one consider locus problems in higer dimensions, the number of
//                 mover coordinates (placed as the last variables) is to be given as an option.
// 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;
  //if(defined(@P)==1){te=1; kill @P; kill @R; kill @RP; }
  //setglobalrings();
  // Options
  list DD=#;
  int moverdim=nvars(R);
  int version=0;
  int nv=nvars(R);
  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]=="movdim"){moverdim=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=locusdim(G[i][2],moverdim);
    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);
     }
  }
  setring(R);
  ideal h;
  if(t1)
  {
    def 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];
        }
      }
    }
  }
  else{list P1;}
  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); // AQUI
      if(d==0)
      {
        L2[i][4]=string("Accumulation",L2[i][4]);
      }
      else{L2[i][4]=string("Degenerate",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];}
  }
  setring(R);
  def L=imap(@P,LN);
  for(i=1;i<=size(L);i++){if(size(L[i][2])==0){L[i][2]=ideal(1);}}
  //if(te==0){kill @R; kill @RP; kill @P;}
  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;
 }
 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 components: for each point in the component,
//               the number of solutions in the variables is finite, and
//               the solutions depend on the point in the component if the component is not 0-dimensional.
//               Special components:  for each point in the component,
//               the number of solutions in the variables is finite,
//               the component is not 0-dimensional, but the solutions do not depend on the
//               values of the parameters in the component.
//               Accumlation points: are 0-dimensional components for which it exist
//               an infinite number of solutions.
//               Degenerate components: are components of dimension greater than 0 for which
//               for every point in the component there exist infinite solutions.
//          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(G)
        The input must be the grobcov  of a parametrical ideal in Q[a][x],
        (a=parameters, x=variables). In fact a must be the tracer coordinates
        and x the mover coordinates and other auxiliary ones.
        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)
        who 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, Degenerate.
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,
         and represents the number of mever variables. they should be given as
         the last variables of the ring.
         \"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 system 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 2. If one consider locus problems in higer dimensions, the number of
         mover coordinates (placed as the last variables) is to be given as an option.
RETURN: The locus. The output is a list of the components ( C_1,.. C_n ) where
        C_i =  (p_i,(p_i1,..p_is_i), type_i,l evel_i )   and type_i can be
        'Normal', 'Special', Accumulation', 'Degenerate'. The 'Special' components
        return more information, namely the antiimage of the component, that
        is 0-dimensional for these kind of components.
        Normal components: for each point in the component,
        the number of solutions in the variables is finite, and
        the solutions depend on the point in the component.
        Special components:  for each point in the component,
        the number of solutions in the variables is finite. The
        antiimage of the component is 0-dimensional.
        Accumlation points: are 0-dimensional components whose
        antiimage is not zero-dimansional.
        Degenerate components: are components of dimension greater than 0
        whose antiimage is not-zero-diemansional.
        The components are given in canonical P-representation.
        The levels of a class of locus are 1,
        because they represent locally closed. sets.
NOTE: It can only be called after computing the grobcov of the
      parametrical ideal in generic representation ('ext',0),
      which is the default.
      The basering R, must be of the form Q[a_1,..,a_m][x_1,..,x_n].
KEYWORDS: geometrical locus, locus, loci.
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=#;
  int moverdim=nvars(R);
  int version=0;
  int nv=nvars(R);
  if(nv<4){version=1;}
  int comment=0;
  ideal Fm;
  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]=="system"){Fm=DD[2*i];}
    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
    if(DD[2*i-1]=="family"){poly F=DD[2*i];}
  }
  int j; int k;
  def B0=GG[1][2];
  def H0=GG[1][3][1][1];
  if (equalideals(B0,ideal(1)) or (equalideals(H0,ideal(0))==0))
  {
    return(locus0(GG,DD));
  }
  else
  {
    int n=nvars(R);
    ideal vmov=var(n-1),var(n);
    ideal N;
    intvec xw; intvec yw;
    for(i=1;i<=n-1;i++){xw[i]=0;}
    xw[n]=1;
    for(i=1;i<=n;i++){yw[i]=0;}
    yw[n-1]=1;
    poly px; poly py;
    int te=1;
    i=1;
    while( te and i<=size(B0))
    {
      if((deg(B0[i],xw)==1) and (deg(B0[i])==1)){px=B0[i]; te=0;}
      i++;
    }
    i=1; te=1;
    while( te and i<=size(B0))
    {
      if((deg(B0[i],yw)==1) and (deg(B0[i])==1)){py=B0[i]; te=0;}
      i++;
    }
    N=px,py;
    te=indepparameters(N);
    if(te)
    {
      string("locus detected that the mover must avoid point (",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{"Warning! Problem with more than one mover. Not able to solve it."; list L; return(L);}
  }
  def LL=locus0(nGP,DD);
  kill @RP; kill @P; kill @R;
  return(LL);
}
example
{
  "EXAMPLE:"; echo = 2;
  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;
  // System S96=
  S96;
  locus(grobcov(S96));
  kill R;
  // ********************************************
  ring R=(0,a,b),(x4,x3,x2,x1),dp;
  short=0;
  ideal S=(x1-3)^2+(x2-1)^2-9,
               (4-x2)*(x3-3)+(x1-3)*(x4-1),
               (3-x1)*(x3-x1)+(4-x2)*(x4-x2),
               (4-x4)*a+(x3-3)*b+3*x4-4*x3,
               (a-x1)^2+(b-x2)^2-(x1-x3)^2-(x2-x4)^2;
  // System S=
  S;
  locus(grobcov(S));
  kill R;
  //********************************************

  ring R=(0,x,y),(x1,x2),dp;
  short=0;
  ideal S=-(x - 5)*(x1 - 1) - (x2 - 2)*(y - 2),
                (x1 - 5)^2 + (x2 - 2)^2 - 4,
                -2*(x - 5)*(x2 - 2) + 2*(x1 - 5)*(y - 2);
  // System S=
  S;
  locus(grobcov(S));
}

//  locusdg(G):  Special routine for determining the locus of points
//                 of  geometrical constructions. Given a parametric ideal J with
//                 parameters (a_1,..a_m) and variables (x_1,..,xn),
//                 representing the system determining
//                 the locus of points (a_1,..,a_m) who verify certain
//                 properties, computing the grobcov G of
//                 J and applying to it locus, determines the different
//                 classes of locus components. The components can be
//                 Normal, Special, Accumulation point, Degenerate.
//                 The output are the components given in P-canonical form
//                 of at most 4 constructible sets: Normal, Special, Accumulation,
//                 Degenerate.
//                 The description of the algorithm and definitions is
//                 given in a forthcoming paper by Abanades, Botana, Montes, Recio.
//                 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 2. If onw consider locus problems in higer dimensions, the number of
//                 mover coordinates (placed as the last variables) is to be given as an option.
//
//  input:      The output of locus(G);
//  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 components: for each point in the component,
//               the number of solutions in the variables is finite, and
//               the solutions depend on the point in the component if the component is not 0-dimensional.
//               Special components:  for each point in the component,
//               the number of solutions in the variables is finite,
//               the component is not 0-dimensional, but the solutions do not depend on the
//               values of the parameters in the component.
//               Accumlation points: are 0-dimensional components for which it exist
//               an infinite number of solutions.
//               Degenerate components: are components of dimension greater than 0 for which
//               for every point in the component there exist infinite solutions.
//          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(L)   The call must be 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).

NOTE: It can only be called after computing the locus.
           Calling sequence:    locusdg(locus(grobcov(S)));
KEYWORDS: geometrical locus, locus, loci, 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;
  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;
  // System S96=
  S96;
  locus(grobcov(S96));
  locusdg(locus(grobcov(S96)));
  kill R;
  //********************************************
  ring R=(0,a,b),(x4,x3,x2,x1),dp;
  ideal S=(x1-3)^2+(x2-1)^2-9,
               (4-x2)*(x3-3)+(x1-3)*(x4-1),
               (3-x1)*(x3-x1)+(4-x2)*(x4-x2),
               (4-x4)*a+(x3-3)*b+3*x4-4*x3,
               (a-x1)^2+(b-x2)^2-(x1-x3)^2-(x2-x4)^2;
  short=0;
  locus(grobcov(S));
  locusdg(locus(grobcov(S)));
  kill R;
  //********************************************

  ring R=(0,x,y),(x1,x2),dp;
  short=0;
  ideal S=-(x - 5)*(x1 - 1) - (x2 - 2)*(y - 2),
                (x1 - 5)^2 + (x2 - 2)^2 - 4,
                -2*(x - 5)*(x2 - 2) + 2*(x1 - 5)*(y - 2);
  locus(grobcov(S));
  locusdg(locus(grobcov(S)));
}

// locusto: Transforms the output of locus, locusdg, envelop and envelopdg
//             into a string that can be reed from different computational systems.
// input:
//     list L: The output of locus
// output:
//     string s: The output of locus converted to a string readable by other programs
proc locusto(list L)
"USAGE: locusto(L);
          The argument must be the output of locus or locusdg or
          envelop or envelopdg.
          It transforms the output into a string in standard form
          readable in many languages (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) )
              -envelopdg(envelop(F,C))      -> locusto( envelopdg(envelop(F,C)) )
KEYWORDS: geometrical locus, locus, loci.
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;
  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))));
  kill R;
  //********************************************

  // 1. Take a fixed line l: x1-y1=0  and consider
  //    the family F of a lines parallel to l passing through the mover point M
  // 2. Consider a circle x1^2+x2^2-25, and a mover point M(x1,x2) on it.
  // 3. Compute the envelop of the family of lines.

  ring R=(0,x,y),(x1,y1),lp;
  poly F=(y-y1)-(x-x1);
  ideal C=x1^2+y1^2-25;
  short=0;

  // Curves Family F=
  F;
  // Conditions C=
  C;

  locusto(envelop(F,C));
  locusto(envelopdg(envelop(F,C)));
  kill R;
  //********************************************

  // 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.

  ring R=(0,x,y),(x1,y1,x2,y2),lp;
  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);

  short=0;

  // Curves Family F=
  F;
  // Conditions C=
  C;

  locusto(envelop(F,C));
  locusto(envelopdg(envelop(F,C)));
}

// Auxiliary routine
// locusdim
// input:
//    B:         ideal, a basis of a segment of the grobcov
//    dgdim: int, the dimension of the mover (for locus)
//                by default dgdim is equal to the number of variables
static proc locusdim(ideal B, list #)
{
  def R=basering;
  int dgdim;
  int nv=nvars(R);
  if (size(#)>0){dgdim=#[1];}
  else {dgdim=nv;}
  int d;
  list v;
  ideal vi;
  int i;
  for(i=1;i<=dgdim;i++)
  {
    v[size(v)+1]=varstr(nv-dgdim+i);
    vi[size(v)+1]=var(nv-dgdim+i);
  }
  ideal B0;
  for(i=1;i<=size(B);i++)
  {
    if(subset(variables(B[i]),vi))
    {
      B0[size(B0)+1]=B[i];
    }
  }
  def RR=ringlist(R);
  def RR0=RR;
  RR0[2]=v;
  def R0=ring(RR0);
  setring(R0);
  def B0r=imap(R,B0);
  B0r=std(B0r);
  d=dim(B0r);
  setring R;
  return(d);
}

static proc norspec(ideal F)
{
  def RR=basering;
  def Rx=ringlist(RR);

   def Rx=ringlist(RR);
  def @P=ring(Rx[1]);
  list Lx;
  Lx[1]=0;
  Lx[2]=Rx[2]+Rx[1][2];
  Lx[3]=Rx[1][3];
  Lx[4]=Rx[1][4];
  Rx[1]=0;
  def D=ring(Rx);
  def @RP=D+@P;
  exportto(Top,@R);      // global ring Q[a][x]
  exportto(Top,@P);      // global ring Q[a]
  exportto(Top,@RP);     // global ring K[x,a] with product order
  setring(RR);
}

// envelop
// Input: n arguments
//   poly F: the polynomial defining the family of curves in ring R=0,(x,y),(x1,..,xn),lp;
//   ideal C=g1,..,g_{n-1}:  the set of constraints
// Output: the components of the envolvent;
proc envelop(poly F, ideal C, list #)
"USAGE: envelop(F,C);
          The first argument F must be the family of curves of which
          on want to compute the envelop.
          The second argument C must be the ideal of conditions
          over the variables, and should contain as polynomials
          as the number of variables -1.
RETURN: The components of the envelop with its taxonomy:
           The taxonomy distinguishes 'Normal',
          'Special', 'Accumulation', 'Degenerate' components.
          In the case of 'Special' components, it also
          outputs the antiimage of the component
          and an integer (0-1). If the integer is 0
          the component is not a curve of the family and is
          not considered as 'Relevant' by the envelopdg routine
          applied to it, but is considered as 'Relevant' if the integer is 1.
NOTE: grobcov is called internally.
          The basering R, must be of the form Q[a][x] (a=parameters, x=variables).
KEYWORDS: geometrical locus, locus, loci, envelop
EXAMPLE:  envelop; shows an example"
{
  int tes=0; int i;   int j;
  def R=basering;
  if(defined(@P)==1){tes=1; kill @P; kill @R; kill @RP;}
  setglobalrings();
  // Options
  list DD=#;
  int moverdim=nvars(R);
  int version=0;
  int nv=nvars(R);
  if(nv<4){version=1;}
  int comment=0;
  ideal Fm;
  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]=="system"){Fm=DD[2*i];}
    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
  }
  int n=nvars(R);
  list v;
  for(i=1;i<=n;i++){v[size(v)+1]=var(i);}
  def MF=jacob(F);
  def TMF=transpose(MF);
  def Mg=MF;
  def TMg=TMF;
  for(i=1;i<=n-1;i++)
  {
    Mg=jacob(C[i]);
    TMg=transpose(Mg);
    TMF=concat(TMF,TMg);
  }
  poly J=det(TMF);
  ideal S=ideal(F)+C+ideal(J);
  DD[size(DD)+1]="family";
  DD[size(DD)+1]=F;
  def G=grobcov(S,DD);
   def L=locus(G, DD);
  return(L);
}
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.

  ring R=(0,x,y),(x1,y1,x2,y2),lp;
  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);

  short=0;

  // Curves Family F=
  F;
  // Conditions C=
  C;

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

// envelopdg
// Input: list L: the output of envelop(poly F, ideal C, list #)
// Output: the relevant components of the envolvent in dynamic geometry;
proc envelopdg(list L)
"USAGE: envelopdg(L);
          The input list L must be the output of the call to
          the routine 'envolop' of the family of curves
RETURN: The relevant components of the envelop in Dynamic Geometry.
           'Normal' and 'Accumulation' components are always considered
           'Relevant'. 'Special' components of the envelop outputs
           three objects in its characterization: 'Special', the antiimage ideal,
           and the integer 0 or 1, that indicates that the given component is
           formed (1) or is not formed (0) by curves of the family. Only if yes,
           'envelopdg' considers the component as 'Relevant' .
NOTE: It must be called to the output of the 'envelop' routine.
          The basering R, must be of the form Q[a,b,..][x,y,..].
KEYWORDS: geometrical locus, locus, loci, envelop.
EXAMPLE:  envelop; shows an example"
{
  list LL;
  list Li;
  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")){Li=L[i]; Li[3]="Relevant"; LL[size(LL)+1]=Li;}
    }
    else
    {
      if(typeof(L[i][3])=="list")
      {
        if(L[i][3][3]==1)
        {
          Li=L[i]; Li[3]="Relevant"; LL[size(LL)+1]=Li;
        }
      }
    }
  }
  return(LL);
}
example
{
  "EXAMPLE:"; echo = 2;

  // 1. Take a fixed line l: x1-y1=0  and consider
  //    the family F of a lines parallel to l passing through the mover point M
  // 2. Consider a circle x1^2+x2^2-25, and a mover point M(x1,x2) on it.
  // 3. Compute the envelop of the family of lines.

  ring R=(0,x,y),(x1,y1),lp;
  short=0;
  poly F=(y-y1)-(x-x1);
  ideal C=x1^2+y1^2-25;
  short=0;

  // Curves Family F=
  F;
  // Conditions C=
  C;

  envelop(F,C);
  envelopdg(envelop(F,C));
}

// 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 the 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)
{
//  int te0=0;
//  def RR=basering;
//  if(defined(@P)){te0=1;  def Rx=@R;  kill @P; setring RR;}
  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];
      }
      else{"ERROR p1 NOT FOUND";}
      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;
    }
  }
  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=PtoCrep(list(L[i]));
    Lc[size(Lc)+1]=Sc;
  }
  list S=ConsLevels(Lc)[1];
  list Sout;
  list Lev;
  for(i=1;i<=size(S);i++)
  {
    Lev=list(i,Prep(S[i][2][1],S[i][2][2]));
    Sout[size(Sout)+1]=Lev;
  }
  return(Sout);
}

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