source: git/Singular/LIB/grobcov.lib @ 20f2239

//////////////////////////////////////////////////////////////////////////////
version="version grobcov.lib 4.1.2.  January_2020 "; // $Id$;
           // version N9;  January 2020;
info="
LIBRARY:  grobcov.lib
          \"Groebner Cover for parametric ideals.\",
          Comprehensive Groebner Systems, Groebner Cover,
          Canonical Forms,  Parametric Polynomial Systems,
          Automatic Deduction of Geometric Theorems,
          Dynamic Geometry, Loci, Envelope, Constructible sets.
          See: A. Montes A, M. Wibmer,
          \"Groebner Bases for Polynomial Systems with parameters\",
          Journal of Symbolic Computation 45 (2010) 1391-1425.
          (https://www.mat.upc.edu//en/people/antonio.montes/).

IMPORTANT: Recently published book:

           A. Montes. \" The Groebner Cover\":
           Springer, Algorithms and Computation in Mathematics 27 (2019)
           ISSN 1431-1550
           ISBN 978-3-030-03903-5
           ISBN 978-3-030-03904-2  (e-Book)
           Springer Nature Switzerland AG 2018

           https://www.springer.com/gp/book/9783030039035

           The book can also be used as a user manual of all
           the  routines included in this library.
           It defines and proves all the theoretic results used
           in the library,  and shows examples of all the routines.
           There are many previous papers related to the subject,
           and the book actualices all the contents.

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

OVERVIEW: In 2010, the library was designed to contain
          Montes-Wibmer's algorithm for computing the
          Canonical Groebner Cover of a  parametric ideal.
          The central  routine is grobcov.
          Given  a  parametric ideal, grobcov outputs
          its Canonical  Groebner Cover, consisting of a set
          of triplets of (lpp, basis,  segment). The basis
          (after normalization) is the reduced  Groebner basis
          for each point of the segment. The segments
          are disjoint, locally closed and  correspond to
          constant lpp (leading power product) of the basis,
          and are represented in canonical representation.
          The segments cover the  whole parameter space.
          The output is canonical, it only depends on the
          given parametric ideal and the monomial order,
          because the segments  have
          different lpph of the homogenized system.
          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. Its existence was proved for the
          first time by Michael Wibmer \"Groebner bases for
          families of affine or projective schemes\",
          JSC, 42,803-834 (2007).

          grobcov uses a first algorithm cgsdr that outputs a
          disjoint  reduced Comprehensive Groebner System
          with constant lpp. For this purpose, in this library,
          the implemented algorithm is Kapur-Sun-Wang
          algorithm, because it is actually the most efficient
          algorithm known for this purpose.
          D. Kapur, Y. Sun, and D.K. Wang \"A New Algorithm
          for  Computing Comprehensive Groebner Systems\".
          Proceedings of ISSAC'2010, ACM Press, (2010), 29-36.

         The library has evolved to include new applications of
          the  Groebner Cover, and new theoretical developments
          have been done.

          A routine locus has been included to compute
          loci of points, and determining the taxonomy of the
          components.
          Additional routines to transform the output to string
          (locusdg,  locusto) are also included and used in the
          Dynamic Geometry software  GeoGebra. They were
          described in:
          M.A. Abanades, F. Botana, A. Montes, T. Recio:
          \''An Algebraic Taxonomy  for Locus Computation in
          Dynamic Geometry\''.
          Computer-Aided Design 56 (2014) 22-33.

          Routines for determining the generalized envelope of a family
          of hypersurfaces (envelop, AssocTanToEnv,
          FamElemsToEnvCompPoints) are also included.

          It also includes procedures for
          Automatic Deduction of Geometric Theorems (ADGT).

          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 determines the canonical
          levels of a constructible set. It is described in:
          J.M. Brunat, A. Montes, \"Computing the canonical
          representation of constructible sets\".
          Math.  Comput. Sci. (2016) 19: 165-178.
          A complementary routine Levels transforms the output
          of ConsLevels into the proper locally closed sets
          forming the levels of the constructible.
          Antoher complementary routine Grob1Levels
          has been included to select the locally closed sets of
          the segments of the grobcov that correspond to basis
          different from 1, add them together and return
          the canonical form of this constructible set.
          More recently (2019) given two locally closed sets
          in canonical form the new routine DifConsLCSets
          determines a set of locally closed sets equivalent to
          the difference them. The description of the
          routine is submitted to the Journal of Symbolic Computation.
          This routine can be also used internally by ADGT
          with the option \"neg\",1 . With this option
          DifConsLCSets is used for the negative
          hypothesis and thesis in ADGT.

          This version was finished on 1/2/2020,


NOTATIONS: Before calling any routine of the library grobcov,
          the user  must define the ideal Q[a][x], and all the
          input polynomials and  ideals defined on it.
          Internally the routines define and use also other
          ideals: Q[a], Q[x,a] and so on.

PROCEDURES:

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

cgsdr(F); Is the procedure for obtaining a first disjoint,
          reduced  Comprehensive Groebner  System that is
          used in grobcov, but can also be used independently
          if only a CGS is required. It is a more efficient routine
          than buildtree (the own routine of 2010 that  is no
          more available).
          Now, Kapur-Sun-Wang (KSW) algorithm is used.

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

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

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

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

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

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

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

locus(G); Special routine for determining the geometrical
          locus of points verifying given conditions. Given
          a parametric ideal J in Q[x1,..,xn][u1,..,um]
          with parameters (x1,..,xn) and variables
          (u1,..,um), representing the system determining
          the n-dimensional locus with tracer point (x1,..,xn)
          verifying certain properties,
          one can apply locus to the system F, for obtaining the locus.

          locus provides all the components of the locus and
          determines their taxonomy, that can be:
           \"Normal\", \"Special\", \"Accumulation\",
          \"Degenerate\".
          The mover  are the u-variables.
          The user can ventually restrict them to a subset of them
          for geometrical reasons but this can change the true
          taxonomy.

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 \"Normal\", and\"Accumulation\"
          components.

envelop(F,C); Special routine for determining the envelop
          of a family of hyper-surfaces F  in
          Q[x1,..,xn][t1,..,tm] depending on an ideal of
          constraints C in Q[t1,..,tm]. It computes the
          locus of the envelop, and determines the
          different components as well as its taxonomy:
          \"Normal\", \"Special\", \"Accumulation\",
          \"Degenerate\". (See help for locus).

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

stdlocus(F); Simple procedure to determine the components
          of the locus, alternative to  locus that uses only
          standard GB computation. Cannot determine the
          taxonomy of  the irreducible components.

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

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

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

WLemma(F,A); Given an ideal F in Q[a][x] and an ideal A
          in Q[a], it returns the list (lpp,B,S)  were B is the
          reduced Groebner basis of the specialized F over
          the segment S, subset of V(A) with top A,
          determined by Wibmer's Lemma.
          S is determined in P-representation
          (or optionally in C-representation). The basis is
          given by I-regular functions.

WLCGS(F);  Given a parametric ideal F in Q[a][x]
          determines a CGS in full-representation using WLemma

intersectpar(L); Auxiliary routine. Given a list of ideals defined on K[a][x]
         it determines the intersection of all of them in K[x,a]

ADGT(H,T,H1,T1); Given 4 ideals H,T,H1,T1 in Q[a][x], corresponding to
        a problem of Automatic Deduction of Geometric Theorems,
        it determines the supplementary conditions over the parameters
        for the Proposition (H and not H1) => (T and not T1) to be a
        Theorem.
        If H1=1 then H1 is not considered, and analogously for T1.

ConsLevels(A); Given a list of locally closed sets, constructs the
        canonical representation of the levels of A an its complement.

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

Grob1Levels(G); From the output of grobcov, Grob1Levels selects the segments
        of G with basis different from 1 (having solutions), and determines
        the levels of the constructible set formed by them.

DifConsLCSets(A,B); given the canonical forms of the constructible sets A and B,
         A=[a1,a2,..,ak],  B=[b1,b2,...,bj], DifConsLCSets returns a list of
         locally closed sets of the set A minus B, that can be transformed into the
         canonical form of A minus B applying ConsLevels.

SEE ALSO: compregb_lib
";

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

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

// Development of the library:
// Library grobcov.lib
// (Groebner Cover):
// Release 0: (public)
// Initial data: 21-1-2008
// Uses buildtree for cgsdr
// Final data: 3-7-2008
// Release 2: (prived)
// Initial data: 6-9-2009
// Last version using buildtree for cgsdr
// Final data: 25-10-2011
// Release B: (prived)
// Initial data: 1-7-2012
// Uses both buildtree and KSW for cgsdr
// Final data: 4-9-2012
// Release D. Includes routine locus
// Release G: (public)
// Initial data: 4-9-2012
// Uses KSW algorithm for cgsdr
// Final data: 21-11-2013
// Release K: Includes routine envelop
// Release L: (public)
// New routine ConsLevels: 25-1-2016
// New routine Levels: 25-1-2016
// New routine Grob1Levels: 25-1-2016
// Updated locus: 10-7-2015 (uses ConsLevels)
// Release M: (public)
// New routines for computing the envelope of a family of
//    hyper-surfaces and associated questions: 22-4-2016: 20-9-2016
// New routine WLemma for computing the result of
//    Wibmer's Lemma:  19-9-2016
// Final data October 2016
// Updated locus (messages) // Updated locus (messages)
// Final data Mars 2017
// Release N4: (public)
// New routine ADGT for Automatic Discovery of Geometric Theorems: 21-1-2018
// Final data February 2018
// Release N8: July 2019. Actualized versions of the routines and options
// Release N9: December 2019. New routine DifConsLCSets,
// Updated also ADGT to use as option DifConsLCSets

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

// elimfromlistel: elimine the ideal J from the list L
// Input: list L;  list of ideals
//           ideal J;  a possible element of L
// Output:ideal K with the elements of L different from J
//        ring Q[x1,..,xm]
static proc elimidealfromlist(list L,ideal J)
{
  int i;
  int j=0;
  list K;
  if (size(L)==0){return(L);}
  for (i=1;i<=size(L);i++)
  {
    if (not(equalideals(J,L[i])))
    {
      j++;
      K[j]=L[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==1))
  {
      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
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,5);
    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);
}

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;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  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);
}

// input ideal J, ideal K
// output 1 if all the polynomials in J are members of K
//            0 if not
 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);
//}

// 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;
  def Rx=ringlist(RR);
  if(size(Rx[1])==4)
  {
    te=1;
    def P=ring(Rx[1]);
    Rx[1]=0;
    def D=ring(Rx);
    def RP=D+P;
    setring(RP);
    def Ja=imap(RR,J);
  }
  else{ def Ja=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]);}}
  if(te==1)
  {
    setring(RR);
    def Jc=imap(RP,Jb);
  }
  else{def Jc=Jb;}
  // if(te){kill @R; kill @RP; kill @P;}
  return(Jc);
};

// simpqcoeffs : simplifies a quotient of two polynomials
// input: two coefficients (or terms), that are considered as a quotient
// output: the two coefficients 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)
//   mu*f=sum(q_i*F_i)+r
//   no monomial of r is divisible by no lpp of F
proc pdivi(poly f,ideal F)
"USAGE: pdivi(poly f,ideal F);
          poly f: the polynomial in Q[a][x] to be divided
          ideal F: the divisor ideal in Q[a][x].
          (a=parameters, x=variables).
RETURN: A list (poly r, ideal q, poly m). r is the remainder
          of the pseudodivision, q is the set of quotients,
          and m is the coefficient factor by which f is to
          be multiplied.
NOTE: pseudodivision of a poly f by an ideal F in Q[a][x].
          Returns a list (r,q,m) such that
          m*f=r+sum(q.F),
          and no lpp of a divisor divides a pp of r.
KEYWORDS: division; reduce
EXAMPLE:  pdivi; shows an example"
{
  F=simplify(F,2);
  int i;
  int j;
  poly v=1;
  for(i=1;i<=nvars(basering);i++){v=v*var(i);}
  poly r=0;
  poly mu=1;
  def p=f;
  ideal q;
  for (i=1; i<=ncols(F); i++){q[i]=0;};
  ideal lpf;
  ideal lcf;
  for (i=1;i<=ncols(F);i++){lpf[i]=leadmonom(F[i]);}
  for (i=1;i<=ncols(F);i++){lcf[i]=leadcoef(F[i]);}
  poly lpp;
  poly lcp;
  poly qlm;
  poly nu;
  poly rho;
  int divoc=0;
  ideal qlc;
  while (p!=0)
  {
    i=1;
    divoc=0;
    lpp=leadmonom(p);
    lcp=leadcoef(p);
    while (divoc==0 and i<=size(F))
    {
      qlm=lpp/lpf[i];
      if (qlm!=0)
      {
        qlc=simpqcoeffs(lcp,lcf[i]);
        //string("T_i=",i,";  qlc=",qlc);
        nu=qlc[2];
        mu=mu*nu;
        rho=qlc[1]*qlm;
        //"T_nu="; nu; "mu="; mu; "rho="; rho;
        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;
        //"T_q[i]="; q[i];
        divoc=1;
      }
      else {i++;}
    }
    if (divoc==0)
    {
      r=r+lcp*lpp;
      p=p-lcp*lpp;
    }
    //"T_r="; r; "p="; p;
  }
  list res=r,q,mu;
  return(res);
}
example
{
  "EXAMPLE:"; echo = 2;
  if(defined(R)){kill R;}
  ring R=(0,a,b,c),(x,y),dp;
  short=0;
  // Divisor
  poly f=(ab-ac)*xy+(ab)*x+(5c);
  // Dividends
  ideal F=ax+b,cy+a;
  // (Remainder, quotients, factor)
  def r=pdivi(f,F);
  r;
  // Verifying the division
  r[3]*f-(r[2][1]*F[1]+r[2][2]*F[2]+r[1]);
}

//*****************************

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

// facvar: Returns all the free-square factors of the elements
//         of ideal J (non repeated). Integer factors are ignored,
//         even 0 is ignored. It can be called from ideal Q[a][x], but
//         the given ideal J must only contain poynomials in the
//         parameters a.
//         Operates in the ring Q[a], but can be called from ring Q[a][x],
// 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.
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;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  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 Q[a][x]
// 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;
  if((f==0) or (equalideals(N,ideal(1)))){ return(f);}
  def v=variables(f);
  int i;
  poly X=1;
  for(i=1;i<=size(v);i++){X=X*v[i];}
  matrix M=coef(f,X);
  list Mc;
  for(i=1;i<=ncols(M);i++){Mc[i]=M[2,i];}
  // "T_M="; M;
  // "T_Mc=";Mc;

  poly g=M[2,1];
  if (size(M)!=2)
  {
    for(i=2;i<=size(M) div 2;i++)
    {
      g=gcd(g,M[2,i]);
    }
  }
  // "T_g="; g;
  if (g==1){ return(f);}
  else
  {
    def wg=factorize(g);
    // "T_wg="; wg;
    if (wg[1][1]==1){ return(f);}
    else
    {
      poly simp=1;
      int te;
      for(i=1;i<=size(wg[1]);i++)
      {
        te=inconsistent(E+wg[1][i],N);
        if(te)
        {
          simp=simp*(wg[1][i])^(wg[2][i]);
        }
      }
    }
    if (simp==1){ return(f);}
    else
    {
      //def simp0=imap(P,simp);
      def f2=f/simp;
      return(f2);
    }
  }
}

//*******************

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

//*******************

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

// sortideal: sorts the polynomials in an ideal by lm in ascending order
static proc sortideal(ideal Fi)
{
  def RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  Rx[1]=0;
  def D=ring(Rx);
  def RP=D+P;
  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;
  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;
  ideal G;
  ideal H;
  int i;
  if (size(F)==0){return(ideal(0));}
  for (i=1;i<=size(F);i++)
  {
    H=delfromideal(F,i);
    G[i]=pdivi(F[i],H)[1];
  }
  return(G);
}


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

// Crep0
// Computes the C-representation of V(N) \ V(M).
// input:
//    ideal N (null ideal) (not necessarily radical nor maximal)
//    ideal M (hole ideal) (not necessarily containing N)
// output:
//    the list (a,b) of the canonical ideals
//    the Crep0 of V(N) \ V(M)
// Assumed to be called in a ring Q[x] (example @P)
static proc Crep0(ideal N, ideal M)
{
  list l;
  ideal Np=std(N);
  if (equalideals(Np,ideal(1)))
  {
    l=ideal(1),ideal(1);
    return(l);
  }
  int i;
  list L;
  ideal Q=Np+M;
  ideal P=ideal(1);
  L=minGTZ(Np);
  for(i=1;i<=size(L);i++)
  {
    L[i]=std(L[i]);
    if(idcontains(L[i],Q)==0)
    {
      P=intersect(P,L[i]);
    }
  }
  P=std(P);
  Q=std(radical(Q+P));
  if(equalideals(P,Q)){return(l);}
  list T=P,Q;
  return(T);
}

// Prep
// Computes the P-representation of V(N) \ V(M).
// input:
//    ideal N (null ideal) (not necessarily radical nor maximal)
//    ideal M (hole ideal) (not necessarily containing N)
// output:
//    the ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r)));
//    the Prep of V(N) \ V(M)
// Assumed to be called in the ring ring Q[a][x]. But the data must only contain parameters.
proc Prep(ideal N, ideal M)
"USAGE: Prep(ideal N,ideal M);
       ideal N (null ideal) (not necessarily radical nor
       maximal) in Q[a]. (a=parameters, x=variables).
       ideal M (hole ideal) (not necessarily containing N)
       in Q[a]. To be called in a ring Q[a][x] or a ring
       Q[a]. But the ideals can contain only the
       parameters in Q[a].
RETURN: The canonical P-representation of the locally closed
       set V(N) \ V(M)
       Output: [Comp_1, .. , Comp_s ] where
       Comp_i=[p_i,[p_i1,..,p_is_i]]
KEYWORDS: locally closed set; canonical form
EXAMPLE:  Prep; shows an example"
{
  int te;
  def RR=basering;
  def Rx=ringlist(RR);
  if(size(Rx[1])==4)
  {
    def P=ring(Rx[1]);
    te=1; setring(P); ideal Np=imap(RR,N); ideal Mp=imap(RR,M);
  }
  else {te=0; def Np=N; def Mp=M;}
  def La=Prep0(Np,Mp);
  if(te==1) {setring(RR); def LL=imap(P,La); }
  if(te==0){def LL=La;}
  return(LL);
}
example
{
  "EXAMPLE:"; echo = 2;
  short=0;
  if(defined(R)){kill R;}
  ring R=0,(a,b,c),lp;
  ideal p=a*b;;
  ideal q=a,b-1;
  // P-representation of V(p) \ V(q)
  Prep(p,q);
}

// Prep0
// Computes the P-representation of V(N) \ V(M).
// input:
//    ideal N (null ideal) (not necessarily radical nor maximal)
//    ideal M (hole ideal) (not necessarily containing N)
// output:
//    the ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r)));
//    the Prep of V(N) \ V(M)
// Assumed to be called in a ring Q[x] (example @P)
static proc Prep0(ideal N, ideal M)
{
  int te;
  if (N[1]==1)
  {
    return(list(list(ideal(1),list(ideal(1)))));
  }
  int i; int j; list L0;
  list Ni=minGTZ(N);
  list prep;
  for(j=1;j<=size(Ni);j++)
  {
    option(redSB);
    Ni[j]=std(Ni[j]);
  }
  list Mij;
  for (i=1;i<=size(Ni);i++)
  {
    Mij=minGTZ(Ni[i]+M);
    for(j=1;j<=size(Mij);j++)
    {
      option(redSB);
      Mij[j]=std(Mij[j]);
    }
    if ((size(Mij)==1) and (equalideals(Ni[i],Mij[1])==1)){;}
    else
    {
        prep[size(prep)+1]=list(Ni[i],Mij);
    }
  }
  //"T_before="; prep;
  if (size(prep)==0){prep=list(list(ideal(1),list(ideal(1))));}
  //"T_Prep="; prep;
  //def Lout=CompleteA(prep,prep);
  //"T_Lout="; Lout;
  return(prep);
}

// PtoCrep
// Computes the C-representation from the P-representation.
// input:
//    list ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r)));
//         the P-representation of V(N) \ V(M)
// output:
//    list (ideal ida, ideal idb)
//    the C-representaion of V(N) \ V(M) = V(ida) \ V(idb)
// Assumed to be called in the ring Q[a] or Q[x]
proc PtoCrep(list L)
"USAGE: PtoCrep(list L)
       list L=  [ Comp_1, .. , Comp_s ] where
       list Comp_i=[p_i,[p_i1,..,p_is_i] ], is the
       P-representation of a locally closed set
       V(N) \ V(M). To be called in a ring Q[a][x]
       or a ring Q[a]. But the ideals can contain
       only the parameters in Q[a].
RETURN:The canonical C-representation [P,Q] of the
       locally closed set. A pair of radical ideals with
       P included in Q, representing the
       set V(P) \ V(Q)
KEYWORDS: locally closed set; canonical form
EXAMPLE:  PtoCrep; shows an example"
{
  int te;
  def RR=basering;
  def Rx=ringlist(RR);
  if(size(Rx[1])==0){return(PtoCrep0(L));}
  else
  {
    def P=ring(Rx[1]);
    setring(P);
    list Lp=imap(RR,L);
    def LLp=PtoCrep0(Lp);
    setring(RR);
    def LL=imap(P,LLp);
    return(LL);
  }
}
example
{
 "EXAMPLE:"; echo = 2;
  if(defined(R)){kill R;}
  ring R=0,(a,b,c),lp;
  short=0;
  ideal p=a*(a^2+b^2+c^2-25);
  ideal q=a*(a-3),b-4;
  // C-representaion of V(p) \ V(q)
  def Cr=Crep(p,q);
  Cr;
  // C-representation of V(p) \ V(q)
  def L=Prep(p,q);
  L;
  PtoCrep(L);
}

// PtoCrep0
// Computes the C-representation from the P-representation.
// input:
//    list ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r)));
//         the P-representation of V(N) \ V(M)
// output:
//    list (ideal ida, ideal idb)
//    the C-representation of V(N) \ V(M) = V(ida) \ V(idb)
// Assumed to be called in a ring Q[x] (example @P)
static proc PtoCrep0(list L)
{
  int te=0;
  def Lp=L;
  int i; int j;
  ideal ida=ideal(1); ideal idb=ideal(1); list Lb; ideal N;
  for (i=1;i<=size(Lp);i++)
  {
    option(returnSB);
    //"T_Lp[i]="; Lp[i];
    N=Lp[i][1];
    Lb=Lp[i][2];
    //"T_Lb="; Lb;
    ida=intersect(ida,N);
    for(j=1;j<=size(Lb);j++)
    {
      idb=intersect(idb,Lb[j]);
    }
  }
  //idb=radical(idb);
  def La=list(ida,idb);
  return(La);
}

// input: F a parametric ideal in Q[a][x]
// output: a disjoint and reduced Groebner System.
//      It uses Kapur-Sun-Wang algorithm, and with the options
//      can compute the homogenization before  (('can',0) or ( 'can',1))
//      and dehomogenize the result.
proc cgsdr(ideal F, list #)
"USAGE: cgsdr(ideal F);
       F: ideal in Q[a][x] (a=parameters, x=variables) to be
       discussed. Computes a disjoint, reduced Comprehensive
       Groebner System (CGS). cgsdr is the starting point of
       the fundamental routine grobcov.
       The basering R, must be of the form Q[a][x],
       (a=parameters, x=variables),
       and should be defined previously.
RETURN: Returns a list T describing a reduced and disjoint
       Comprehensive Groebner System (CGS). The output
       is a list of  (full,hole,basis), where the  ideals
       full and hole represent the segment V(full) \ V(hole).
       With option (\"out\",0) the segments are grouped
       by leading power products (lpp) of the reduced
       Groebner basis and given in P-representation.
       The returned list is of the form:
        [ [lpp, [num,basis,segment],...,
                [num,basis,segment],lpph],  ... ,
          [lpp, [num,basis,segment],...,
                [num,basis,segment],lpph] ].
       The bases are the reduced Groebner bases (after
       normalization) for each point of the corresponding
       segment. The third element lpph of each lpp
       segment is the lpp of the homogenized ideal
       used ideal in the CGS  as a string, that
       is shown only when option  (\"can\",1) is used.
       With option (\"can\",0) the homogenized basis is used.
       With option (\"can\",1) the homogenized ideal is used.
       With option (\"can\",2) the given basis is used.
       With option (\"out\",1) (default) only KSW is applied and
       segments are given as difference of varieties and are
       not grouped The returned list is of the form:
       [[E,N,B],..[E,N,B]]
       E is the top variety
       N is the hole variety.
       Segment = V(E) \ V(N)
       B is the reduced Groebner basis
OPTIONS:  An option is a pair of arguments: string, integer.
       To modify the default options, pairs of arguments
       -option name, value- of valid options must be
       added to the call. Inside grobcov the default option
       is \"can\",1. It can be used also with option
       \"can\",0 but then the output is not the canonical
       Groebner Cover. grobcov cannot be used with
       option \"can\",2.
       When cgsdr is called directly, the options are
       \"can\",0-1-2: The default value is \"can\",2.
       In this case no homogenization is done. With option
       (\"can\",0) the given basis is homogenized,
       and with option (\"can\",1) the whole given ideal
       is homogenized before computing the cgs
       and dehomogenized after.
       With option (\"can\",0) the homogenized basis is used.
       With option (\"can\",1) the homogenized ideal is used.
       With option (\"can\",2) the given basis is used.
       \"null\",ideal E: The default is (\"null\",ideal(0)).
       \"nonnull\",ideal N: The default (\"nonnull\",ideal(1)).
       When options \"null\" and/or \"nonnull\" are given,
       then the parameter space is restricted to V(E) \ V(N).
       \"comment\",0-1: The default is (\"comment\",0).
       Setting (\"comment\",1) will provide information
       about the development of the computation.
       \"out\",0-1: (default is 1) the output segments are
       given as as difference of varieties.
       With option \"out\",0 the output segments are
       given in P-representation and  the segments
       grouped by lpp.
       With options (\"can\",0) and (\"can\",1) the option
       (\"out\",1) is set to (\"out\",0) because
       it is not compatible.
       One can give none or whatever of these options.
       With the default options (\"can\",2,\"out\",1),
       only the Kapur-Sun-Wang algorithm is computed.
       The algorithm used is:
          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.
       It is very efficient but is only the starting point
       for the computation of grobcov.
       When grobcov is computed, the call to cgsdr
       inside uses specific options that are
       more expensive (\"can\",0-1,\"out\",0).
KEYWORDS: CGS; disjoint; reduced; Comprehensive Groebner System
EXAMPLE:  cgsdr; shows an example"
{
  int te;
  def RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  // if(size(Rx[1])==4){te=1; Rx[1]=0; def D=ring(Rx); def RP=D+P;}
  // else{te=0;} //setglobalrings();}
  // INITIALIZING OPTIONS
  int i; int j;
  def E=ideal(0);
  def N=ideal(1);
  int comment=0;
  int can=2;
  int out=1;
  poly f;
  ideal B;
  int start=timer;
  list L=#;
  for(i=1;i<=size(L) div 2;i++)
  {
    if(L[2*i-1]=="null"){E=L[2*i];}
    else
    {
      if(L[2*i-1]=="nonnull"){N=L[2*i];}
      else
      {
        if(L[2*i-1]=="comment"){comment=L[2*i];}
        else
        {
          if(L[2*i-1]=="can"){can=L[2*i];}
          else
          {
            if(L[2*i-1]=="out"){out=L[2*i];}
          }
        }
      }
    }
  }
  //if(can==2){out=1;}
  B=F;
  if ((printlevel) and (comment==0)){comment=printlevel;}
  if((can<2) and (out>0)){"Option out,1 is not compatible with can,0,1"; out=0;}
  // DEFINING OPTIONS
  list LL;
  LL[1]="can";     LL[2]=can;
  LL[3]="comment"; LL[4]=comment;
  LL[5]="out";     LL[6]=out;
  LL[7]="null";    LL[8]=E;
  LL[9]="nonnull"; LL[10]=N;
  if(comment>=1)
  {
    " "; string("Begin cgsdr with options: ",LL);
  }
  int ish;
  for (i=1;i<=size(B);i++){ish=ishomog(B[i]); if(ish==0){break;};}
  if (ish)
  {
    if(comment>0){" "; string("The given system is homogneous");}
    def GS=KSW(B,LL);
    //can=0;
  }
  else
  {
  // ACTING DEPENDING ON OPTIONS
  if(can==2)
  {
    // WITHOUT HOMOHGENIZING
    if(comment>0){" "; string("Option of cgsdr: do not homogenize");}
    def GS=KSW(B,LL);
    // setglobalrings();
  }
  else
  {
    if(can==1)
    {
      // COMPUTING THE HOMOGOENIZED IDEAL
      if(comment>0){" "; string("Homogenizing the whole ideal: option can=1");}
      list RRL=ringlist(RR);
      RRL[3][1][1]="dp";
      def Pa=ring(RRL[1]);
      list Lx;
      Lx[1]=0;
      Lx[2]=RRL[2]+RRL[1][2];
      Lx[3]=RRL[1][3];
      Lx[4]=RRL[1][4];
      RRL[1]=0;
      def D=ring(RRL);
      def RP=D+Pa;
      setring(RP);
      def B1=imap(RR,B);
      option(redSB);
      if(comment>0){" ";string("Basis before computing its std basis="); B1;}
      B1=std(B1);
      if(comment>0){" ";string("Basis after computing its std basis="); B1;}
      setring(RR);
      def B2=imap(RP,B1);
    }
    else
    { // (can=0)
       if(comment>0){" "; string( "Homogenizing the basis: option can=0");}
      def B2=B;
    }
    // COMPUTING HOMOGENIZED CGS
    poly @t;
    ring H=0,@t,dp;
    def RH=RR+H;
    setring(RH);
    // setglobalrings();
    def BH=imap(RR,B2);
    def LH=imap(RR,LL);
    //"T_BH="; BH;
    //"T_LH="; LH;
    for (i=1;i<=size(BH);i++)
    {
      BH[i]=homog(BH[i],@t);
    }
    if (comment>0){" "; string("Homogenized system = "); BH;}
    def RHx=ringlist(RH);
    def PH=ring(RHx[1]);
    RHx[1]=0;
    def DH=ring(RHx);
    def RPH=DH+PH;
    def GSH=KSW(BH,LH);
    //"T_GSH="; GSH;
    //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]));}
        }
      }
    }
  }
  return(GS);
}
example
{
  "EXAMPLE:"; echo = 2;
  // Casas conjecture for degree 4:
  if(defined(R)){kill R;}
  ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp;
  short=0;
  ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0),
          x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1),
          x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0),
          x2^2+(2*a3)*x2+(a2),
          x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0),
          x3+(a3);
  cgsdr(F);
}

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

// LCUnion
// Given a list of the P-representations of locally closed segments
// for which we know that the union is also locally closed
// it returns the P-representation of its union
// input:  L list of segments in P-representation
//      ((p_j^i,(p_j1^i,...,p_jk_j^i | j=1..t_i)) | i=1..s )
//      where i represents a segment
// output: P-representation of the union
//       ((P_j,(P_j1,...,P_jk_j | j=1..t)))
static proc LCUnion(list LL)
{
  def RR=basering;
  def Rx=ringlist(RR);
  def PP=ring(Rx[1]);
  setring(PP);
  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]];
  }
  //"T_P="; P;
  // 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(PP,T);
  return(TT);
}

// LCUnionN
// 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 LCUnionN(list L)
{
  int i; int j; int k; list H; list C; list T;
  list L0; list P0; list P; list Q0; list Q;
  //"T_L="; L;
  for (i=1;i<=size(L);i++)
  {
    P0[size(P0)+1]=L[i][1];
    for (j=1;j<=size(L[i]);j++)
    {
      L0[size(L0)+1]=intvec(i,j);
    }
  }
  //"T_P0="; P0;
  Q0=selectminideals(P0);
  //"T_Q0="; Q0;
  for (i=1;i<=size(Q0);i++)
  {
    //Q[i]=L0[Q0[i]];
    P[i]=L[Q0[i][1]];// [Q[i][2]];
  }
  //"T_P="; P;
  // 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))
  // list C;
  for (k=1;k<=size(Q0);k++)
  {
    kill C; list C;
    H=P[k][2]; // holes of P[k][1]
    for (i=1;i<=size(L);i++)
    {
      if (i!=Q0[k]) // (i!=Q0[k])
      {
        //for (j=1;j<=size(L[i]);j++)
        //{
          C[size(C)+1]=L[i];
        //}
      }
    }
    T[size(T)+1]=list(P[k][1],addpart(H,C)); // Q0[k],
  }
  return(T);
}


// 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 Q[a] 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;
  //          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]))
        {
          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];
            }
          }
          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 Q[a] ring.
static proc addpartfine(list H, list C0)
{
  //"T_H="; H;
  int i; int j; int k; int te; intvec notQ; int l; list sel;
  intvec jtesC;
  if ((size(H)==1) and (equalideals(H[1],ideal(1)))){return(H);}
  if (size(C0)==0){return(H);}
  list newQ; list nQ; list Q; list nQ1; list Q0;
  def Q1=H;
  //Q1=sortlistideals(Q1,idbefid);
  def C=C0;
  while(equallistideals(Q0,Q1)==0)
  {
    Q0=Q1;
    i=0;
    Q=Q1;
    kill notQ; intvec notQ;
    while(i<size(Q))
    {
      i++;
      for(j=1;j<=size(C);j++)
      {
        te=idcontains(Q[i],C[j][1]);
        if(te)
        {
          for(k=1;k<=size(C[j][2]);k++)
          {
            if(idcontains(Q[i],C[j][2][k]))
            {
              te=0; break;
            }
          }
          if (te)
          {
            if ((size(notQ)==1) and (notQ[1]==0)){notQ[1]=i;}
            else{notQ[size(notQ)+1]=i;}
            kill newQ; list newQ;
            for(k=1;k<=size(C[j][2]);k++)
            {
              nQ=minGTZ(Q[i]+C[j][2][k]);
              for(l=1;l<=size(nQ);l++)
              {
                option(redSB);
                nQ[l]=std(nQ[l]);
                newQ[size(newQ)+1]=nQ[l];
              }
            }
            sel=selectminideals(newQ);
            kill nQ1; list nQ1;
            for(l=1;l<=size(sel);l++)
            {
              nQ1[l]=newQ[sel[l]];
            }
            newQ=nQ1;
            for(l=1;l<=size(newQ);l++)
            {
              Q[size(Q)+1]=newQ[l];
            }
            break;
          }
        }
      }
    }
    kill Q1; list Q1;
    for(i=1;i<=size(Q);i++)
    {
      if(memberpos(i,notQ)[1]==0)
      {
        Q1[size(Q1)+1]=Q[i];
      }
    }
    sel=selectminideals(Q1);
    kill nQ1; list nQ1;
    for(l=1;l<=size(sel);l++)
    {
      nQ1[l]=Q1[sel[l]];
    }
    Q1=nQ1;
  }
  if(size(Q1)==0){Q1=ideal(1),ideal(1);}
  return(Q1);
}

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

// Auxiliary routine
// precombine
// input:  L: list of ideals (works in @P)
// output: F0: ideal of polys. F0[i] is a poly in the intersection of
//             all ideals in L except in the ith one, where it is not.
//             L=(p1,..,ps);  F0=(f1,..,fs);
//             F0[i] \in intersect_{j#i} p_i
static proc precombine(list L)
{
  int i; int j; int tes;
  def RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  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],5)!=0){tes=0; F[i]=L3[i][j];}
    }
    if (tes){"ERROR a polynomial in all p_j except p_i was not found";}
  }
  setring(RR);
  def F0=imap(P,F);
  return(F0);
}

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

// Central routine for grobcov: ideal F is assumed to be homogeneous
// gcover
// input: ideal F: a generating set of a homogeneous ideal in Q[a][x]
//    list #: optional
// output: the list
//   S=((lpp, generic basis, Prep, Crep),..,(lpp, generic basis, Prep, Crep))
//      where a Prep is ( (p1,(p11,..,p1k_1)),..,(pj,(pj1,..,p1k_j)) )
//            a Crep is ( ida, idb )
static proc gcover(ideal F,list #)
{
  int i; int j; int k; ideal lpp; list GPi2; list pairspP; ideal B; int ti;
  int i1; int tes; int j1; int selind; int i2; //int m;
  list prep; list crep; list LCU; poly p; poly lcp; ideal FF;
  list lpi;
  def RR=basering;
  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;
  GS=cgsdr(F,L); // "null",NW[1],"nonnull",NW[2],"cgs",CGS,"comment",comment);
  int start=timer;
  GP=GS;
  ideal lppr;
  list LL;
  list S;
  poly sp;
  for (i=1;i<=size(GP);i++)
  {
    kill LL;
    list LL;
    lpp=GP[i][1];
    GPi2=GP[i][2];
    lpph=GP[i][3];
    kill pairspP; list pairspP;
    for(j=1;j<=size(GPi2);j++)
    {
      pairspP[size(pairspP)+1]=GPi2[j][3];
    }
    LCU=LCUnion(pairspP);
    kill prep; list prep;
    kill crep; list crep;
    for(k=1;k<=size(LCU);k++)
    {
      prep[k]=list(LCU[k][2],LCU[k][3]);
      B=GPi2[LCU[k][1][1]][2]; // ATENTION last 1 has been changed to [2]
      LCU[k][1]=B;
    }
    //"Deciding if combine is needed";
    crep=PtoCrep(prep);
    if(size(LCU)>1){tes=1;}
    else
    {
      tes=0;
      for(k=1;k<=size(B);k++){B[k]=pnormalf(B[k],crep[1],crep[2]);}
    }
    if(tes)
    {
      // combine is needed
      kill B; ideal B;
      for (j=1;j<=size(LCU);j++)
      {
        LL[j]=LCU[j][2];
      }
      FF=precombine(LL);
      for (k=1;k<=size(lpp);k++)
      {
        kill L; list L;
        for (j=1;j<=size(LCU);j++)
        {
          L[j]=list(LCU[j][2],LCU[j][1][k]);
        }
        B[k]=combine(L,FF);
      }
    }
    for(j=1;j<=size(B);j++)
    {
      B[j]=pnormalf(B[j],crep[1],crep[2]);
    }
    S[i]=list(lpp,B,prep,crep,lpph);
    if(comment>=1)
    {
      lpi[size(lpi)+1]=string("[",i,"]");
      lpi[size(lpi)+1]=S[i][1];
    }
  }
  if(comment>=1)
  {
    string("Time in LCUnion + combine = ",timer-start);
    if(comment>=2){string("lpp=",lpi)};
  }
  return(S);
}

// grobcov
// input:
//    ideal F: a parametric ideal in Q[a][x], (a=parameters, x=variables).
//    list #: (options) list("null",N,"nonnull",W,"can",0-1,ext",0-1, "rep",0-1-2)
//            where
//            N is the null conditions ideal (if desired)
//            W is the ideal of non-null conditions (if desired)
//            The value of \"can\" is 1 by default and can be set to 0 if we do not
//            need to obtain the canonical GC, but only a GC.
//            The value of \"ext\" is 0 by default and so the generic representation
//             of the bases is given. It can be set to 1, and then the full
//             representation of the bases is given.
//            The value of \"rep\" is 0 by default, and then the segments
//            are given in canonical P-representation. It can be set to 1
//            and then they are given in canonical C-representation.
//            If it is set to 2, then both representations are given.
// output:
//    list S: ((lpp,basis,(idp_1,(idp_11,..,idp_1s_1))), ..
//             (lpp,basis,(idp_r,(idp_r1,..,idp_rs_r))) ) where
//            each element of S corresponds to a lpp-segment
//            given by the lpp, the basis, and the P-representation of the segment
proc grobcov(ideal F,list #)
"USAGE: grobcov(ideal F[,options]);
       F: ideal in Q[a][x] (a=parameters, x=variables) to be
       discussed.This is the fundamental routine of the
       library. It computes the Groebner Cover of a parametric
       ideal F in Q[a][x]. See
       A. Montes , M. Wibmer, \"Groebner Bases for Polynomial
       Systems with parameters\".
       JSC 45 (2010) 1391-1425.)
       or the not yet published book
       A. Montes. \" The Groebner Cover\" (Discussing
       Parametric Polynomial Systems).
       The Groebner Cover of a parametric ideal F consist
       of a set of pairs(S_i,B_i), where the S_i are disjoint
       locally closed segments of the parameter space,
       and the B_i are the reducedGroebner bases of the
       ideal on every point of S_i. The ideal F must be
       defined on a parametric ring Q[a][x] (a=parameters,
       x=variables).
RETURN: The list  [[lpp_1,basis_1,segment_1],  ...,
       [lpp_s,basis_s,segment_s]]
       optionally  [[ lpp_1,basis_1,segment_1,lpph_1],  ...,
       [lpp_s,basis_s,segment_s,lpph_s]]
       The lpp are constant over a segment and
       correspond to the set of lpp of the reduced
       Groebner basis for each point of the segment.
       With option (\"showhom\",1) the lpph will be
       shown: The lpph corresponds to the lpp of the
       homogenized ideal and is different for each
       segment. It is given as a string, and shown
       only for information. With the default option
       \"can\",1, the segments have different lpph.
       Basis: to each element of lpp corresponds
       an I-regular function given in full
       representation (by option (\"ext\",1)) or
       in generic representation (default option (\"ext\",0)).
       The I-regular function is the corresponding
       element of the reduced Groebner basis for
       each point of the segment with the given lpp.
       For each point in the segment, the polynomial
       or the set of polynomials  representing it,
       if they do not specialize to 0, then after
       normalization, specializes to the corresponding
       element of the reduced Groebner basis.
       In the full representation at least one of the
       polynomials representing the I-regular function
       specializes to non-zero.
       With the default option (\"rep\",0) the
       representation of the segment is the
       P-representation.
       With option (\"rep\",1) the representation
       of the segment is the C-representation.
       With option (\"rep\",2) both representations
       of the segment are given.
       The P-representation of a segment is of the form
       [[p_1,[p_11,..,p_1k1]],..,[p_r,[p_r1,..,p_rkr]]]
       representing the segment
       Union_i ( V(p_i) \ ( Union_j V(p_ij) ) ),
       where the p's are prime ideals.
       The C-representation of a segment is of the form
       (E,N) representing V(E) \ V(N), and the ideals E
       and N are radical and N contains E.
OPTIONS: An option is a pair of arguments: string,
       integer. To modify the default options, pairs
       of arguments -option name, value- of valid options
       must be added to the call.
       \"null\",ideal E: The default is (\"null\",ideal(0)).
       \"nonnull\",ideal N: The default is
       (\"nonnull\",ideal(1)).
       When options \"null\" and/or \"nonnull\" are given,
       then the parameter space is restricted to V(E) \ V(N).
       \"can\",0-1: The default is (\"can\",1).
       With the default option the homogenized
       ideal is computed before obtaining the Groebner
       Cover, so that the result is the canonical Groebner
       Cover. Setting (\"can\",0) only homogenizes the
       basis so the result is not exactly canonical,
       but the computation is shorter.
       \"ext\",0-1: The default is (\"ext\",0).
       With the default (\"ext\",0), only the generic
       representation of the bases is computed
       (single polynomials, but not specializing
       to non-zero for every point of the segment.
       With option (\"ext\",1) the full representation
       of the bases is computed (possible sheaves)
       and sometimes a simpler result is obtained,
       but the computation is more time consuming.
       \"rep\",0-1-2: The default is (\"rep\",0)
       and then the segments are given in canonical
       P-representation.
       Option (\"rep\",1) represents the segments
       in canonical C-representation, and
       option (\"rep\",2) gives both representations.
       \"comment\",0-3: The default is (\"comment\",0).
       Setting \"comment\" higher will provide
       information about the development of the
       computation.
       \"showhom\",0-1: The default is (\"showhom\",0).
       Setting \"showhom\",1 will output the set
       of lpp of the homogenized ideal of each segment
       as last element. One can give none or whatever
       of these options.
NOTE:    The basering R, must be of the form Q[a][x],
       (a=parameters, x=variables), and
       should be defined previously. The ideal
       must be defined on R.
KEYWORDS: Groebner cover; parametric ideal; canonical; discussion of parametric ideal
EXAMPLE:  grobcov; shows an example"
{
  list S; int i; int ish=1; list GBR; list BR; int j; int k;
  ideal idp; ideal idq; int s; ideal ext; list SS;
  ideal E; ideal N; int canop;  int extop; int repop;
  int comment=0; int m;
  def RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  Rx[1]=0;
  def D=ring(Rx);
  def RP=D+P;
  list L0=#;
  list Se;
  int out=0;
  int showhom=0;
  int hom;
  L0[size(L0)+1]="res"; L0[size(L0)+1]=ideal(1);
  // default options
  int start=timer;
  E=ideal(0);
  N=ideal(1);
  canop=1; // canop=0 for homogenizing the basis but not the ideal (not canonical)
           // canop=1 for working with the homogenized ideal
  repop=0; // repop=0 for representing the segments in Prep
           // repop=1 for representing the segments in Crep
           // repop=2 for representing the segments in Prep and Crep
  extop=0; // extop=0 if only generic representation of the bases are to be computed
           // extop=1 if the full representation of the bases are to be computed
  for(i=1;i<=size(L0) div 2;i++)
  {
    if(L0[2*i-1]=="can"){canop=L0[2*i];}
    else
    {
      if(L0[2*i-1]=="ext"){extop=L0[2*i];}
      else
      {
        if(L0[2*i-1]=="rep"){repop=L0[2*i];}
        else
        {
          if(L0[2*i-1]=="null"){E=L0[2*i];}
          else
          {
            if(L0[2*i-1]=="nonnull"){N=L0[2*i];}
            else
            {
              if (L0[2*i-1]=="comment"){comment=L0[2*i];}
              else
              {
                if (L0[2*i-1]=="showhom"){showhom=L0[2*i];}
              }
            }
          }
        }
      }
    }
  }
  if(not((canop==0) or (canop==1)))
  {
    string("Option can = ",canop," is not supported. It is changed to can = 1");
    canop=1;
  }
  for(i=1;i<=size(L0) div 2;i++)
  {
    if(L0[2*i-1]=="can"){L0[2*i]=canop;}
  }
  if ((printlevel) and (comment==0)){comment=printlevel;}
  list LL;
  LL[1]="can";     LL[2]=canop;
  LL[3]="comment"; LL[4]=comment;
  LL[5]="out";     LL[6]=0;
  LL[7]="null";    LL[8]=E;
  LL[9]="nonnull"; LL[10]=N;
  LL[11]="ext";    LL[12]=extop;
  LL[13]="rep";    LL[14]=repop;
  LL[15]="showhom";    LL[16]=showhom;
  if (comment>=1)
  {
    string("Begin grobcov with options: ",LL);
  }
  kill S;
  def S=gcover(F,LL);
  // NOW extendGC
  if(extop)
  {
    S=extendGC(S,LL);
  }
  // NOW repop and showhom
  list Si; list nS;
  for(i=1;i<=size(S);i++)
  {
    if(repop==0){Si=list(S[i][1],S[i][2],S[i][3]);}
    if(repop==1){Si=list(S[i][1],S[i][2],S[i][4]);}
    if(repop==2){Si=list(S[i][1],S[i][2],S[i][3],S[i][4]);}
    if(showhom==1){Si[size(Si)+1]=S[i][5];}
    nS[size(nS)+1]=Si;
  }
  S=nS;
  if (comment>=1)
  {
    string("Time in grobcov = ", timer-start);
    string("Number of segments of grobcov = ", size(S));
  }
  return(S);
}
example
{
"EXAMPLE:"; echo = 2;
// Casas conjecture for degree 4:
  if(defined(R)){kill R;}
  ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp;
  short=0;
  ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0),
            x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1),
            x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0),
            x2^2+(2*a3)*x2+(a2),
            x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0),
            x3+(a3);
  grobcov(F);

// EXAMPLE:
// M. Rychlik robot;
// Complexity and Applications of Parametric Algorithms of
//    Computational Algebraic Geometry.;
// In: Dynamics of Algorithms, R. de la Llave, L. Petzold and J. Lorenz eds.;
// IMA Volumes in Mathematics and its Applications,
//    Springer-Verlag 118: 1-29 (2000).;
//    (18. Mathematical robotics: Problem 4, two-arm robot)."
  if (defined(R)){kill R;}
  ring R=(0,a,b,l2,l3),(c3,s3,c1,s1), dp;
  short=0;
  ideal S12=a-l3*c3-l2*c1,b-l3*s3-l2*s1,c1^2+s1^2-1,c3^2+s3^2-1;
  S12;
  grobcov(S12);
}

// Auxiliary routine called by extendGC
// extendpoly
// input:
//   poly f: a generic polynomial in the basis
//   ideal idp: such that ideal(S)=idp
//   ideal idq: such that S=V(idp) \ V(idq)
////   NW the list of ((N1,W1),..,(Ns,Ws)) of red-rep of the grouped
////      segments in the lpp-segment  NO MORE USED
// output:
proc extendpoly(poly f, ideal idp, ideal idq)
"USAGE: extendGC(poly f,ideal p,ideal q);
       f is a polynomial in Q[a][x] in generic representation
       of an I-regular function F defined on the locally
       closed segment S=V(p) \ V(q).
       p,q are ideals in Q[a], representing the Crep of
       segment S.
RETURN: the extended representation of F in S.
       It can consist of a single polynomial or a set of
       polynomials when needed.
NOTE: The basering R, must be of the form Q[a][x],
       (a=parameters,x=variables), and should be
       defined previously. The ideals must be defined on R.
KEYWORDS: Groebner cover; parametric ideal; locally closed set;
       parametric ideal; generic representation; full representation;
       I-regular function
EXAMPLE:  extendpoly; shows an example"
{
  int te=0;
  def RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  Rx[1]=0;
  def D=ring(Rx);
  def RP=D+P;
  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);
    }
  }
}
example
{
  echo = 2; "EXAMPLE:";
  if(defined(R)){kill R;}
  ring R=(0,a1,a2),(x),lp;
  short=0;
  poly f=(a1^2-4*a1+a2^2-a2)*x+(a1^4-16*a1+a2^3-4*a2);
  ideal p=a1*a2;
  ideal q=a2^2-a2,a1*a2,a1^2-4*a1;
  extendpoly(f,p,q);" ";
// EXAMPLE;
if (defined(R)){kill R;}
ring R=(0,a0,b0,c0,a1,b1,c1,a2,b2,c2),(x), dp;
short=0;
poly f=(b1*a2*c2-c1*a2*b2)*x+(-a1*c2^2+b1*b2*c2+c1*a2*c2-c1*b2^2);
ideal p=
  (-a0*b1*c2+a0*c1*b2+b0*a1*c2-b0*c1*a2-c0*a1*b2+c0*b1*a2),
  (a1^2*c2^2-a1*b1*b2*c2-2*a1*c1*a2*c2+a1*c1*b2^2+b1^2*a2*c2-b1*c1*a2*b2+c1^2*a2^2),
  (a0*a1*c2^2-a0*b1*b2*c2-a0*c1*a2*c2+a0*c1*b2^2+b0*b1*a2*c2-b0*c1*a2*b2
  - c0*a1*a2*c2+c0*c1*a2^2),
  (a0^2*c2^2-a0*b0*b2*c2-2*a0*c0*a2*c2+a0*c0*b2^2+b0^2*a2*c2-b0*c0*a2*b2+c0^2*a2^2),
  (a0*a1*c1*c2-a0*b1^2*c2+a0*b1*c1*b2-a0*c1^2*a2+b0*a1*b1*c2-b0*a1*c1*b2
  -c0*a1^2*c2+c0*a1*c1*a2),
  (a0^2*c1*c2-a0*b0*b1*c2-a0*c0*a1*c2+a0*c0*b1*b2-a0*c0*c1*a2+b0^2*a1*c2
  -b0*c0*a1*b2+c0^2*a1*a2),
  (a0^2*c1^2-a0*b0*b1*c1-2*a0*c0*a1*c1+a0*c0*b1^2+b0^2*a1*c1-b0*c0*a1*b1+c0^2*a1^2),
  (2*a0*a1*b1*c1*c2-a0*a1*c1^2*b2-a0*b1^3*c2+a0*b1^2*c1*b2-a0*b1*c1^2*a2
  -b0*a1^2*c1*c2+b0*a1*b1^2*c2-b0*a1*b1*c1*b2+b0*a1*c1^2*a2-c0*a1^2*b1*c2+c0*a1^2*c1*b2);

ideal q=
  (-a1*c2+c1*a2),
  (-a1*b2+b1*a2),
  (-a0*c2+c0*a2),
  (-a0*b2+b0*a2),
  (-a0*c1+c0*a1),
  (-a0*b1+b0*a1),
  (-a1*b1*c2+a1*c1*b2),
  (-a0*b1*c2+a0*c1*b2),
  (-a0*b0*c2+a0*c0*b2),
  (-a0*b0*c1+a0*c0*b1);

extendpoly(f,p,q);
}

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

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

def FGC=extendGC(GCS,"rep",1);
  // 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;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  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,5)==0){J[size(J)+1]=i;}
  }
  for (i=1;i<=size(J);i++)
  {
    F=ideal(1);
    for (k=1;k<=size(P);k++)
    {
      if (k!=J[i])
      {
        F=idint(F,P[k]);
      }
    }
    f=f+F[1];
  }
  setring(RR);
  def fr=imap(P,f);
  return(fr);
}

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

// Auxiliary routine
// monoms
// Input: A polynomial f
// Output: The list of leading terms
static proc monoms(poly f)
{
  list L;
  poly lm; poly lc; poly lp; poly Q; poly mQ;
  def p=f;
  int i=1;
  while (p!=0)
  {
    lm=lead(p);
    p=p-lm;
    lc=leadcoef(lm);
    lp=leadmonom(lm);
    L[size(L)+1]=list(lc,lp);
    i++;
  }
  return(L);
}


// Auxiliary routine
// findindexpolys
// input:
//   list coefs=( (q11,..,q1r_1),..,(qs1,..,qsr_1) )
//               of denominators of the monoms
// output:
//   list ind=(v_1,..,v_t) of intvec
//        each intvec v=(i_1,..,is) corresponds to a polynomial in the sheaf
//        that will be built from it in extend procedures.
static proc findindexpolys(list coefs)
{
  int i; int j; intvec numdens;
  for(i=1;i<=size(coefs);i++)
  {
    numdens[i]=size(coefs[i]);
  }
//  def RR=basering;
//  def Rx=ringlist(RR);
//  def P=ring(Rx[1]);
//  setring(P);
//  def coefsp=imap(RR,coefs);
  def coefsp=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);
  //"T_combpolys="; combpolys;
  //setring(RR);
  //def combpolysT=imap(P,combpolys);
 // return(combpolysT);
 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 fP, poly fQ, ideal idp, ideal idq)
{
  def RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  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,fP);
  poly Q0=imap(RR,fQ);
  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 C, ideal N, ideal M)
{
  int numcombused;
//   def RR=basering;
//   def Rx=ringlist(RR);
//   def P=ring(Rx[1]);
//   setring(P);
//   def C=imap(RR,CC);
//   def N=imap(RR,NN);
//   def M=imap(RR,MM);
  if (ncols(C)==1){return(C);}

  int i; int j; int k; list c; intvec ci; intvec c0; intvec c1;
  list T; list T0; list T1; list LL; ideal N1;ideal M1; int te=0;
  for(i=1;i<=ncols(C);i++)
  {
    if((C[1,i]!=0) and (C[2,i]!=0))
    {
      if(c0==intvec(0)){c0[1]=i;}
      else{c0[size(c0)+1]=i;}
    }
  }
  def C1=submat(C,1..2,c0);
  for (i=1;i<=ncols(C1);i++)
  {
    c=comb(ncols(C1),i);
    for(j=1;j<=size(c);j++)
    {
      ci=c[j];
      numcombused++;
      if(i==1){N1=N+C1[2,j]; M1=M;}
      if(i>1)
      {
        kill c0; intvec c0 ; kill c1; intvec c1;
        c1=ci[size(ci)];
        for(k=1;k<size(ci);k++){c0[k]=ci[k];}
        T0=searchinlist(c0,LL);
        T1=searchinlist(c1,LL);
        N1=T0[1]+T1[1];
        M1=intersect(T0[2],T1[2]);
      }
      T=list(ci,PtoCrep0(Prep0(N1,M1)));
      LL[size(LL)+1]=T;
      if(equalideals(T[2][1],ideal(1))){te=1; break;}
    }
    if(te){break;}
  }
  ci=T[1];
  def Cs=submat(C1,1..2,ci);
//  setring(RR);
//  return(imap(P,Cs));
  return(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
//    chosen 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
// intersp: computes the intersection of the ideals in S in @P
static proc intersp(list S)
{
  def RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  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;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  setring(P);
  def fp=imap(RR,f);
  def idap=imap(RR,ida);
  poly @t;
  ring H=0,@t,dp;
  def PH=P+H;
  setring(PH);
  def fH=imap(P,fp);
  def idaH=imap(P,idap);
  idaH[size(idaH)+1]=1-@t*fH;
  option(redSB);
  def G=std(idaH);
  if (G==1){te=1;} else {te=0;}
  setring(RR);
  return(te);
}

// Auxiliary routine
// selectextendcoef
// input:
//    matrix CC: CC=(p_a1 .. p_ar_a)
//                  (q_a1 .. q_ar_a)
//            the matrix of elements of a coefficient in oo[a].
//    (ideal ida, ideal idb): the canonical representation of the segment S.
// output:
//    list caout
//            the minimum set of elements of CC needed such that at least one
//            of the q's is non-null on S, as well as the C-rep of of the
//            points where the q's are null on S.
//            The elements of caout are of the form (p,q,prep);
static proc selectextendcoef(matrix CC, ideal ida, ideal idb)
{
  def RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  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;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  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;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  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],5)!=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
//   Uses Rabinowiitz trick
static proc inconsistent(ideal E, ideal N)
{
  int j;
  int te=1;
  int tt;
  def RR=basering;
  def Rx=ringlist(RR);
  if(size(Rx[1])==4)
  {
    tt=1;
    def P=ring(Rx[1]);
    setring(P);
    def EP=imap(RR,E);
    def NP=imap(RR,N);
  }
  else
  {
    def EP=E;
    def NP=N;
  }
  poly @t;
  ring H=0,@t,dp;
  if(tt==1)
  {
    def RH=P+H;
   setring(RH);
   def EH=imap(P,EP);
   def NH=imap(P,NP);
  }
  else
  {
    def RH=RR+H;
    setring(RH);
    def EH=imap(RR,EP);
    def NH=imap(RR,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(RR);
  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 RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  Rx[1]=0;
  def D=ring(Rx);
  def RP=D+P;
  setring(RP);
  def HH=imap(RR,H);
  poly h=1;
  for (i=1;i<=size(HH);i++)
  {
    h=lcm(h,HH[i]);
  }
  setring(RR);
  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 #)
{
//   def RR=basering;
//   def Rx=ringlist(RR);
//   def P=ring(Rx[1]);
//   Rx[1]=0;
//   def D=ring(Rx);
//   def RP=D+p;
//   // 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);
    //"T_CG="; CG;
    if( size(CG)>0)
    {
      CG=groupKSWsegments(CG);
      if (comment>0)
      {
        string("Number of lpp segments = ",size(CG));
        string("Time in KSW + group + Prep = ",timer-start);
      }
    }
  }
  return(CG);
}

// Auxiliary routine
// sqf
static proc sqf(poly f)
{
  def RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  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 RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  Rx[1]=0;
  def D=ring(Rx);
  def RP=D+P;
  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(RR,E0);
  E1=std(E1);
  ideal N1=imap(RR,N0);
  N1=std(N1);
  setring(RR);
  E0=imap(P,E1);
  N0=imap(P,N1);
  if (inconsistent(E0,N0)==1)
  {
    return(emp);
  }
  setring(RP);
  def FRP=imap(RR,F);
  def ERP=imap(RR,E);
  FRP=FRP+ERP;
  option(redSB);
  def GRP=std(FRP);
  setring(RR);
  def G=imap(RP,GRP);
  if (memberpos(1,G)[1]==1)
  {
    if(comment>1){"Basis 1 is found"; E; N;}
    list KK; KK[1]=list(E0,N0,ideal(1));
    return(KK);
   }
  ideal Gr; ideal Gm; ideal GM;
  for (i=1;i<=size(G);i++)
  {
    if (variables(G[i])[1]==0){Gr[size(Gr)+1]=G[i];}
    else{Gm[size(Gm)+1]=G[i];}
  }
  ideal Gr0;
  for (i=1;i<=size(Gr);i++)
  {
    Gr0[i]=sqf(Gr[i]);
  }


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

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

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

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

//********************* End KapurSunWang *************************

//********************* Begin ConsLevels ***************************

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

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

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

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

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


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

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

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

def LL=ConsLevels(L);
LL;
def LLL=Levels(LL);
LLL;
}

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

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

def LL=ConsLevels(L);
LL;
def LLL=Levels(LL);
LLL;
}

proc DifConsLCSets(list A, list B)
"USAGE: DifConsLCSets(list A,list B);
       Input: The input lists A and B must be each one
       the canonical representations of the respective constructible sets,
       i.e. outputs of the routine ConsLevels for a constructible set,
       or from the routine Grob1Levels applied to the
       output of grobcov.
         A=[a1,a2,..,ak],
         B=[b1,b2,..,bj],
       where the a's and the b's are the closures
       of the levels of the constructible and the complements
       determined by ConsLevels (or GrobLevels)

       To be called in a ring Q[a][x]
       or a ring Q[a]. But the ideals can contain
       only the parameters in Q[a].
RETURN: A list of locally closed sets equivalent to the difference  S= A "\" B.
       Lc=[ [1][p1,q1]] [[2][p2,q2]]..],
       For obtaining the canonical representation into levels of
       the constructible A "\" B one have to apply ConsLevels and
       then optatively Levels.

KEYWORDS: constructible sets; canonical form
EXAMPLE:   DifConsLCSets shows an example"
{
  int n; int m; int t; int i; int j; int k;
  ideal ABup;
  ideal ABdw;
  if (size(B) mod 2==1){B[size(B)+1]=ideal(1);}
  if (size(A) mod 2==1){A[size(A)+1]=ideal(1);}
  //"T_A=";A;
 // "T_B="; B;
  n=size(A) div 2;
  m=(size(B) div 2)-1;
  //string("T_n=",n);
  //string("T_m=",m);
  list L;
  list M;
  list ABupC;
  //list LL;
  for(i=1;i<=n;i++)
  {
    //string("T_i=",i);
    t=1;
    j=0;
    //list L;
    while (t==1 and j<=m)
    {
      //string("T_j=",j);
      ABdw=intersectpar(list(A[2*i],B[2*j+1]));
      //"T_ABdw="; ABdw;
      ABup=A[2*i-1];
      //"T_ABup1="; ABup;
      if(j>0)
      {
        for(k=1;k<=size(B[2*j]);k++)
        {
          ABup[size(ABup)+1]=B[2*j][k];
        }
      }
      //"T_ABup2="; ABup;
      ABupC=Crep(ABup,ideal(1));
      //"T_ABupC="; ABupC;
      ABup=ABupC[1];
       //"T_ABup="; ABup;
      if(ABup==1){t=0;}
     //if(equalideals(ABup,ideal(1))){t=0;}
      else
      {
        M=Crep(ABup,ABdw);
        //"T_M="; M;
        //if(not(equalideals(M[1],ideal(1)))) {L[size(L)+1]=M;}
        if(not(size(M)==0)) {L[size(L)+1]=M;}
      }
      //"L="; L;
      j++;
    }
    //LL[size(LL)+1]=L;
  }
  return(L);
}
example
{
"EXAMPLE:"; echo = 2;
if(defined(R)){kill R;}
ring R=(0,x,y,z,t),(x1,y1),lp;
ideal a1=x;
ideal a2=x,y;
ideal a3=x,y,z;
ideal a4=x,y,z,t;

ideal b1=y;
ideal b2=y,z;
ideal b3=y,z,t;
ideal b4=1;

list L1=a1,a2,a3,a4;
list L2=b1,b2,b3,b4;

L1;
L2;

def LL=DifConsLCSets(L1,L2);
LL;
def LLL=ConsLevels(LL);
LLL;
def LLLL=Levels(LLL);
LLLL;
}

//**************************** End ConstrLevels ******************

//******************** Begin locus and envelop ******************************

// indepparameters
// Auxiliary routine to detect 'Special' components of the locus
// Input: ideal B
// Output:
//   1 if the ideal does not depend on the parameters
//   0 if they depend
static proc indepparameters(ideal B)
{
  def RR=basering;
  list Rx=ringlist(RR);
  def P=ring(Rx[1]);
  Rx[1]=0;
  def D=ring(Rx);
  def RP=D+P;
    // if(defined(@P)){kill @P; kill @RP; kill @R;}
   //  setglobalrings();
   ideal v=variables(B);
  setring RP;
  def BP=imap(RR,B);
  def vp=imap(RR,v);
  ideal varpar=variables(BP);
  int te;
  te=equalideals(vp,varpar);
  setring(RR);
  // kill @P; kill @RP; kill @R;
  if(te){return(1);}
  else{return(0);}
}

// indepparameterspoly
// Auxiliary routine to detect 'Special' components of the locus
// Input: ideal B
// Output:
//   1 if the solutions of the ideal (or poly) does not depend on the parameters
//   0 if they depend
static proc indepparameterspoly(B)
{
  def RR=basering;
  list Rx=ringlist(RR);
  def P=ring(Rx[1]);
  Rx[1]=0;
  def D=ring(Rx);
  def RP=D+P;
    // if(defined(@P)){kill @P; kill @RP; kill @R;}
   //  setglobalrings();
   ideal v=variables(B);
  setring RP;
  def BP=imap(RR,B);
  def vp=imap(RR,v);
  ideal varpar=variables(BP);
  int te;
  te=equalideals(vp,varpar);
  setring(RR);
  // kill @P; kill @RP; kill @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 returns 1
static proc dimP0(ideal N)
{
  def RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  setring P;
  // if(defined(@P)){ kill @P; kill @RP; kill @R;}
  // setglobalrings();
  // setring(@P);
  int te=1;
  def NP=imap(RR,N);
  attrib(NP,"IsSB",1);
  int d=dim(std(NP));
  //"T_d="; d;
  if(d==0){te=0;}
  setring(RR);
  return(te);
}

//  DimPar
//  Auxiliary 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 P=ring(RHx[1]);
   setring(P);
   def E2=std(imap(RRH,E));
   attrib(E2,"IsSB",1);
   def d=dim(E2);
   setring RRH;
   return(d);
 }

//  DimPar
//  Auxiliary routine to locus2 determining the dimension of a parametric ideal
//  it is identical to DimPar but adds information about the character of the component
static proc dimComp(ideal PA)
{
  def RR=basering;
  list Rx=ringlist(RR);
  int n=size(Rx[1][2]);
  def P=ring(Rx[1]);
  setring(P);
  list Lout;
  def CP=imap(RR,PA);
  attrib(CP,"IsSB",1);
  int d=dim(std(CP));
  if(d==n-1){Lout[1]=d;Lout[2]="Degenerate"; }
  else {Lout[1]=d; Lout[2]="Accumulation";}
  setring RR;
  return(Lout);
}

// 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 (anti-image).
//     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;
    ideal vmov;
    int version=0;
    int nv=nvars(RR);
    int moverdim=nv; //npars(RR);
    if(nv<4){version=1;}
    int d;
    poly F;
    for(i=1;i<=(size(DD) div 2);i++)
    {
      if(DD[2*i-1]=="moverdim"){moverdim=DD[2*i];}
      if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
      if(DD[2*i-1]=="movdim"){moverdim=DD[2*i];}
      if(DD[2*i-1]=="version"){version=DD[2*i];}
      if(DD[2*i-1]=="family"){F=DD[2*i];}
    }
    if(F!=0){Fenv=1;}
    list L0;
    if(dimP0(E)==0){L0=2,"Normal";} //  2 es fals pero ha de ser >0 encara que sigui 0
    else
    {
      if(version==0)
      {
          // Computing std(B+E,plex(x,y,x1,..xn)) one can detect if there is a first part
          // independent of parameters giving the variables with dimension 0
       dd=indepparameters(B);
       //"T_dd="; dd;
        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 antiimi=string(GGG);
        L0=d,antiimi;
        if(d==0)  // d<dim(V(E) \ V(N))?
        {
           " ";string("Anti-image of Special component = ", GGG);
        }
        else
        {
          L0[2]="Normal";
        }
      }
    }
     //"T_L0="; L0;
    return(L0);
  }

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

        //" ";string("Anti-image 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,5);
   //"T_AA="; AA; "T_FH="; FH; "T_EH="; EH; "T_rh="; rh;
   if(rh==0){env=1;} else{env=0;}
   setring RR;
          //L0[3]=env;
    //"T_env="; env;
    return(env);
}


// locus0(G): Private routine used by locus (the public routine), that
//               builds the different components.
// input:     The output G of the grobcov (in generic representation, which is the default option)
//       Options: The algorithm allows the following options as pair of arguments:
//               "moverdim", d  : by default moverdim is m=number of variables,
//                  but it can be set to less 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 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 dim @P, but as option it can be set to other values.
// output:
//         list, the canonical P-representation of the Normal and Non-Normal locus:
//              The Normal locus has two kind of components: Normal and Special.
//              The Non-normal locus has two kind of components: Accumulation and Degenerate.
//              This routine is complemented 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 RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  def RX=Rx;
  RX[1]=0;
  def D=ring(RX);
  def RP=D+P;
  //"T_locus0 ha estat cridat";
  // if(defined(@P)==1){te=1; kill @P; kill @R; kill @RP; }
  // setglobalrings();
  //Options
  list DD=#;
  ideal vmov;
  int dimpar=size(Rx[1][2]);
  int dimvar=size(Rx[2]);
  int nv=nvars(RR);
  int moverdim=nv; //size(ringlist(RR)[1][2]);
   for(i=1;i<=nv;i++)
  {
    //vmov[size(vmov)+1]=var(i+nv-moverdim);
    vmov[size(vmov)+1]=var(i);
  }
  //string("T_moverdim=",moverdim,"; dimvar=",dimvar,"; dimpar=",dimpar,"; vmov=",vmov);
  int version=0;
  //if(nv<4){version=1;}
  int comment=0;
  ideal Fm;
  poly F;
  for(i=1;i<=(size(DD) div 2);i++)
  {
    if(DD[2*i-1]=="moverdim"){moverdim=DD[2*i];}
    if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    if(DD[2*i-1]=="version"){version=DD[2*i];}
    if(DD[2*i-1]=="system"){Fm=DD[2*i];}
    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
    if(DD[2*i-1]=="family"){F=DD[2*i];}
  }
  list HHH;
  if (GG[1][1][1]==1 and GG[1][2][1]==1 and GG[1][3][1][1][1]==0 and GG[1][3][1][2][1]==1)
    {return(HHH);}
  list G1; list G2;
  def G=GG;
  list Q1; list Q2;
  int d; int j; int k;
  t1=1;
  for(i=1;i<=size(G);i++)
  {
    attrib(G[i][1],"IsSB",1);
    d=dim(std(G[i][1]));
    if(d==0){G1[size(G1)+1]=G[i];}
    else
    {
      if(d>0){G2[size(G2)+1]=G[i];}
    }
  }
  if(size(G1)==0){t1=0;}  // normal point components
  if(size(G2)==0){t2=0;}  // non-normal point components
  setring(RP);
  if(t1)
  {
    list G1RP=imap(RR,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]);
       }
       B=std(B);
       P1RP[size(P1RP)+1]=list(G1RP[i][3][k][1],G1RP[i][3][k][2],B);
     }
  }  //In P1RP the basis has been reduced wrt the top
   setring(RR);
  ideal h;
  list P1;
  if(t1)
  {
    P1=imap(RP,P1RP);
    for(i=1;i<=size(P1);i++)
    {
      for(j=1;j<=size(P1[i][3]);j++)
      {
        h=factorize(P1[i][3][j],1);
        P1[i][3][j]=h[1];
        for(k=2;k<=size(h);k++)
        {
          P1[i][3][j]=P1[i][3][j]*h[k];
        }
      }
    }
  } //In P1 factors in the basis independent of the parameters have been obtained
  ideal BB; int dd; list NS;
  for(i=1;i<=size(P1);i++)
  {
     //"T_P1="; P1;
     NS=NorSing(P1[i][3],P1[i][1],P1[i][2][1],DD);
     //"T_NS="; NS;
     dd=NS[1];
     if(dd==0){P1[i][3]=NS;}  //"Special"; //dd==0
     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;

  // if(defined(@P)){te=1; kill @P; kill @RP; kill @R;}
  // setglobalrings();

  setring(P);
  ideal J;
  if(t1==1)
  {
    def C1=imap(RR,P1);
    def L1=AddLocus(C1);
   }
  else{list C1; list L1; kill P1; list P1;}
  if(t2==1)
  {
    def C2=imap(RR,P2);
    def L2=AddLocus(C2);
  }
  else{list L2; list C2; kill P2; list P2;}
  for(i=1;i<=size(L2);i++)
  {
    J=std(L2[i][2]);
    d=dim(J);
    if(d+1==dimpar)
    {
      L2[i][4]=string("Degenerate",L2[i][4]);
    }
    else{L2[i][4]=string("Accumulation",L2[i][4]);}
  }

  list LN;
  if(t1==1)
  {
    for(i=1;i<=size(L1);i++){LN[size(LN)+1]=L1[i];}
  }
  if(t2==1)
  {
    for(i=1;i<=size(L2);i++){LN[size(LN)+1]=L2[i];}
  }
  int tLN=1;
  if(size(LN)==0){tLN=0;}
  setring(RR);
  if(tLN)
  {
   def L=imap(P,LN);
    for(i=1;i<=size(L);i++){if(size(L[i][2])==0){L[i][2]=ideal(1);}}
    //list LL;
    for(i=1;i<=size(L);i++)
    {
      if(typeof(L[i][4])=="list") {L[i][4][1]="Special";}
      l[1]=L[i][2];
      l[2]=L[i][3];
      l[3]=L[i][4];
      l[4]=L[i][5];
      L[i]=l;
    }
  }
  else{list L;}
  // if(te==1){kill @P; kill @RP; kill @R;}
  return(L);
}

// locus2(G): Private routine used by locus (the public routine), that
//                builds the different components.
// input:      The output G of the grobcov (in generic representation, which is the default option)
//                The ideal F defining the locus problem (G is the grobcov of F and has been
//                already determined by locus.
//       Options: The algorithm allows the following options as pair of arguments:
//                "moverdim", d  : by default moverdim is the number of variables but it can
//                 be set to minor 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 equal to the number of parameters of tracer variables, but as option
//                 it can be set to other values.
// output:
//         list, the canonical P-representation of the Normal and Non-Normal locus:
//              The Normal locus has two kind of components: Normal and Special.
//              The Non-normal locus has two kind of components: Accumulation and Degenerate.
//              This routine is compemented by locus that calls it in order to eliminate problems
//              with degenerate points of the mover.
//         The output components are given as
//              ((p1,(p11,..p1s_1),type_1,level_1),..,(pk,(pk1,..pks_k),type_k,level_k)
//         The components are given in canonical P-representation of the subset.
//              If all levels of a class of locus are 1, then the set is locally closed. Otherwise the level
//              gives the depth of the component.
static proc locus2(list G, ideal F, list #)
{
 int st=timer;
 list Snor; list Snonor;
 int d; int i; int j;
 def RR=basering;
 def Rx=ringlist(RR);
 def RP=ring(Rx[1]);
 int tax=1;
 list DD=#;
 list GG=G;
 int n=size(Rx[1][2]);
 int moverdim;
 for(i=1;i<=(size(DD) div 2);i++)
  {
     if(DD[2*i-1]=="moverdim"){moverdim=DD[2*i];}
     if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
//    if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
//    if(DD[2*i-1]=="version"){version=DD[2*i];}
//    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
//    if(DD[2*i-1]=="grobcov"){GG=DD[2*i];}
    if(DD[2*i-1]=="anti-image"){tax=DD[2*i];}
  }
 for(i=1;i<=size(GG);i++)
 {
   attrib(GG[i][1],"IsSB",1);
   GG[i][1]=std(GG[i][1]);
   d=dim(GG[i][1]);
   if(d==0)
   {
     for(j=1;j<=size(GG[i][3]);j++)
     {
       Snor[size(Snor)+1]=GG[i][3][j];
     }
   }
   else
   {
     if(d>0)
     {
       for(j=1;j<=size(GG[i][3]);j++)
       {
         Snonor[size(Snonor)+1]=GG[i][3][j];
      }
     }
   }
 }
 //"T_Snor="; Snor;
 //"T_Snonor="; Snonor;
 int tnor=size(Snor); int tnonor=size(Snonor);
 setring RP;
 list SnorP;
 list SnonorP;
 if(tnor)
 {
   SnorP=imap(RR,Snor);
   st=timer;
   SnorP=LCUnionN(SnorP);
 }
 if(tnonor)
 {
   SnonorP=imap(RR,Snonor);
   SnonorP=LCUnionN(SnonorP);
 }
 //"T_SnorP after LCUnion="; SnorP;
// "T_SnonorP  after LCUnion="; SnonorP;
 setring RR;
 ideal C;  list N;  list BAC; list AI;
 list NSC; list DAC;
 list L;
 ideal B;
 int k;
 int j0; int k0; int te;
 poly kkk=1;
 ideal AI0;
 int dimP;

 if(tnor)
 {
   Snor=imap(RP,SnorP);
   for(i=1;i<=size(Snor);i++)
   {
     C=Snor[i][1];
     N=Snor[i][2];
     dimP=DimPar(C);
      //"T_G="; G;
     AI=NS(F,G,C,N,DD);
     //"T_AI="; AI;
     AI0=AI[2];
     if(size(AI)==3){AI[2]="normal"; AI[3]=ideal(0);}
     else
     {
     if(AI[1]>=dimP)
     {
       AI[2]="Normal";
       if(tax>0){AI[3]=AI0;}
     }
     else
        {
           //"T_Al=";
           if(AI[1]<n-1){AI[2]="Special"; AI[3]=AI0;}
        }
     }
     Snor[i][size(Snor[i])+1]=AI;
   }
   for(i=1;i<=size(Snor);i++)
   {
     L[size(L)+1]=Snor[i];
   }
  }
  ideal AINN;
 if(tnonor)
 {
   Snonor=imap(RP,SnonorP);
   //"T_Snonor="; Snonor;
   //"T_G="; G;
   for(i=1;i<=size(Snonor);i++)
   {
     DAC=dimComp(Snonor[i][1]);
     Snonor[i][size(Snonor[i])+1]=DAC;
     if(tax==1)
     {
       AINN=AISnonor(Snonor[i][1],GG);
       Snonor[i][3][3]=AINN;
      // AQUI ES ON AFEGIM LA ANTIMATGE DE LES NONOR
     }
   }
   for(i=1;i<=size(Snonor);i++)
   {
     L[size(L)+1]=Snonor[i];
   }
 }
 return(L);
}

// Auxiliary algorithm of locus2.
//           The algorithm searches the basis corresponding to C, in the grobcov.
//           It reduces the basis modulo the component.
//           The result is the reduced basis BR.
//           For each hole of the component
//              it searches the segment where the hole is included
//              and select form its basis the polynomials
//              only dependent on the variables.
//              These polynomials are non-null in an open set of
//              the component, and are included in the list NoNul of non-null factors
// input: F: the ideal of the locus problem
//           G the grobcov of F
//           C the top of a component of normal points
//           N the holes of the component
// output: (d,tax,a)
//           where d is the dimension of the anti-image
//           a is the anti-image of the component and
//           tax is the taxonomy \"Normal\" if d is equal to the dimension of C
//           and \"Special\" if it is smaller.
//           When a normal point component has degree greater than 9, then the
//           taxonomy is not determined, and (n,'normal', 0) is returned as third
//           element of the component. (n is the dimension of the space).
static proc NS(ideal F,list G, ideal C, list N, list #)
{
  // Initializing and defining rings
  int i; int j; int k; int te; int j0;int k0; int m;
  list DD=#;
  def RR=basering;
  def Rx=ringlist(RR);
  def Lx=Rx;
  def P=ring(Rx[1]);               // ring Q[bx]
  Lx[1]=0;
  def D=ring(Lx);                   // ring Q[bu]
  def PR0=P+D;
  def PRx=ringlist(PR0);
  PRx[3][2][1]="lp";
  def PR=ring(PRx);               // ring Q[bx,bu]  in lex order
  int n=size(Rx[1][2]);           // number of parameters
  def nv=nvars(RR);               // number of variables
  int moverdim=nv;               // number of mover variables
  ideal vmov;                         // mover variables   (bw)
  for(i=1;i<=(size(DD) div 2);i++)
  {
    if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    if(DD[2*i-1]=="moverdim"){moverdim=DD[2*i];}
//    if(DD[2*i-1]=="version"){version=DD[2*i];}
//    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
//    if(DD[2*i-1]=="grobcov"){GG=DD[2*i];}
//    if(DD[2*i-1]=="anti-image"){tax=DD[2*i];}
  }

  //string("T_moverdim= ", moverdim);
  // vmov = the last vmov variables = mover variables
  for(i=1;i<=moverdim;i++)
  {
    vmov[size(vmov)+1]=var(i+nv-moverdim);
  }
  //string("T_nv=",nv,"  moverdim=",moverdim);
//  "T_aixo es nou";
//  for(i=1;i<=nv;i++)
//  {
//    vmov[size(vmov)+1]=var(i);
//  }

  int ddeg; int dp;
  list LK;
  ideal bu;                               // ideal of all variables   (u)
  for(i=1;i<=nv;i++){bu[i]=var(i);}
  ideal mv;                            // ideal of mover variables
  for(i=1;i<=nv;i++){mv[size(mv)+1]=var(i);}

   // Searching the basis associated to C
   j=2; te=1;
   while((te) and (j<=size(G)))
   {
      k=1;
      while((te) and (k<=size(G[j][3])))
      {
        if (equalideals(C,G[j][3][k][1])){j0=j; k0=k; te=0;}
        k++;
      }
      j++;
   }
   if(te==1){"ERROR";}
   def B=G[j0][2];  // Aixo aniria be per les nonor
   //"T_B=G[j0][2]="; B;
   //string("T_k0=",k0," G[",j0,"]="); G[j0];

   // Searching the elements in Q[v_m]  on basis different from B that are nul there
   // and cannot become 0 on the antiimage of B. They are placed on NoNul
   list NoNul;
   ideal BNoNul;
   ideal covertop;           // basis of the segment where a hole of C is the top
   int te1;
   for(i=1;i<=size(N);i++)
   {
     j=2; te=1;
     while(te and j<=size(G))
     {
       if(j!=j0)
       {
         k=1;
         while(te and k<=size(G[j][3]))
         {
           covertop=G[j][3][k][1];
           if(equalideals(covertop,N[i]))
           {
             te=0; te1=1; BNoNul=G[j][2];
            }
           else
           {
             if(redPbasis(covertop,N[i]))
             {
               te=0; te1=1; m=1;
               while( te1 and m<=size(G[j][3][k][2]) )
               {
                 if(equalideals(G[j][3][k][2][m] ,N[i] )==1){te1=0;}
                 m++;
               }
             }
             if(te1==1){ BNoNul=G[j][2];}
           }
           k++;
         }

         if((te==0) and (te1==1))
         {
          // Selecting the elements independent of the parameters,
          // They will be non null on the segment
           for(m=1;m<=size(BNoNul);m++)
           {
             if(indepparameterspoly(BNoNul[m]))
             {
                NoNul[size(NoNul)+1]=BNoNul[m];
             }
           }
         }
       }
       j++;
    }
  }

   // Adding F to B
   for(i=1;i<=size(F);i++)
   {
     B[size(B)+1]=F[i];
   }

  def E=NoNul;
  poly kkk=1;
  if(size(E)==0){E[1]=kkk;}

  // Avoiding computations that are too expensive for obtaining
  //   the anti-image of normal point components
  setring(P);
  def CP=imap(RR,C);
  ddeg=deg(CP);
  setring(RR);
  // if(n+nv>10 or ddeg>=10){LK=n,ideal(0),"normal"; return(LK);}
  if(ddeg>=10){LK=n,ideal(0),"normal"; return(LK);}  // 8 instead of 10 ?

  // Reducing  basis B modulo C  in the ring PR      lex(x,u)
  setring(PR);
  def buPR=imap(RR,bu);
  def EE=imap(RR,E);
  ideal BR;
  def BB=imap(RR,B);
  def CC=imap(RR,C);
  attrib(CC,"IsSB",1);
  CC=std(CC);
  for(i=1;i<=size(BB);i++)
  {
    BR[i]=reduce(BB[i],CC);
  }

  // Computing the GB(B) in the ring PR     lex(x,u)
  BR=std(BR);
  ideal BRP;

 // Eliminating the polynomials in B containing variables v_t
  for(i=1;i<=size(BR);i++)
  {
    if(subset(variables(BR[i]),buPR )){BRP[size(BRP)+1]=BR[i];}
  }
  BR=BRP;

  // Factoring and eliminating the factors in N
  list H;
  ideal BP;
  poly f;
  int tj;
  //ideal vv=imap(RR,bu);
  for(i=1;i<=size(BR);i++)
  {
    f=1;
    H=factorize(BR[i]);
    if((subset(variables(H[1][2]),buPR)) and (size(H[1]))==2){BP[size(BP)+1]=H[1][2];}
    else
    {
      for(j=1;j<=size(H[1]);j++)
      {
        tj=subset(H[1][j],EE);
        if((subset(variables(H[1][j]),buPR)) and (tj==0)){f=f*H[1][j];}
      }
      if(leadcoef(f)!=f){BP[size(BP)+1]=f;}
   }
  }

  // Determining the GB basis of B in Q[u]
  BP=std(BP);
  // Determining the dimension of the anti-image
  //   in the ring D
  setring D;
  def KK=imap(PR,BP);
  KK=std(KK);
  ideal KKM;
  ideal vmovM=imap(RR,vmov);

  // selecting the element in Q[w]
  for(i=1;i<=size(KK);i++)
  {
    if(subset(variables(KK[i]),vmovM))
    {
      KKM[size(KKM)+1]=KK[i];
    }
  }

  // Determining the dimensions
  list AIM=dimM(KKM,nv,moverdim);
  //int d2=AIM[1];
  def AAIM=AIM[2];
  setring(RR);
  //def BBR=imap(D,KK);
  def LKK=imap(D,AIM);
  //LKK=d2,string(BBR);
  return(LKK);
}

// Auxiliary algorithm of locus2.
//           The algorithm searches the basis corresponding to C, in the grobcov.
//           It reduces the basis modulo the component.
//           The result is the reduced basis BR.
//           for each hole of the component
//              it searches the segment where the hole is included
//              and select form its basis the polynomials
//              only dependent on the variables.
//              These polynomials are non-null in an open set of
//              the component, and are included in the list NoNul of non-null factors
// input: C: The top of the the segment of a non-normal locus
//           G the grobcov of F
// output: B: the anti-image of the input segment
static proc AISnonor(ideal C,list G, list #)
{
  // Initializing and defining rings
  int i; int j; int k; int te; int j0;int k0; int m;
  list DD=#;
  def RR=basering;
  def Rx=ringlist(RR);
  def Lx=Rx;
  def P=ring(Rx[1]);               // ring Q[bx]
  Lx[1]=0;
  def D=ring(Lx);                   // ring Q[bu]
  def PR0=P+D;
  def PRx=ringlist(PR0);
  PRx[3][2][1]="lp";
  def PR=ring(PRx);               // ring Q[bx,bu]  in lex order
  int n=size(Rx[1][2]);           // number of parameters
  def nv=nvars(RR);               // number of variables
  int moverdim=nv;                     // number of mover variables
  ideal vmov;                         // mover variables   (bw)

  for(i=1;i<=(size(DD) div 2);i++)
  {
    if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    if(DD[2*i-1]=="moverdim"){moverdim=DD[2*i];}
//    if(DD[2*i-1]=="version"){version=DD[2*i];}
//    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
//    if(DD[2*i-1]=="grobcov"){GG=DD[2*i];}
//    if(DD[2*i-1]=="anti-image"){tax=DD[2*i];}
  }
  for(i=1;i<=moverdim;i++)
  {
    vmov[size(vmov)+1]=var(i+nv-moverdim);
  }

  int ddeg; int dp;
  list LK;
  ideal bu;                               // ideal of all variables   (u)
  for(i=1;i<=nv;i++){bu[i]=var(i);}
  ideal mv;                            // ideal of mover variables
  for(i=1;i<=nv;i++){mv[size(mv)+1]=var(i);}

   // Searching the basis associated to Nonor
   j=2; te=1;
   while((te) and (j<=size(G)))
   {
      k=1;
      while((te) and (k<=size(G[j][3])))
      {
        if (equalideals(C,G[j][3][k][1])){j0=j; k0=k; te=0;}
        k++;
      }
      j++;
   }
   if(te==1){"ERROR";}
   def B=G[j0][2];
   //"T_B=G[j0][2]="; B;
   //string("T_k0=",k0," G[",j0,"]="); G[j0];
   return(B);
}

// NOT USED
// Auxiliary algorithm of locus2.
//           The algorithm searches the basis corresponding to C, in the grobcov.
//           It reduces the basis modulo the component.
//           The result is the reduced basis BR.
//           for each hole of the component
//              it searches the segment where the hole is included
//              and select form its basis the polynomials
//              only dependent on the variables.
//              These polynomials are non-null in an open set of
//              the component, and are included in the list NoNul of non-null factors
// input: F: the ideal of the locus problem
//           G the grobcov of F
//           C the top of a component of normal points
//           N the holes of the component
// output: (d,tax,a)
//           where d is the dimension of the anti-image
//           a is the anti-image of the component and
//           tax is the taxonomy \"Normal\" if d is equal to the dimension of C
//           and \"Special\" if it is smaller.
static proc NS5(ideal F,list G, ideal C, list N)
{
  // Initializing and defining rings
  int i; int j; int k; int te; int j0;int k0; int m;
  def RR=basering;
  def Rx=ringlist(RR);
  def Lx=Rx;
  def P=ring(Rx[1]);
  Lx[1]=0;
  def D=ring(Lx);
  def RP=D+P;
  def PR0=P+D;
  def PRx=ringlist(PR0);
  PRx[3][2][1]="lp";
  def PR=ring(PRx);
  int n=size(Rx[1][2]);                                     // dimension of the tracer space
  def nv=nvars(RR);                                        // number of variables
  ideal v; for(i=1;i<=nv;i++){v[i]=var(i);}        // variables
  int nm;                                                        // number of mover variables
  ideal vmov;                                                  // mover variables
  if(nv>=n){nm=n;}
  else{nm=nv;}
  for(i=1;i<=nm;i++)
  {
    vmov[size(vmov)+1]=var(i+nv-nm);
  }

  // Searching the basis associated to C
   j=2; te=1;
   while((te) and (j<=size(G)))
   {
      k=1;
      while((te) and (k<=size(G[j][3])))
      {
        if (equalideals(C,G[j][3][k][1])){j0=j; k0=k; te=0;}
        k++;
      }
      j++;
   }
   if(te==1){"ERROR";}
   def B=G[j0][2];     // basis associated to C
   //"T_B="; B;

  // Searching the elements in Q[vmov]  on bases different from B that are nul there
  // and cannot become 0 on the antiimage of C. They are placed on NoNul
  list NoNul;          // elements in Q[vmov]  on bases different from B
  ideal BNoNul;
  ideal covertop;   // basis of the segment whose top is a hole of our segment
  int te1;
  for(i=1;i<=size(N);i++)
  {
    j=2; te=1;
    while(te and j<=size(G))
    {
      if(j!=j0)
      {
        k=1;
        while(te and k<=size(G[j][3]))
        {
          covertop=G[j][3][k][1];
          if(equalideals(covertop,N[i]))
          {
            te=0; te1=1; BNoNul=G[j][2];
          }
          else
          {
            if(redPbasis(covertop,N[i]))
            {
              te=0; te1=1; m=1;
              while( te1 and m<=size(G[j][3][k][2]) )
              {
                if(equalideals(G[j][3][k][2][m] ,N[i] )==1){te1=0;}
                m++;
              }
            }
            if(te1==1){ BNoNul=G[j][2];}
          }
          k++;
        }

        if((te==0) and (te1==1))
        {
          // Selecting the elements independent of the parameters,
          // They will be non null on the segment
          for(m=1;m<=size(BNoNul);m++)
          {
            if(indepparameterspoly(BNoNul[m]))
            {
               NoNul[size(NoNul)+1]=BNoNul[m];
            }
         }
       }
     }
     j++;
   }
 }
  def NN=NoNul;
  poly kkk=1;
  if(size(NN)==0){NN[1]=kkk;}
  //"T_NN="; NN;

 // Union of F and B
 for(i=1;i<=size(F);i++)
 {
   B[size(B)+1]=F[i];
 }
 //"T_B+F="; B;

  // Reducing  basis B modulo C  in the ring PR
  setring(PR);
  def vv=imap(RR,v);
  ideal BBR;
  def BB=imap(RR,B);
  def CC=imap(RR,C);
  attrib(CC,"IsSB",1);
  CC=std(CC);
  for(i=1;i<=size(BB);i++)
  {
    BBR[i]=reduce(BB[i],CC);
  }
  //"T_BBR="; BBR;

  // Selecting the elements of the ideal in Q[v]
  // setring RP;
  ideal Bv;
  // def Bv=imap(PR,BBR);
  //ideal vvv=imap(RR,v);
  // Bv=std(Bv);
  ideal BR;
  for(i=1;i<=size(BBR);i++)
  {
    if(subset(variables(BBR[i]),vv)){BR[size(BR)+1]=BBR[i];}
  }
  //"T_BR="; BR;
  // setring PR;
  // def BR=imap(RP,BRv);
  def NNN=imap(RR,NN);

  // Factoring and eliminating the factors in N
  list H;
  ideal BP;
  poly f;
  int tj;
  // ideal vv=imap(RR,v);
  for(i=1;i<=size(BR);i++)
  {
    f=1;
    H=factorize(BR[i]);
    if((subset(variables(H[1][2]),vv)) and (size(H[1]))==2){BP[size(BP)+1]=H[1][2];}
    else
    {
      for(j=1;j<=size(H[1]);j++)
      {
        tj=subset(H[1][j],NNN);
        if((subset(variables(H[1][j]),vv)) and (tj==0)){f=f*H[1][j];}
      }
      if(leadcoef(f)!=f){BP[size(BP)+1]=f;}
   }
  }
  //"T_BP="; BP;

  //  Obtaining Am, the anti-image of C on vmov
  setring D;
  ideal A=imap(PR,BP);
  ideal vm=imap(RR,vmov);
  ideal Am;
  for(i=1;i<=size(A);i++)
  {
    if (subset(variables(A[i]),vm)){Am[size(Am)+1]=A[i];}
  }
  //"T_Am="; Am;


  list AIM=dimM(Am,nv,nm);
  setring(RR);
  def LAM=imap(D,AIM);
  return(LAM);
}

static proc dimM(ideal KKM, int nv, int nm)
{
  def RR=basering;
  list L;
  int i;
  def Rx=ringlist(RR);
  for(i=1;i<=nm;i++)
  {
    L[i]=Rx[2][nv-nm+i];
  }
  Rx[2]=L;
  intvec iv;
  for(i=1;i<=nm;i++){iv[i]=1;}
  Rx[3][1][2]=iv;
   def DM=ring(Rx);
  //"Rx="; Rx;
  setring(DM);
  ideal KKMD=imap(RR,KKM);
  attrib(KKMD,"IsSB",1);
  KKMD=std(KKMD);
  int d=dim(KKMD);
  setring(RR);
  def KAIM=imap(DM,KKMD);
  list LAIM=d,KAIM;
  return(LAIM);
}

// Procedure using only standard GB in lex(x,a) order to obtain the
//    component of the locus.
//    It is not so fine as locus and cannot evaluate the taxonomy, but
//    it is much simpler and efficient.
// input: ideal S for determining the locus
// output: the irreducible components of the locus
//    Data must be given in Q[a][x]
proc stdlocus(ideal F)
"USAGE: stdlocus(ideal F)
       The input ideal must be the set equations defining the locus.
       Calling sequence: locus(F);
       The input ring must be a parametrical ideal in Q[x][u],
       (x=tracer variables, u=remaining variables).
       (Inverts the concept of parameters and variables of the ring).
       Special routine for determining the locus of points of  a geometrical construction.
       Given a parametric ideal F representing the system determining the locus of points (x)
       which verify certain properties, the call to stdlocus(F)
       determines the different irreducible components of the locus.
       This is a simple routine, using only standard Groebner basis computation,
       elimination and prime decomposition instead of using grobcov.
       It does not determine the taxonomy, nor the holes of the components
RETURN:The output is a list of the tops of the components [C_1, .. , C_n] of the locus.
       Each component is given its top ideal p_i.
NOTE: The input must be the locus system.
KEYWORDS: geometrical locus; locus
EXAMPLE: stdlocus; shows an example"
{
  int i; int te;
  def RR=basering;
  list Rx=ringlist(RR);
  int n=npars(RR);  // size(Rx[1][2]);
  int nv=nvars(RR);
  ideal vpar;
  ideal vvar;
  //"T_n="; n;
   //"T_nv="; nv;
  for(i=1;i<=n;i++){vpar[size(vpar)+1]=par(i);}
  for(i=1;i<=nv;i++){vvar[size(vvar)+1]=var(i);}
  //string("T_vpar = ", vpar," vvar = ",vvar);
  def P=ring(Rx[1]);
  Rx[1]=0;
  def D=ring(Rx);
  def RP=D+P;
  list Lx=ringlist(RP);
  setring(RP);
  def FF=imap(RR,F);
  def vvpar=imap(RR,vpar);
  //string("T_vvpar = ",vvpar);
  ideal B=std(FF);
  //"T_B="; B;
  ideal Bel;
  //"T_vvpar="; vvpar;
  for(i=1;i<=size(B);i++)
  {
    if(subset(variables(B[i]),vvpar)) {Bel[size(Bel)+1]=B[i];}
  }
  //"T_Bel="; Bel;
  list H;
  list FH;
  H=minAssGTZ(Bel);
  int t1;
  if(size(H)==0){t1=1;}
  setring RR;
  list empt;
  if(t1==1){return(empt);}
  else
  {
    def HH=imap(RP,H);
    return(HH);
  }
}
example
{
  "EXAMPLE:"; echo = 2;
if(defined(R)){kill R;}
ring R=(0,x,y),(x1,y1),dp;
short=0;

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

stdlocus(S96);
}

//  locus(F):  Special routine for determining the locus of points
//                 of  geometrical constructions.
//  input:      The ideal of the locus equations
//  output:
//          The output components are given as
//               ((p1,(p11,..p1s_1),tax_1),..,(pk,(pk1,..pks_k),tax_k)
//               Elements 1 and 2 represent the P-canonical form of the component.
//               The third element tax is:
//                 for normal point components, tax=(d,taxonomy,anti-image)
//                    being d=dimension of the anti-image on the mover variables,
//                           taxonomy='Normal'  or  'Special', and
//                           anti-image=ideal of the anti-image over the mover variables
//                                which are taken to be the last n variables.
//                 for non-normal point components, tax =(d,taxonomy)
//                    being d=dimension of the component  and
//                           taxonomy='Accumulation' or 'Degenerate'.
//          The components are given in canonical P-representation of the subset.         l
//          The normal locus has two kind of components: Normal and Special.
//            Normal component:
//            - each point in the component has 0-dimensional anti-image.
//            - the anti-image in the mover coordinates is equal to the dimension of the component
//            Special component:
//            - each point in the component has 0-dimensional anti-image.
//            - the anti-image in the mover coordinates is smaller than the dimension of the component
//          The non-normal locus has two kind of components: Accumulation and Degenerate.
//           Accumulation points:
//             - each point in the component has anti-image of dimension greater than 0.
//             - the component has dimension less than n-1.
//           Degenerate components:
//             - each point in the component has anti-image of dimension greater than 0.
//             - the component has dimension n-1.
//           When a normal point component has degree greater than 9, then the
//           taxonomy is not determined, and (n,'normal', 0) is returned as third
//           element of the component. (n is the dimension of the space).

proc locus(ideal F, list #)
"USAGE:  locus(ideal F[,options])
        Special routine for determining the locus of points of
        a geometrical construction.
INPUT:  The input ideal must be the set equations
        defining the locus. Calling sequence:
        locus(F[,options]);
        The input ring must be a parametrical ideal
        in Q[x1,..,xn][u1,..,um],
        (x=tracer variables, u=mover variables).
        The number of tracer variables is n=dim(space).
        The mover variables can be restricted by the user to the last
        k variaboes by adding the option \"moverdim\",k
        (Inverts the concept of parameters and variables of
        the ring).
OPTIONS: An option is a pair of arguments: string, integer.
        To modify the default options, pairs of arguments
        -option name, value- of valid options must be added to
        the call.The algorithm allows the following options as
        pair of arguments:

        \"moverdim\", k  to restrict the mover-variables of the
        anti-image to the last k variables.
        By default k is equal to the number of u-variables,
        The restriction to a smaller value
        can produce an error in the character
         \"Normal\" or \"Special\" of a locus component.

        \"grobcov\", list of the previous computed grobcov(F).
        It is to be used when we modify externally the grobcov,
        for example to obtain the real grobcov.

        \"anti-image\",0  if the anti-image is not to be
        shown for \"Normal\" components.

        \"comments\", c: by default it is 0, but it can be set
        to 1.
RETURN: The output is a list of the components:
        ((p1,(p11,..p1s_1),tax_1),..,(pk,(pk1,..pks_k),tax_k)
        Elements 1 and 2 of a component represent the
        P-canonical form of the component.
        The third element tax is:
          for normal point components,
            tax=(d,taxonomy,anti-image) being
             d=dimension of the anti-image on the mover variables,
             taxonomy=\"Normal\"  or  \"Special\" and
             anti-image=ideal of the anti-image over the mover
             variables.
          for non-normal point components,
            tax =(d,taxonomy) being
             d=dimension of the component  and
             taxonomy=\"Accumulation\" or \"Degenerate\".
        The components are given in canonical P-representation.
        The normal locus has two kind of components:
          Normal and Special.
          Normal component:
           - each point in the component has 0-dimensional
              anti-image.
           - the anti-image in the mover coordinates is equal
              to the dimension of the component
          Special component:
           - each point in the component has 0-dimensional
              anti-image.
           - the anti-image in the mover coordinates is smaller
              than the dimension of the component
        The non-normal locus has two kind of components:
          Accumulation and Degenerate.
          Accumulation component:
           - each point in the component has anti-image of
              dimension greater than 0.
           - the component has dimension less than n-1.
         Degenerate components:
           - each point in the component has anti-image
              of dimension greater than 0.
           - the component has dimension n-1.
       When a normal point component has degree greater than 9,
         then the taxonomy is not determined, and (n,'normal', 0)
         is returned as third element of the component. (n is the
         dimension of the space).

       Given a parametric ideal F representing the system F
       determining the locus of points (x) which verify certain
       properties, the call to locus(F) determines the different
       classes of locus components, following the taxonomy
       defined in the book:
       A. Montes. \"The Groebner Cover\"
       A previous paper gives particular definitions
       for loci in 2d.
       M. Abanades, F. Botana, A. Montes, T. Recio,
       \"An Algebraic Taxonomy for Locus Computation
       in Dynamic Geometry\",
       Computer-Aided Design 56 (2014) 22-33.
NOTE: The input must be the locus system.
KEYWORDS: geometrical locus; locus; dynamic geometry
EXAMPLE: locus; shows an example"
{
  int tes=0; int i;  int m; int mm; int n;
  def RR=basering;
  list GG;
  //Options
  list DD=#;
  ideal vmov;
  int nv=nvars(RR);                                  // number of variables
  int moverdim=nv; //size(ringlist(RR)[1][2]);   // number of parameters
  if(moverdim>nv){moverdim=nv;}
  // for(i=1;i<=nv;i++){vmov[size(vmov)+1]=var(i);}
  int version=2;
  // if(nv<4){version=2;}
  int comment=0;
  int tax=1;
  ideal Fm;
  for(i=1;i<=(size(DD) div 2);i++)
  {
    if(DD[2*i-1]=="moverdim"){moverdim=DD[2*i];}
    if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    if(DD[2*i-1]=="version"){version=DD[2*i];}
    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
    if(DD[2*i-1]=="grobcov"){GG=DD[2*i];}
    if(DD[2*i-1]=="anti-image"){tax=DD[2*i];}
  }
  for(i=1;i<=moverdim;i++){vmov[size(vmov)+1]=var(i+nv-moverdim);}
  if(size(GG)==0){GG=grobcov(F);}
  //"T_GG=grobcov(F)="; GG;
  int j; int k; int te;
  def B0=GG[1][2];
  def H0=GG[1][3][1][1];
  list nGP;
  if (equalideals(B0,ideal(1)) )
  {
    if(version==2){return(locus2(GG,F,DD));}
    else{return(locus0(GG,DD));}
  }
  else
  {
    n=nvars(R);
    ideal vB;
    ideal N;
    for(i=1;i<=size(B0);i++)
    {
      if(subset(variables(B0[i]),vmov)){N[size(N)+1]=B0[i];}
    }
    attrib(N,"IsSB",1);
    N=std(N);
   if((size(N))>=2)
    {
       //def dN=dim(N);
       te=indepparameters(N);
       if(te)
       {
         string("locus detected that the mover must avoid points (",N,") in order to obtain the correct locus");" ";
         //eliminates segments of GG where N is contained in the basis
         nGP[1]=GG[1];
         nGP[1][1]=ideal(1);
         nGP[1][2]=ideal(1);
         def GP=GG;
         ideal BP;
         ideal fBP;
         for(j=2;j<=size(GP);j++)
         {
           te=1; k=1;
           BP=GP[j][2];
          // eliminating multiple factors in the polynomials of BP
           for(mm=1;mm<=size(BP);mm++)
           {
             fBP=factorize(BP[mm],1);
             BP[mm]=1;
             for(m=1;m<=size(fBP);m++)
             {
               BP[mm]=BP[mm]*fBP[m];
             }
           }
           // end eliminating multiple factors
           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
    {
      nGP=GG;
      " ";string("Unavoidable ",moverdim,"-dimensional locus");
      //string("Try option 'vmov',ideal(of mover variables) to avoid some point of the mover");
      //" ";"Elements of the basis of the generic segment in mover variables=";  N;" ";
      list L; return(L);
    }
  }
  if(comment>0){"Input for locus2 GB="; nGP; "input for locus  F="; F;}
  if(version==2)
  {
    //"T_nGP="; nGP;
    def LL=locus2(nGP,F,DD);
  }
  else{def LL=locus0(nGP,DD);}
  // kill @RP; kill @P; kill @R;
  return(LL);
}
example
{
"EXAMPLE:"; echo = 2;
if(defined(R)){kill R;}
ring R=(0,x,y),(x1,y1),dp;
short=0;

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

locus(S96);

// EXAMPLE
// Consider two parallel planes z1=-1 and z1=1, and two orthogonal parabolas on them.
// Determine the locus generated by the lines that rely the two parabolas
// through the points having parallel tangent vectors.

if(defined(R)){kill R;}
ring R=(0,x,y,z),(lam,x2,y2,z2,x1,y1,z1), lp;
short=0;

ideal L=z1+1,
            x1^2-y1,
            z2-1,
            y2^2-x2,
            4*x1*y2-1,
            x-x1-lam*(x2-x1),
            y-y1-lam*(y2-y1),
            z-z1-lam*(z2-z1);

locus(L);  // uses "moverdim",7

// Now observe the effect on the antiimage of option moverdim:
// It does not take into account that lam is free.
locus(L,"moverdim",3);

// Observe that with the option "mpverdim",3 the taxonomy becomes Special because considering only x1,y1,z1 as mover variables does not take into account that lam is free in the global anti-image.
}

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

def L96=locus(S96);
L96;

locusdg(L96);
}

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

locusto(locusdg(L));
}

// envelop
// Input:
//   poly F: the polynomial defining the family of hypersurfaces in ring R=0,(x_1,..,x_n),(u_1,..,u_m),lp;
//   ideal C=g1,..,g_{n-1}:  the set of constraints;
//   options.
// Output: the components of the envolvent;
proc envelop(poly F, ideal C, list #)
"USAGE: envelop(poly F,ideal C[,options]);
       poly F must represent the family of hyper-surfaces for
       which on want to compute its envelop. ideal C must be
       the ideal of restrictions on the variables defining the
       family, and should contain less polynomials than the
       number of variables. (x_1,..,x_n) are the variables of
       the hyper-surfaces of F, that are considered as
       parameters of the parametric ring. (u_1,..,u_m) are
       the parameteres of the hyper-surfaces, that are
       considered as variables of the parametric ring.
       Calling sequence:
       ring R=(0,x_1,..,x_n),(u_1,..,u_m),lp;
       poly F=F(x_1,..,x_n,u_1,..,u_m);
       ideal C=g_1(u_1,..u_m),..,g_s(u_1,..u_m);
       envelop(F,C[,options]);   where s<m.
RETURN: The output is a list of the components [C_1, .. , C_n]
       of the locus. Each component is given by
       Ci=[pi,[pi1,..pi_s_i],tax] where
       pi,[pi1,..pi_s_i] is the canonical P-representation of
       the component.
       Concerning tax: (see help for locus)
       For normal-point components is
         tax=[d,taxonomy,anti-image], being
         d=dimension of the anti-image
         taxonomy=\"Normal\" or \"Special\"
         anti-image=values of the mover corresponding
         to the component
       For non-normal-point components is
         tax=[d,taxonomy]
         d=dimension of the component
         taxonomy=\"Accumulation\" or \"Degenerate\".
OPTIONS: An option is a pair of arguments: string, integer.
       To modify the default options,
       pairs of arguments -option name, value- of valid options
       must be added to the call. The algorithm allows the
       following option as pair of arguments:
       \"comments\", c: by default it is 0, but it can be set to 1.
       \"anti-image\", a: by default a=1 and the anti-image is
       shown also for \"Normal\" components.
       For a=0, it is not shown.
NOTE: grobcov and locus are called internally.
       The basering R, must be of the form Q[x][u]
       (x=parameters, u=variables).
       This routine uses the generalized definition of envelop
       introduced in the book
       A. Montes. \"The Groebner Cover\" (Discussing Parametric
       Polynomial Systems) not yet published.
KEYWORDS: geometrical locus; locus; envelop
EXAMPLE:  envelop; shows an example"
{
  def RR=basering;
  int tes=0; int i;   int j;  int k; int m;
  int d;
  int dp;
  ideal BBB;
  // new
  ideal CC=C;
  CC[size(CC)+1]=F;
  int movreal=size(variables(CC));
  //Options
  list DD=#;

  list Rx=ringlist(RR);
  ideal vmov;
  int moverdim=movreal;         //size(Rx[1][2]);
  int nv=nvars(RR);
  if(moverdim>nv){moverdim=nv;}
  for(i=1;i<=moverdim;i++){vmov[size(vmov)+1]=var(i+nv-moverdim);}
  //for(i=1;i<=nv;i++)  {vmov[size(vmov)+1]=var(i);}
  int version=2;
  int comment=0;
  ideal Fm;
  int tax=1;

  int familyinfo;
  for(i=1;i<=(size(DD) div 2);i++)
  {
    if(DD[2*i-1]=="anti-image"){tax=DD[2*i];}
    if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    if(DD[2*i-1]=="version"){version=2;DD[2*i]=2;}
    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
    if(DD[2*i-1]=="familyinfo"){familyinfo=0;DD[2*i]=0;}
  };
  int ng=size(C);
  ideal S=F;
  //"T_C="; C;
  for(i=1;i<=size(C);i++){S[size(S)+1]=C[i];}
  //"T_S="; S;
  int s=nv-ng;
  if(s>0)
  {
    matrix M[ng+1][ng+1];
    def cc=comb(nv,ng+1);
    //string("nv=",nv," ng=",ng," s=",s);
    //"T_cc="; cc;
    poly J;
    for(k=1;k<=size(cc);k++)
    {
      for(j=1;j<=ng+1;j++)
      {
        M[1,j]=diff(F,var(cc[k][j]));
      }
      for(i=1;i<=ng;i++)
      {
        for(j=1;j<=ng+1;j++)
        {
          M[i+1,j]=diff(C[i],var(cc[k][j]));
        }
      }
      J=det(M);
      S[size(S)+1]=J;
    }
 }
 if(comment>0){"System S before grobcov ="; S;}
  //def G=grobcov(S,DD);
  list HHH;
  DD[size(DD)+1]="moverdim";   // new option
  DD[size(DD)+1]=moverdim;
  // string("moverdim=",moverdim);
  def L=locus(S,DD);
  // "L="; L;
  list GL;
  ideal fam; ideal env;

  def P=ring(Rx[1]);
  Rx[1]=0;
  def D=ring(Rx);
  def RP=P+D;
  list LL;
  list NormalComp;
  ideal Gi;
  ideal BBBB;
  poly B0;
  return(L);
}
example
{
  "EXAMPLE:"; echo = 2;
if(defined(R)){kill R;}
ring R=(0,x,y),(r,s,y1,x1),lp;
poly F=(x-x1)^2+(y-y1)^2-r;
ideal g=(x1-2*(s+r))^2+(y1-s)^2-s;

def E=envelop(F,g);
E;

E[1][1];

string("Taxonomy=", E[1][3][2]);


// EXAMPLE
// Steiner Deltoid
// 1. Consider the circle x1^2+y1^2-1=0, and a mover point M(x1,y1) on it.
// 2. Consider the triangle A(0,1), B(-1,0), C(1,0).
// 3. Consider lines passing through M perpendicular to two sides of ABC triangle.
// 4. Obtain the envelop of the lines above.
if(defined(R)){kill R;}
ring R=(0,x,y),(x1,y1,x2,y2),lp;
short=0;
ideal C=(x1)^2+(y1)^2-1,
             x2+y2-1,
             x2-y2-x1+y1;
matrix M[3][3]=x,y,1,x2,y2,1,x1,0,1;
poly F=det(M);
// Curves Family F
F;
// Conditions C=
C;
envelop(F,C);
}

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

if(defined(R)){kill R;}
ring R=(0,x,y),(r,s,y1,x1),lp;

poly F=(x-x1)^2+(y-y1)^2-r;
ideal g=(x1-2*(s+r))^2+(y1-s)^2-s;

def E=envelop(F,g);
E;

def A=AssocTanToEnv(F,g,E[1][1][1]);
A;

def M1=coef(A[2][1],x1);
def M2=coef(A[2][2],y1);
def M3=coef(A[2][3],s);
def M4=coef(A[2][4],r);

"x1=";-M1[2,2]/M1[2,1];

"y1=";-M2[2,2]/M2[2,1];

"s=";-M3[2,2]/M3[2,1];

"r=";-M4[2,2]/M4[2,1];
}

proc FamElemsAtEnvCompPoints(poly F,ideal C, ideal E)
"USAGE: FamElemsAtEnvCompPoints(poly F,ideal C,ideal E);
       poly F must be the family of hyper-surfaces whose
       envelope is analyzed. It must be defined in the ring
       R=Q[x_1.,,x_n][u_1,..,u_m],
       ideal C must be the ideal of restrictions on the
       variables u1,..um.
       Must contain less polynomials than m.
       ideal E must be a component of
       envelop(F,C), previously computed.
       After computing the envelop of a family of
       hyper-surfaces F, with constraints C,
       Consider a component with top E. The call to
       FamElemsAtEnvCompPoints(F,C,E)
       returns the parameter values of the
       set of all hyper-surfaces of the family passing at
       one point of the envelop component E.
       Calling sequence:
       ring R=(0,x_1,..,x_n),(u_1,..,u_m),lp;
       poly F=F(x_1,..,x_n,u_1,..,u_m);
       ideal C=g_1(u_1,..u_m),..,g_s(u_1,..u_m);
       poly E(x_1,..,x_n);
       FamElemsAtEnvCompPoints(F,C,E[,options]);
RETURN: list [lpp,basis,segment]. The basis determines
       the parameter values of the of hyper-surfaces that
       pass at a fixed point of the envelop component E.
       The lpp determines the dimension of the set.
       The segment is the component and is given in Prep.
       Fixing the values of (x_1,..,x_n) inside E, the basis
       allows to determine the values of the parameters
       (u_1,..u_m), of the hyper-surfaces passing at a point
       of E. See the book
       A. Montes. \"The Groebner Cover\" (Discussing
       Parametric Polynomial Systems).
NOTE: grobcov is called internally.
       The basering R, must be of the form Q[a][x]
       (a=parameters, x=variables).
KEYWORDS: geometrical locus; locus; envelop; associated tangent
EXAMPLE:  FamElemsAtEnvCompPoints; shows an example"
{
  int i;
  int moreinfo=0;
  ideal S=C;
  ideal EE=E;
  S[size(S)+1]=F;
  //S[size(S)+1]=E;
  for(i=1;i<=size(E);i++){S[size(S)+1]=E[i];}
  def G=grobcov(S);
  list GG;
  for(i=2; i<=size(G); i++)
  {
    GG[size(GG)+1]=G[i];
  }


  if(moreinfo>0){return(GG);}
  else
  {
    int t=0;
    list HH;
    i=1;
    while(t==0 and i<=size(G))
    {
      //string("T_G[",i,"][3][1][1][1]="); G[i][3][1][1][1];
      //string("T_EE="); EE;
      if(G[i][3][1][1][1]==EE)
      {
         t=1;
         HH=G[i];
      }
      i++;
    }
    return(HH);
  }
}
example
{
  "EXAMPLE:"; echo = 2;

if(defined(R)){kill R;}
ring R=(0,x,y),(t),dp;
short=0;
poly F=(x-5*t)^2+y^2-9*t^2;
ideal C;

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

// E is a component of the envelope:
ideal E=Env[1][1];

E;

def A=AssocTanToEnv(F,C,E);
A;

// The basis of the parameter values of the associated
//    tangent component is
A[2][1];

// Thus t=-(5/12)*y, and  the associated tangent family
//    element at (x,y) is

subst(F,t,-(5/12)*y);

def FE=FamElemsAtEnvCompPoints(F,C,E);
FE;

factorize(FE[2][1]);

// Thus the unique family element passing through the envelope point (x,y)
// corresponds to the value of t of the Associated Tangent

// EXAMPLE:
if(defined(R)){kill R;}
ring R=(0,x,y),(r,s,y1,x1),lp;

poly F=(x-x1)^2+(y-y1)^2-r;
ideal g=(x1-2*(s+r))^2+(y1-s)^2-s;

def E=envelop(F,g);
E;

def A=AssocTanToEnv(F,g,E[1][1][1]);
A;

def M1=coef(A[2][1],x1);
def M2=coef(A[2][2],y1);
def M3=coef(A[2][3],s);
def M4=coef(A[2][4],r);

// The parameter values corresponding to the family
//    element tangent at point (x,y) of the envelope are:
"x1=";-M1[2,2]/M1[2,1];

"y1=";-M2[2,2]/M2[2,1];

"s=";-M3[2,2]/M3[2,1];

"r=";-M4[2,2]/M4[2,1];

// Now detect if there are other family elements passing at this point:
def FE=FamElemsAtEnvCompPoints(F,g,E[1][1][1]);
FE;

// FE[1] is the set of lpp. It has dimension 4-2=2.
//     Thus there are points of the envelope at which
//     they pass infinitely many circles of the family.
//     To separe the points of the envelope further analysis must be done.
}

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

discrim(f,x);
}

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

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

// used in NS
// returns 0 if E does not reduce modulo N
// returns 1 if it reduces
static proc redPbasis(ideal E, ideal N)
{
  int i;
  def RR=basering;
  def Rx=ringlist(RR);
  def Lx=Rx;
  def P=ring(Rx[1]);
  setring P;
  def EP=imap(RR,E);
  def NP=imap(RR,N);
  NP=std(NP);
  list L;
  int red=1;
  i=1;
  while(red and (i<=size(EP)))
  {
    if(reduce(EP[i],NP,5)!=0){red=0;}
    i++;
  }
  setring RR;
  return(red);
}


//******************** End locus and envelop ******************************

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

// input ideal F in @R
//          ideal a in @R but only depending on parameters
//          F is a generating ideal in V(a);
// output:  ideal b in @R but depending only on parameters
//              ideal G=GBasis(F) in V(a) \ V(b)
proc WLemma(ideal F,ideal a, list #)
"USAGE: WLemma(F,A[,options]);
       The first argument ideal F in Q[x_1,..,x_n][u_1,..,u_m];
       The second argument ideal A in Q[x_1,..,x_n].
       Calling sequence:
       ring R=(0,x_1,..,x_n),(u_1,..,u_m),lp;
       ideal  F=f_1(x_1,..,x_n,u_1,..,u_m),..,
           f_s(x_1,..,x_n,u_1,..,u_m);
       ideal A=g_1(u_1,..u_m),..,g_s(u_1,..u_m);
       list # : Options
       Calling sequence:
       WLemma(F,A[,options]);

       Given the ideal F  and ideal A
       it returns the list (lpp,B,S)  were B is the
       reduced Groebner basis of the specialized F over
       the segment S, subset of V(A) with top A,
       determined by Wibmer's Lemma.
       S is determined in P-representation
       (or optionally in C-representation). The basis is
       given by I-regular functions.
OPTIONS: either (\"rep\", 0) or (\"rep\",1) the representation of
       the resulting segment, by default is
       0 =P-representation, (default) but can be set to
       1=C-representation.
RETURN: list of [lpp,B,S] =
       [leading power product, basis,segment],
       being B the reduced Groebner Basis given by
       I-regular functions in full representation, of
       the specialized ideal F on the segment S,
       subset of V(A) with top A.
       given in P- or C-representation.
       It is the result of Wibmer's Lemma. See
       A. Montes , M. Wibmer, \"Groebner Bases for
       Polynomial Systems with parameters\".
       JSC 45 (2010) 1391-1425.)
       or the book
       A. Montes. \"The Groebner Cover\" (Discussing
       Parametric Polynomial Systems).

NOTE: The basering R, must be of the form Q[a][x]
      (a=parameters, x=variables).
KEYWORDS: Wibmer's Lemma
EXAMPLE:  WLemma; shows an example"
{
  list L=#;
  int rep=0;
  int i; int j;
  if(size(L)>0)
  {
    for(i=1;i<=size(L) div 2;i++)
    {
      if(L[2*i-1]=="rep"){rep=L[2*i];}
    }
  }
  def RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  Rx[1]=0;
  def D=ring(Rx);
  def RP=D+P;
  setring(RP);
  ideal FF=imap(RR,F);
  FF=std(FF);
  ideal AA=imap(RR,a);
  AA=std(AA);
  FF=FF,AA;
  FF=std(FF);
  ideal FFa;
  poly r;
  for(i=1; i<=size(FF);i++)
  {
    r=reduce(FF[i],AA);
    if(r!=0){FFa[size(FFa)+1]=r;}
  }
  // FFa is GB(F+a,>xa)
  setring RR;
  ideal Fa=imap(RP,FFa);
  ideal AAA=imap(RP,AA);
  ideal lppFa;
  ideal lcFa;
  for(i=1;i<=size(Fa);i++)
  {
    lppFa[size(lppFa)+1]=leadmonom(Fa[i]);
    lcFa[size(lcFa)+1]=leadcoef(Fa[i]);
  }
  // "T_lppFa="; lppFa;
  // "T_lcFa="; lcFa;
  setring RP;
  ideal lccr=imap(RR,lppFa);
  lccr=std(lccr);
  setring RR;
  ideal lcc=imap(RP,lccr);
  list J; list Jx;
  ideal Jci;
  ideal Jxi;
  list B;
  // "T_lcc="; lcc;
  for(i=1;i<=size(lcc);i++)
  {
    kill Jci; ideal Jci; kill Jxi; ideal Jxi;
    for(j=1;j<=size(Fa);j++)
    {
      if(lppFa[j]==lcc[i])
      {
        Jci[size(Jci)+1]=lcFa[j];
        Jxi[size(Jxi)+1]=Fa[j];
      }
    }
    J[size(J)+1]=Jci;
    B[size(B)+1]=Jxi;
  }
 // "T_J="; J;
  if(size(J)>0)
  {
    setring P;
    list Jp=imap(RR,J);
    ideal JL=product(Jp);
    // JL=prod(lc(Fa))
    def AAAA=imap(RR,AAA);
    // "T_AAA="; AAA;
    // "T_JLA="; JLA;
    def CPR=Crep(AAAA, JL);
    def PPR=Prep(AAAA,JL);
  }
  setring RR;
  if(size(J)>0)
  {
    def JLA=imap(P,JL);
    def PR=imap(P,PPR);
    def CR=imap(P,CPR);
    // PR=Prep(a,b)
    // CR=Crep(a,b)
    for(i=1;i<=size(B);i++)
    {
      for(j=1;j<=size(B[i]);j++)
      {
        B[i][j]=pnormalf(B[i][j],CR[1],CR[2]);
      }
      B[i]=elimrepeated(B[i]);
    }
     // B=reduced basis on CR
     //"T_PR="; PR;
     //"T_CR="; CR;
     //"T_B="; B;
    if(rep==1){return(list(lcc,B,CR));}
    else{return(list(lcc,B,PR));}
  }
  else
  {
    "PIP";
    lcc=ideal(0);
    B=ideal(0);
    list NN;
    NN[1]=list(AAA,ideal(1));
    return(list(lcc,B,NN));
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  if(defined(RE)){kill RE;}
  ring RE=(0,a,b,c,d,e,f),(x,y),lp;
  ideal F=a*x^2+b*x*y+c*y^2,d*x^2+e*x*y+f*y^2;
  ideal A=a*e-b*d;

WLemma(F,A);

WLemma(F,A,"rep",1);
}

// Detect if ideal J is in the list of ideals L
// Input ideal J,  list L
// Output:  1 if J is in L, and 0 if not
static proc idinlist(ideal J,list L)
{
  int i=0;
  int te=0;
  while(te==0 and i<=size(L)-1)
  {
    i++;
    if(equalideals(J,L[i])){te=1;}
  }
  return(te);
}

// input ideal F in Kˆa][x]
// output:  a disjoint CGS in full representation of the ideal F using Wibmer's Lemma WLemma
proc WLcgs(ideal F)
USAGE: WLcgs(ideal F)
       // WLemma(F,A[,options]);
       ideal F in Q[x_1,..,x_n][u_1,..,u_m];
       ring R=(0,x_1,..,x_n),(u_1,..,u_m),lp;
       ideal  F=f_1(x_1,..,x_n,u_1,..,u_m),..,
           f_s(x_1,..,x_n,u_1,..,u_m);
       list # : Options
       Calling sequence:
       WLcgs(ideal F)

       Given the ideal F
       it returns the list of (lpp,B,S)[i] of the grobcov.
       B[i] is the reduced Groebner basis in full-representation of the specialized F over
       the segments S[i],
       S is determined in P-representation
       (or optionally in C-representation). The basis is
       given by I-regular functions in full-representation
OPTIONS: either (\"rep\", 0) or (\"rep\",1) the representation of
       the resulting segment, by default is
       0 =P-representation, (default) but can be set to
       1=C-representation.
RETURN: list of [lpp,B,S][i] =
       [leading power product, basis,segment],
       being B[i] the reduced Groebner Basis given by
       I-regular functions in full representation, of
       the specialized ideal F on the segment S[i],
       given in P-representation.
       It is the result of Wibmer's Lemma. See
       A. Montes , M. Wibmer, \"Groebner Bases for
       Polynomial Systems with parameters\".
       JSC 45 (2010) 1391-1425.)
       or the book
       A. Montes. \"The Groebner Cover\" (Discussing
       Parametric Polynomial Systems).

NOTE: The basering R, must be of the form Q[a][x]
      (a=parameters, x=variables).
KEYWORDS: Wibmer's Lemma
EXAMPLE:  WLcgs; shows an example"
{
  int i,j;
  list Etot;
  Etot[1]=ideal(0);
  list Epend=Etot;
  list G;
  list G0;
  list N;
  ideal a;
  while (size(Epend)>0)
  {
    a=Epend[1];
    Epend=elimidealfromlist(Epend,a);
    G0=WLemma(F,a);
    if(size(G0)>0)
    {
      G[size(G)+1]=G0;
      //"T_G0="; G0;
      N=G0[3][1][2];
      //"T_N="; N;
      for(i=1;i<=size(N);i++)
      {
        if(not(equalideals(N[i],ideal(1)) or idinlist(N[i],Etot)))
        {
          Etot[size(Etot)+1]=N[i];
          Epend[size(Epend)+1]=N[i];
        }
        //"T_i="; i;
      }
      //"T_Etot="; Etot;
      //"T_Epend=";Epend;
    }
  }
  return(G);
}
example
{
  "EXAMPLE:"; echo = 2;
  if(defined(RRR)){kill RRR;}
  ring RRR=(0,b,c,d,e,f),(x,y,t),lp;
  short=0;
  ideal S=x^2+2*c*x*y+2*d*x*t+b*y^2+2*e*y*t+f*t^2,
            x+c*y+d*t,c*x+b*y+e*t;
  grobcov(S);
  WLcgs(S);
}


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


// Not used
static proc redbasis(ideal B, ideal C)
{
  int i;
  def RR=basering;
  def Rx=ringlist(RR);
  def Lx=Rx;
  def P=ring(Rx[1]);
  Lx[1]=0;
  def D=ring(Lx);
  def RP=D+P;
  setring RP;
  ideal BB=imap(RR,B);
  ideal CC=imap(RR,C);
  attrib(CC,"IsSB",1);
  CC=std(CC);
  for(i=1;i<=size(BB);i++)
  {
    BB[i]=reduce(BB[i],CC);
  }
  setring(RR);
  def BBB=imap(RP,BB);
  return(BBB);
}


// not used
// Input: ideals E, N
// Output: the ideal N without the polynomials in E
//   Works in any kind of ideal
static proc idminusid(ideal E,ideal N)
{
  int i; int j;
  ideal h;
  int te=1;
  for(i=1; i<=size(N);i++)
  {
    te=1;
    for(j=1;j<=size(E);j++)
    {
      if(N[i]==E[j]){te=0;}
    }
    if(te==1){h[size(h)+1]=N[i];}
  }
  return(h);
}

// not used
// eliminar els factors de cada polinomi de F que estiguin a N\ E
static proc simpB(ideal F,ideal E,ideal N)
{
  ideal FF;
  poly ff;
  int i; int j;
  ideal J=idminusid(E,N);
  //"T_J="; J;
  for(i=1;i<=size(F);i++)
  {
    for(j=1;j<=size(J);j++)
    {
      ff=elimfacsinP(F[i],J[j]);
    }
    FF[size(FF)+1]=ff;
  }
  return(FF);
}

// used in simpB that is not used
static proc elimfacsinP(poly f,poly g)
{
  def RR=basering;
  def Rx=ringlist(RR);
  int i; int j;
  int n=size(Rx[1][2]);
  def Lx=Rx;
  Lx[1]=0;
  def D=ring(Lx);
  def P=ring(Rx[1]);
  def RP=D+P;
  setring P;
  ideal vp;
  for(i=1;i<=n;i++)
  {
    vp[size(vp)+1]=var(i);
  }
  setring RP;
  def gg=imap(RR,g);
  ideal vpr=imap(P,vp);
  poly ff=imap(RR,f);
  def L=factorize(ff);
  def L1=L[1];
  poly p=1;
  for(i=1;i<=size(L1);i++)
  {
    if(L1[i]==gg){;}
    else{p=p*L1[i];}
  }
  setring RR;
  def pp=imap(RP,p);
  return(pp);
}

//****************************** Begin ADGT *************************

// used in ADGT
// Given G=grobcov(F,"rep",1) to have the GC in C-representation,
// Grob1Levels determines the canonical levels of the constructible subset
// of the parameter space for which there exist solutions of F
// To be called in Q[a][x]
proc Grob1Levels(list G)
"USAGE: Grob1Levels(list G);
       G is the output of grobcov(F,\"rep\",1)
       for obtaining the segments in C-rep.
       Then Grob!Levels, selects the set of segments S of G having solutions
       (i.e. with basis different from 1), and determines the canonical levels
       of this constructible set.
       To be called in a ring Q[a][x].
RETURN: The list of ideals
       [a1,a2,...,at]
       representing the closures of the canonical levels of S
       and its complement C wrt to the closure of S.

       The levels of S and C are
       Levels of S:  [a1,a2],[a3,a4],...
       Levels of C:  [a2,a3],[a4,a5],...
       S=V(a1) \ V(a2) u V(a3) \ V(a4) u ...
       C=V(a2 \ V(a3) u V(a4) \ V(a5) u ...
       The expression of S can be obtained from the
       output of Grob1Levels by
       the call to Levels.
NOTE: The algorithm was described in
       J.M. Brunat, A. Montes. \"Computing the canonical
       representation of constructible sets.\"
       Math.  Comput. Sci. (2016) 19: 165-178.
KEYWORDS: constructible set; locally closed set; canonical form
EXAMPLE:  Grob1Levels; shows an example"
{
  int i;
  list S;
  def RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  setring(RR);
  for(i=1;i<=size(G);i++)     // select the segments with solutions
  {
    if(not(G[i][1][1][1]==1))
    {
      S[size(S)+1]=G[i][3];
    }
  }
  if(size(S)==0)
  {
    list L;
    L[1]=1;
  }
  else
  {
    setring P;
    def SP=imap(RR,S);
    list LP=ConsLevels(SP);   // construct the levels of the constructible
    setring RR;
    def L=imap(P,LP);
  }
  return(L);
}
example
{
"EXAMPLE:"; echo = 2;
  if (defined(R)) {kill R;}
  ring R=(0,x,y),(x1,y1,x2,y2),lp;
  ideal F=-y*x1+(x-1)*y1+y,
            (x-1)*(x1+1)+y*y1,
            -y*x2+(x+1)*y2-y,
            (x+1)*(x2-1)+y*y2,
            (x1-x)^2+y1^2-(x1-x)^2-y2^2;

def G=grobcov(F,"rep",1);

G;

def L=Grob1Levels(G);

L;

def LL=Levels(L);

LL;
}

// Auxiliary rutine for intersecting ideal only in the parameters a
// when the ideals are defined in Q[a][x]
proc intersectpar(list S)
"USAGE: interectpar(list of ideals S)
    list S=ideal I1,...,ideal Ik;
RETURN: The intersection of the ideals I1 ...  Ik  in Q[x,a]
NOTE: The routine is called in Q[a][x],
    The ideals I1,..,Ik can be ideals depending only on [a] or on [x,a]
EXAMPLE: intersectpar shows an example"
{
  def RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  Rx[1]=0;
  def D=ring(Rx);
  def RP=D+P;
  setring(RP);
  def SP=imap(RR,S);
  //"T_SP="; SP;
  //"T_typeof(SP[1])="; typeof(SP[1]);
  ideal EP;
  EP=SP[1];
  int i;
  for(i=2;i<=size(SP);i++)
  {
    EP=intersect(EP,SP[i]);
  }
  //def EP=intersect(SPL);
  setring RR;
  def E=imap(RP,EP);
  return(E);
}
example
{
"EXAMPLE:"; echo = 2;
if(defined(R)){kill R;}
ring R=(0,x,y,z),(x1,y1,z1),lp;
ideal I1=x+y*z*x1;
ideal I2=x-y*z*y1;
ideal I3=x+y+z*z1;
list S=I1,I2,I3;
S;

intersectpar(S);
}


proc ADGT(ideal H,ideal T,ideal H1,ideal T1,list #)
"USAGE: ADGT(ideal H, ideal T, ideal H1,ideal T1[,options]);
      H: ideal in Q[a][x] hypothesis
      T: ideal in Q[a][x] thesis
      H1: ideal in Q[a][x] negative hypothesis, only dependent on [a]
      T1: ideal in Q[a][x] negative thesis
           of the Proposition (H and not H1) => (T and not T1)
RETURN: The list  [[1,S1],[2,S2],..],
       S1, S2, .. represent the set of parameter values
       that must be verified as supplementary
       conditions for the Proposition to become a Theorem.
       They are given by default in Prep.
       If the proposition is generally true,
       (the proposition is already a theorem), then
       the generic segment of the internal grobcov
       called is also returned to provide information
       about the values of the variables determined
       for every value of the parameters.
       If the propsition is false for every values of the
       parameters, then the emply list is returned.
OPTIONS: An option is a pair of arguments: string,
       integer. To modify the default options, pairs
       of arguments -option name, value- of valid options
       must be added to the call.
       Option \"rep\",0-1: 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,
       Option \"gseg\",0-1: The default is \"gseg\",1
       and then when the proposition is generally true,
       ADGT outputs a second element which is the
       \"generic segment\" to provide supplementary information.
       With option \"gseg\",0 this is avoided.
       Option \"neg\", 0,1:  The default is \"neg\",0
       With option  \"neg\",0  Rabinovitch trick is used for negative hypothesis and thesis
       With option  \"neg\",1  Difference of constructible sets is used instead.
NOTE:  The basering R, must be of the form Q[a][x],
       (a=parameters, x=variables), and
       should be defined previously. The ideals
       must be defined on R.
KEYWORDS: Automatic Deduction; Automatic Demonstration; Geometric Theorems
EXAMPLE:  ADGT shows an example"
{
  int i; int j;
  def RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  Rx[1]=0;
  def D=ring(Rx);
  def RP=D+P;
  setring RR;
  list Lopt=#;
  int start=timer;
  // options
  int rep=0;
  int gseg=1;
  int neg=0;

  for(i=1;i<=size(Lopt) div 2;i++)
  {
    if(Lopt[2*i-1]=="rep"){rep=Lopt[2*i];}
    if(Lopt[2*i-1]=="gseg"){gseg=Lopt[2*i];}
    if(Lopt[2*i-1]=="neg"){neg=Lopt[2*i];}
  }
  // begin proc

  if(neg==0)
  {
  // Default option "neg",0  uses Rabinovitch for negative hyothesis and thesis
  // Option "neg",1  uses Diference of constructive sets for negative hyothesis and thesis
  if(equalideals(T1,ideal(1))) //nonnullT==0)
  {
    def F=H;
    for(i=1;i<=size(T);i++)
    {
      F[size(F)+1]=T[i];
    }
    list G2;
    if(not(equalideals(H1,ideal(1))))  //nonnullH==1)
    {
      G2=grobcov(F,"nonnull",H1,"rep",1);
    }
    else
    {
      G2=grobcov(F,"rep",1);
    }
  }
  else
  {
    def F=H;
    for(i=1;i<=size(T1);i++)
    {
      F[size(F)+1]=T1[i];
    }
    list G3=grobcov(F,"rep",1);
    def L0=Grob1Levels(G3);
    setring P;
    def LP=imap(RR,L0);
    def H1P=imap(RR,H1);
    if(not(equalideals(H1P,ideal(1))))  //nonnullH==1)
    {
      i=1;
      while(i<=size(LP))
      {
        if(i mod 2==1)
        {
          LP[i]=intersect(LP[i],H1P);
        }
        i++;
      }
    }
    setring RR;
    L0=imap(P,LP);
    F=H;
    for(i=1;i<=size(T);i++)
    {
      F[size(F)+1]=T[i];
    }
    def G2=grobcov(F,"nonnull",L0[1],"rep",1);
    list G1;
    int m=size(L0);
    int r=m div 2;
    if(m mod 2==0){r=r-1;}
    //"L0="; L0;
     for(i=1;i<=r;i++)
    {
      G1=grobcov(F,"null",L0[2*i],"nonnull",L0[2*i+1],"rep",1);
      for(j=1;j<=size(G1);j++)
      {
        if(not(equalideals(G1[j][1],ideal(1))))
        {
          G2[size(G2)+1]=G1[j];
        }
      }
    }
  }
  if(size(G2)==0)
  {
    list L;
    // L[1]=1;
  }
  else
  {
    list L=Levels(Grob1Levels(G2));
    if(rep==0 and size(L)>0)
    {
      setring P;
      def LFP=imap(RR,L);
      // list LFP1=LFP;
      for(i=1;i<=size(LFP);i++)
      {
        LFP[i][2]=Prep(LFP[i][2][1],LFP[i][2][2]);
      }
      setring RR;
      L=imap(P,LFP);
    }
    if(equalideals(G2[1][3][1][1],0) and not(equalideals(G2[1][1],ideal(1))))
    {
      list GL=G2[1];
      list LL=GL[3];
      if(rep==0)
      {
        setring P;
        def LLP=imap(RR,LL);
        LLP=Prep(LLP[1],LLP[2]);
        setring RR;
        LL=imap(P,LLP);
        GL[3]=LL;
      }
      if(gseg==1)
      {
        L[size(L)+1]=list("Generic segment",GL);
      }
    }
  }
  return(L);
  }
  else
  {
    if (neg==1)
    {
      def LL=ADGTDif(H,T,H1,T1,#);
      return(LL);
    }
  }
}
example
{
"EXAMPLE:"; echo = 2;

// Determine the supplementary conditions
// for the non-degenerate  triangle A(-1,0), B(1,0), C(x,y)
// to have an ortic non-degenerate isosceles triangle

if(defined(R)){kill R;}
ring R=(0,x,y),(x2,y2,x1,y1),lp;

// Hypothesis H: the triangle  A1(x1,y1), B1(x2,y2), C1(x,0), is the
// orthic triangle of ABC

ideal H=-y*x1+(x-1)*y1+y,
         (x-1)*(x1+1)+y*y1,
         -y*x2+(x+1)*y2-y,
         (x+1)*(x2-1)+y*y2;

// Thesis T: the orthic triangle is isosceles
ideal T=(x1-x)^2+y1^2-(x2-x)^2-y2^2;

// Negative Hypothesis H1: ABC is non-degenerate
ideal H1=y;

 // Negative Thesis T1: the orthic triangle is non-degenerate
ideal T1=x*(y1-y2)-y*(x1-x2)+x1*y2-x2*y1;

// Complementary conditions for the
// Proposition (H and not H1) => (T and not T1)
// to be true

ADGT(H,T,H1,T1);

// Now using diference of constructible sets for negative hypothesis and thesis

ADGT(H,T,H1,T1,"neg",1);

// The results are identical using both methods for the negative propositions
//  - Rabinovitch or
//  - DifConsLCSets

// EXAMPLE

// Automatic Theorem Proving.
// The nine points circle theorem.

// Vertices of the triangle: A(-2,0), B(2,0), C(2a,2b)

// Height foot: A1(x1,y1),
// Height foot: B1(x2,y2),
// Height foot: C1(2a,0)

// Middle point BC:  A2(a+1,b)
// Middle point CA:  B2 (a-1,b)
// Middle point AB:  C2(0,0)

// Ortocenter:   O(2x0,2y0)
// Middle point of A and O:  A3(x0-1,y0)
// Middle point of B and O:  B3(x0+1,y0)
// Middle point of C and O:  C3(x0+a,y0+b)

// Nine points circle: c:=(X-x3)^2+(Y-y3)^2-r2

if (defined(R1)){kill R1;}

ring R1=(0,a,b),(x1,y1,x2,y2,x0,y0,x3,y3,r2),dp;

short=0;

ideal H=-x1*b+(a-1)*y1+2*b,
      (a-1)*x1+b*y1+2*a-2,
      -x2*b+(a+1)*y2-2*b,
      (a+1)*x2+b*y2-2*a-2,
      -x0*y1+(x1+2)*y0-y1,
      -x0*y2+(x2-2)*y0+y2;

ideal T=(x1-2*x3)^2+(y1-2*y3)^2-r2,
       (a+1-2*x3)^2+(b-2*y3)^2-r2,
       (x0-1-2*x3)^2+(y0-2*y3)^2-r2,
       (x2-2*x3)^2+(y2-2*y3)^2-r2,
       (a-1-2*x3)^2+(b-2*y3)^2-r2,
       (x0+1-2*x3)^2+(y0-2*y3)^2-r2,
       (2*a-2*x3)^2+4*y3^2-r2,
       4*x3^2+4*y3^2-r2,
       (x0+a-2*x3)^2+(y0+b-2*y3)^2-r2;

 ADGT(H,T,b,1);

// Thus the nine point circle theorem is true for all real points excepting b=0.
}

static proc ADGTDif(ideal H,ideal T,ideal H1,ideal T1,list #)
{
  int i; int j;
  def RR=basering;
  def Rx=ringlist(RR);
  def P=ring(Rx[1]);
  Rx[1]=0;
  def D=ring(Rx);
  def RP=D+P;
  setring RR;
  list Lopt=#;
  int start=timer;
  // options
  int rep=0;
  int gseg=1;
  list LLL;
  for(i=1;i<=size(Lopt) div 2;i++)
  {
    if(Lopt[2*i-1]=="rep"){rep=Lopt[2*i];}
    else{if(Lopt[2*i-1]=="gseg"){gseg=Lopt[2*i];}}
  }
  // begin proc
    def F=H;
    for(i=1;i<=size(T);i++)
    {
      F[size(F)+1]=T[i];
    }
    list G2=grobcov(F,"rep",1);
    def L2=Grob1Levels(G2);
    //"L2="; L2;
    ideal FN=T1;
    if(not(equalideals(H1,ideal(1))))
    {
      list FNH1=FN,H1;
      FN=intersectpar(FNH1);
    }
    if(not(equalideals(FN,ideal(1))))
    {
      for(i=1;i<=size(FN);i++)
      {
        F[size(F)+1]=FN[i];
      }
      list G3=grobcov(F,"rep",1);
      def L3=Grob1Levels(G3);
      //"L3="; L3;
    }
    def LL=DifConsLCSets(L2,L3);
    //"LL="; LL;
    LL=ConsLevels(LL);
    LL=Levels(LL);
     //"LL="; LL;
    if(rep==1){return(LL);}
    else
    {
      for(i=1;i<=size(LL);i++)
      {
         LL[i][2]=Prep(LL[i][2][1],LL[i][2][2]);
      }
    }
    return(LL);
}

//****************************** End ADGT *************************
;