Changeset 7023ce in git

Timestamp:

Feb 4, 2015, 10:08:19 AM (9 years ago)

Author:

Stephan Oberfranz <oberfran@…>

Branches:

(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')

Children:

dbf60932968d7dd0d5f400c4209b39e40b442f58

Parents:

ca38640ca6928ad815ae7b2bd1f3516f383315af48d66afcfd578627347de972483894aa987b9f60

Message:

---

Merge branch 'spielwiese' of github.com:Singular/Sources into spielwiese

Files:

• .gdbinit (modified) (1 diff)
• Singular/LIB/JMSConst.lib (modified) (1 diff)
• Singular/LIB/algemodstd.lib (modified) (2 diffs)
• Singular/LIB/bimodules.lib (modified) (1 diff)
• Singular/LIB/dmodvar.lib (modified) (1 diff)
• Singular/LIB/ffsolve.lib (modified) (1 diff)
• Singular/LIB/gradedModules.lib (modified) (25 diffs)
• Singular/LIB/grobcov.lib (modified) (127 diffs)
• Singular/LIB/grwalk.lib (modified) (9 diffs)
• Singular/LIB/homolog.lib (modified) (1 diff)
• Singular/LIB/modwalk.lib (modified) (22 diffs, 1 prop)
• Singular/LIB/noether.lib (modified) (2 diffs)
• Singular/LIB/normal.lib (modified) (11 diffs)
• Singular/LIB/presolve.lib (modified) (6 diffs)
• Singular/LIB/primdec.lib (modified) (122 diffs)
• Singular/LIB/resolve.lib (modified) (7 diffs)
• Singular/LIB/ring.lib (modified) (1 diff)
• Singular/LIB/rinvar.lib (modified) (4 diffs)
• Singular/LIB/rwalk.lib (modified) (11 diffs, 1 prop)
• Singular/LIB/schreyer.lib (modified) (1 diff)
• Singular/LIB/swalk.lib (modified) (1 prop)
• Singular/Makefile.am (modified) (4 diffs)
• Singular/blackbox.cc (modified) (1 diff)
• Singular/dyn_modules/polymake/Makefile.am (modified) (2 diffs)
• Singular/extra.cc (modified) (4 diffs)
• Singular/iparith.cc (modified) (1 diff)
• Singular/ipshell.cc (modified) (6 diffs)
• Singular/subexpr.cc (modified) (1 diff)
• Singular/table.h (modified) (7 diffs)
• Singular/tok.h (modified) (1 diff)
• Singular/walk.cc (modified) (142 diffs, 1 prop)
• Singular/walk.h (modified) (1 diff, 1 prop)
• Tst/Manual/findvars.res.gz.uu (modified) (1 diff)
• Tst/Manual/findvars.tst (modified) (1 diff)
• Tst/Short.lst (modified) (1 diff)
• Tst/Short/bug_529.res.gz.uu (modified) (1 diff)
• Tst/Short/bug_529.stat (modified) (1 diff)
• Tst/Short/bug_tr488.res.gz.uu (added)
• Tst/Short/bug_tr488.stat (added)
• Tst/Short/bug_tr488.tst (added)
• Tst/Short/bug_tr690.res.gz.uu (added)
• Tst/Short/bug_tr690.stat (added)
• Tst/Short/bug_tr690.tst (added)
• Tst/Short/ok_s.lst (modified) (4 diffs)
• Tst/Short/test_unit_ideal_decomposition.res.gz.uu (added)
• Tst/Short/test_unit_ideal_decomposition.stat (added)
• Tst/Short/test_unit_ideal_decomposition.tst (added)
• Tst/Short/test_zero_ideal_decomposition.res.gz.uu (added)
• Tst/Short/test_zero_ideal_decomposition.tst (added)
• doc/NEWS.texi (modified) (2 diffs)
• doc/general.doc (modified) (previous)
• doc/reference.doc (modified) (previous)
• dox/Doxyfile.in (modified) (1 diff)
• factory/configure.ac (modified) (2 diffs)
• kernel/GBEngine/kstd1.cc (modified) (3 diffs)
• kernel/GBEngine/kutil.cc (modified) (9 diffs)
• kernel/GBEngine/syz.cc (modified) (1 diff)
• kernel/GBEngine/syz0.cc (modified) (5 diffs)
• kernel/GBEngine/syz3.cc (modified) (2 diffs)
• kernel/ideals.cc (modified) (5 diffs)
• libpolys/coeffs/OPAE.cc (modified) (2 diffs)
• libpolys/coeffs/OPAEQ.cc (modified) (2 diffs)
• libpolys/coeffs/OPAEp.cc (modified) (2 diffs)
• libpolys/polys/clapsing.cc (modified) (1 diff)
• libpolys/polys/simpleideals.cc (modified) (2 diffs)
• libpolys/polys/simpleideals.h (modified) (1 diff)
• m4/flags.m4 (modified) (1 diff)
• m4/gfanlib-check.m4 (modified) (1 diff)

.gdbinit

@@ -1 +1 @@
 break dErrorBreak
-break omReportError
-break dPolyReportError
 
 

Singular/LIB/JMSConst.lib

@@ -947 +947 @@
 poly h=Minimus(variables(A.h));
 //print(h);
-int l=findvars(h,1)[2][1];
+int l=findvars(h)[2][1];
 if(l!=nvars(basering))
  {

Singular/LIB/algemodstd.lib

@@ -293 +293 @@
 ////////////////////////////////////////////////////////////////////////////////
 
+static proc check_leadmonom_and_size(list L)
+{
+    /*
+     * compare the size of ideals in the list and 
+     * check the corresponding leading monomials
+     * size(L)>=2
+     */
+    ideal J=L[1];
+    int i=size(L);
+    int sc=ncols(J);
+    int j,k;
+    poly g=leadmonom(J[1]);
+    for(j=1;j<=i;j++)
+    {
+        if(ncols(L[j])!=sc)
+        {
+            return(0);
+        }
+    }
+    for(k=2;k<=i;k++)
+    {
+        for(j=1;j<=sc;j++)
+        {
+            if(leadmonom(J[j])!=leadmonom(L[k][j]))
+            {
+                return(0);
+            }
+        }
+    }
+    return(1);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
 static proc LiftPolyCRT(ideal I)
 {
@@ -324 +358 @@
         }
     }
-    // apply CRT for polynomials
-    II =chinrempoly(Lk,LL),f;
+    if(check_leadmonom_and_size(Lk))
+    {
+        // apply CRT for polynomials
+        II =chinrempoly(Lk,LL),f;
+    }
+    else
+    {
+        II=0;
+    }
     return(II);
 }

Singular/LIB/bimodules.lib

@@ -722 +722 @@
   ideal J;
   // determine whether I is a left Groebner basis
-  if (attrib(I,"isSB") == 1)
+  if (attrib(I,"isSB"))
   {
     J = I;

Singular/LIB/dmodvar.lib

@@ -698 +698 @@
   int n = nvars(basering);
   int c;
-  if (attrib(I,"isSB") == 1)
+  if (attrib(I,"isSB"))
   {
     c = n - dim(I);

Singular/LIB/ffsolve.lib

@@ -383 +383 @@
         }
       }
-      varinfo = findvars(univar_part,1);
+      varinfo = varaibles(univar_part);
       unsolved_vars = varinfo[3];
       unsolved_var_nums = varinfo[4];

Singular/LIB/gradedModules.lib

@@ -32 +32 @@
 ";
 
+
+
+
 //////////////////////////////////////////////////////////////////////////////////////////////////////////
 // . view graded module/map
@@ -140 +143 @@
 static proc issorted( intvec g, int s )
 {
+  g = s * g; //  "g: ", g;
   int i = size(g);
 
   for(; i > 1; i--)
   {
-    if( s*(g[i] - g[i-1]) < 0 )
+    if( (g[i] - g[i-1]) < 0 )
     {
       return (0);
@@ -158 +162 @@
 "
 {
+  gr = s * gr;
   int m = size(gr);
-
   intvec pivot; 
 
@@ -176 +180 @@
     for(Bn=1; Bn <= (m-Bi); Bn++)
     {
-      D = s*(gr[pivot[Bn]] - gr[pivot[Bn+1]]); // sort gr 
+      D = (gr[pivot[Bn]] - gr[pivot[Bn+1]]); // sort gr 
       P = (D > 0) or ( (D == 0) && ((pivot[Bn]-pivot[Bn+1]) > 0) ); // stability!?
       if(P)
@@ -189 +193 @@
   }
 
-  ASSUME(1, issorted(gr, s));
 
   if( flag && 0 ) // output details in case of non-identical permutation?
   {
     // grades & ordering permutation : gr[pivot] should be sorted:
-    "";
+    "s: ", s;
+    "gr: ", gr;
     "pivot: ",    pivot;
-    "";
-  }
+  }
+
+  ASSUME(1, issorted(intvec(gr[pivot]), 1));
 
   return (pivot);
 }
 
+// Q@Wolfram: what should be done with zero gens!?
 proc grdeg( def M, list # )
 "graded degrees of gens(i.e. columns)"
-
 {
  // use initial grading if given
@@ -219 +224 @@
     kill w;
   }
-
-  // try to homogenize otherwise 
-  if( typeof(attrib( M, "isHomog")) != "intvec" )
-  { // no such attribute? //   ERROR("No grading!");
-    ASSUME(0, /* input must be graded: cannot compute homogenizing grading */ homog(M) ); 
-  }
-
-  
-  ASSUME(0, /* input must be graded! */ typeof(attrib(M, "isHomog")) == "intvec" );
-  // input should be graded:
+  
+  ASSUME(1, grtest(M) );
+
   def w = attrib(M, "isHomog"); // grading weights?
-  ASSUME(0, /* input must be correctly graded! */ nrows(M) == size(w) );   
-
-  if( size(M) == 0 )  {    return (w);  } // TODO???
+
+  if( size(M) == 0 ){ return (w); } // TODO: Q@Wolfram!???
 
   int m = ncols(M); // m > 0 in Singular!
@@ -249 +246 @@
 }
 
-
-static proc order( def M, int s, list #)
-"helper for reorder: compute graded degrees and the permutation to sort them"
-{
-  if( size(#) > 0 ) { intvec gr = grdeg( M, #[1] ); }
-  else              { intvec gr = grdeg( M       ); }
-  
+static proc reorder(def M, int s)
+"
+Reorder gens of M: compute graded degrees and the permutation to sort them
+"
+{
+  // input should be graded:
+  ASSUME(1, grtest(M) );
+
+  intvec w = attrib(M, "isHomog"); // grading weights
+
+  intvec gr = grdeg( M );
   
 //  intvec d = deg(M[1..nocls(M)]); // no need to deal with un-weighted degrees??!
@@ -261 +262 @@
   intvec pivot = mysort(gr, s);
 
+  // grades & ordering permutation for N.  gr[pivot] should be sorted!
   ASSUME(1, issorted(gr[pivot], s));
   
-  return (gr, pivot);
-}
-
-
-static proc reorder(def M, int s, list #)
-"
-helper for    Reorder gens of M 
-"
-{
-  // use initial grading if given
-  if( size(#) == 1 )
-  {
-    def w = #[1];
-    if( typeof(w) == "intvec" )
-    {
-      if( nrows(M) == size(w) )
-      {
-        attrib(M, "isHomog", w);
-      }
-    }
-    kill w;
-  }
-
-  // try to homogenize otherwise 
-  if( typeof(attrib( M, "isHomog")) != "intvec" )
-  { // no such attribute? //   ERROR("No grading!");
-    ASSUME(0, /* input must be graded: cannot compute homogenizing grading */ homog(M) ); 
-  }
-
-  ASSUME(0, /* input must be graded! */ typeof(attrib(M, "isHomog")) == "intvec" );
-  // input should be graded:
-  def w = attrib(M, "isHomog"); // grading weights?
-  ASSUME(0, /* input must be correctly graded! */ nrows(M) == size(w) );   
- 
-
-  intvec d, p;
-
-  if( size(#) > 0 ) { (d, p) = order( M, s, #[1] ); }
-  else              { (d, p) = order( M, s       ); }  
-
-  // grades & ordering permutation for N.  d[p] should be sorted!
-  
-  def N = grobj(module(M[p]), w);  // reorder the starting ideal/module
-
-  d = d[p];
+  module N = grobj(module(M[pivot]), w);  // reorder the starting ideal/module
 
 //  "reorder: "; grview(N);
   
-  return (N, d);
-}
-
-/*
-// only for whole resolutions, please use grorder instead!
-static proc ordres( list L, list # )
-"
-reorder a resolution given by the list (and optionnaly a grading for its 1st entry)
-"
-{  
-  int k = 1; // "k: ", k; 
-
-  // check 1st syzygy property?
-//  if( 0 != size( module( matrix(transpose(L[k+1]))*matrix(transpose(L[k-1+1])) ) ) ){L[k];"";L[k+1];module( matrix(transpose(L[k+1]))*matrix(transpose(L[k-1+1])) );ERROR( "Exactness test failed!!!!" );}
-
-
-  // TODO: rewrite with reorder
-  intvec d, p, pp;
-  module N = L[1];
-  (d, p) = order( N, 1, # ); // grades & ordering permutation for N.  d[p] should be sorted!
-  N = module(N[p]);  // reorder the starting ideal/module
-  attrib( N, "isHomog", w ); // set the grading
-  L[1] = N; // put it back into the list of modules
-
-// grview(N);
- kill w, N;
- 
- ///////////////////////////////////
- k++; 
- while( k <= size(L) )
- {
-  // "k: ", k;
-  module N = L[k];
-  
-  if( size(N) == 0 ) { break; } // resolution should finish with 0 module?
-  
-  if( typeof(attrib( N, "isHomog")) != "intvec" )  
-  { 
-    attrib( N, "isHomog", d); // ERROR("Non-graded input!");
-  }
-  
-//  grview(N);
-
-  intvec w = attrib( N, "isHomog");  
-  (d, pp) = order( N, w, 1 ); // grades & ordering permutation : d[p] should be sorted!
-
-  // reorder k-th syzygy (columns):
-  module T = module(N[pp]);
-  kill N;
-  module N = transpose(T); 
-  kill T;
-
-  // reorder k-th syzygy (rows):
-  module T = module(N[p]);
-  kill N;
-  module N = transpose(T);
-  kill T;
-
-  w = intvec(w[p]); attrib( N, "isHomog", w); // corresponding ordered grading
-
-  L[k] = N; //  grview(N);
-
-  // check syz. property?
-//  if( 0 != size( module( matrix(transpose(N))*matrix(transpose(L[k-1])) ) ) ) { N; L[k-1];  module( matrix(transpose(N))*matrix(transpose(L[k-1])) ); ERROR( "Exactness test failed!!!!" );  } 
-  kill N, w;
-  p = pp;
-  k++;
- }
-
- return (list( L[1 .. (k-1)] ));
-}
-*/
-
-
-
-
-proc grtranspose(def M, list #)
+  return (N, intvec(gr[pivot]));
+}
+
+proc grtranspose(def M)
 "
 transpose graded module/map or a chain complex?...
@@ -417 +302 @@
 //      grview(M[i]);
       L[j] = grtranspose(M[i]);
-      if( i > 1 )
+//      grview(L[j]);
+
+      if( (i > 1) && (j > 0) )
       {
-        ASSUME(2, size( module( matrix(L[j+1])*matrix(L[j]) ) ) == 0 );
+//        grview(L[j+1]);
+        ASSUME(2, size( module( matrix(L[j])*matrix(L[j+1]) ) ) == 0 );
       };
 //      grview(L[j]);
@@ -429 +317 @@
 //////
 // "a";  grview(M);
+  ASSUME(1, grtest(M) );
 
   // TODO: something is wrong here:
@@ -434 +323 @@
  intvec d; module N;  
 
- if( size(#) > 0 ) { (N,d) = reorder(M, -1, #[1]); }
- else              { (N,d) = reorder(M, -1); }
+ (N,d) = reorder(M, -1);
 
  kill M; module M = grobj(transpose(N), -d);
@@ -456 +344 @@
 
 
-proc grorder(def M, list #)
+proc grorder(def M)
 "
 reorder graded module/map or a chain complex?...
@@ -493 +381 @@
     return (L); // ?
   }
+  
+  ASSUME(1, grtest(M) );
 
 // "a";  grview(M);
@@ -498 +388 @@
   intvec d; module N;
 
-  if( size(#) > 0 ) { (N,d) = reorder(M, 1, #[1]); }
-  else              { (N,d) = reorder(M, 1); }  
-  kill M; module M = grobj(transpose(N), -d);
+  (N,d) = reorder(M, 1); kill M; 
+  
+  module M = grobj(transpose(N), -d); kill N,d; 
 
 // "b";  grview(M);
  
-  kill N,d; module N; intvec d;
+  module N; intvec d;
   // reverse order:
   (N,d) = reorder(M, -1); kill M;
@@ -528 +418 @@
   " = Ring: ", string(basering);
 
-  I = groebner(I); attrib(I, "isHomog", intvec(0));
+  I = groebner(I); attrib(I, "isHomog", intvec(0)); 
+  ASSUME(0, grtest(I));
 //  " = Input degrees: "; grview(I);
 
@@ -553 +444 @@
 }
 example
-{ "EXAMPLE:"; echo = 2; int assumeLevel = 5;
+{ "EXAMPLE:"; echo = 2; 
+  if( defined(assumeLevel) ){ int assumeLevel0 = assumeLevel; } else { int assumeLevel; export(assumeLevel); }  // store the state of aL
+  assumeLevel = 5;
 
   string Name = "castelnuovo"; int @p=31991; ring R = (@p),(x,y,z,u,v), dp;ideal I = 5153xy2-98/23y3-101/51xyz+33/41y2z+99/79xz2+7136yz2-106/111z3+119/53xyu+34/57y2u-77/92xzu+84/73yzu-109/78z2u-27/56xu2+10023yu2+82/103zu2-34/25u3+3/2xyv-68/25y2v+12721xzv+4/63yzv-73/21z2v-7291xuv-91/53yuv-4/79zuv-34/91u2v-122/53xv2+123/70yv2-64/73zv2+44/65uv2+14/31v3,xy2-15202y3+10613xyz+13640y2z-107/103xz2+5292yz2+19/119z3-10042xyu+2770y2u+7957xzu+14008yzu+92/121z2u-92/51xu2+1178yu2+1/117zu2-12726u3+82/101xyv-92/17y2v-107/56xzv+14233yzv+79/28z2v+51/50xuv-31/5yuv+95/91zuv+19/108u2v+12151xv2-69/110yv2+37/89zv2-63/116uv2-88/23v3,-5153x2+37/23xy+8706y2-13160xz+68/115yz+5548z2-22/61xu-113/98yu+11818zu+2114u2-101/97xv+89/22yv-3355zv-113/5uv-5521v2;TestGRRes(Name, 2, I); kill R, Name, @p;  "";
@@ -574 +467 @@
 
   string Name = "rat.d10.g9.quart2"; int @p=31991; ring R = (@p),(x,y,z,u,v), dp;ideal I = x3yu2-48/11x2y2u2-8356xy3u2+35/121y4u2+31/66x3zu2-54/83x2yzu2-61/18xy2zu2+11526y3zu2+7372x2z2u2-91/60xyz2u2-95/97y2z2u2-45/71xz3u2+71/115yz3u2+25/54z4u2-61/102x3u3-12668x2yu3+6653xy2u3+41/54y3u3+87/50x2zu3-5004xyzu3+13924y2zu3+2310xz2u3-93/14yz2u3-2/93z3u3-97/125x2u4-58/11xyu4+46/73y2u4-4417xzu4+60/101yzu4+56/75z2u4-113/118xu5+115/4yu5-40zu5-8554u6-54/83x3yuv-9770x2y2uv-590xy3uv+15/49y4uv+94/69x3zuv+121/105x2yzuv+95/88xy2zuv+3186y3zuv+11/6x2z2uv-44/81xyz2uv+637y2z2uv+109/121xz3uv-33yz3uv-94/115z4uv-49/95x3u2v-11/109x2yu2v+45/113xy2u2v+97/84y3u2v+5257x2zu2v+99/49xyzu2v+12584y2zu2v-4294xz2u2v+1137yz2u2v-58/69z3u2v-4749x2u3v+120/97xyu3v-31/103y2u3v+62/97xzu3v-107/74yzu3v+53/59z2u3v+91/33xu4v+1291yu4v+23/34zu4v+58/77u5v+16/17x3yv2-750x2y2v2+86/89xy3v2+123/46y4v2+53/123x3zv2-61/99x2yzv2+12389xy2zv2+10419y3zv2+43/11x2z2v2-146xyz2v2-116/51y2z2v2+13/62xz3v2-5524yz3v2-111/118z4v2-56/55x3uv2-3038x2yuv2+14/27xy2uv2-43/64y3uv2+3385x2zuv2+25/11xyzuv2+92/41y2zuv2+28/113xz2uv2-2049yz2uv2+89/37z3uv2-13094x2u2v2-2774xyu2v2+15474y2u2v2-15791xzu2v2-71/116yzu2v2+77/41z2u2v2-83/68xu3v2-33/106yu3v2+71/37zu3v2-41/17u4v2+12052x3v3+1906x2yv3+13825xy2v3+80/7y3v3-125/96x2zv3-9661xyzv3+85/116y2zv3-72/91xz2v3+13/112yz2v3-126/97z3v3-1637x2uv3+34/103xyuv3+3844y2uv3+77/10xzuv3+6359yzuv3-11185z2uv3-124/121xu2v3+66/91yu2v3-14636zu2v3-1051u3v3+9/64x2v4-12924xyv4-119/41y2v4+74/23xzv4+1622yzv4+73/37z2v4-60/101xuv4+111/22yuv4-45/124zuv4+59/37u2v4-66/37xv5-71/99yv5+12409zv5-113/64uv5-5267v6,-x4y-22/79x3y2-125/42x2y3-116/7xy4+98/111y5-31/66x4z-118/75x3yz+110/93x2y2z-43/92xy3z-788y4z-7372x3z2-2701x2yz2-67/124xy2z2-117/62y3z2+45/71x2z3-8396xyz3-10343y2z3-25/54xz4+30/59yz4+61/102x4u+11736x3yu+12726x2y2u+41/118xy3u-15832y4u-87/50x3zu-130x2yzu+41/8xy2zu-10300y3zu-2310x2z2u-101/5xyz2u+6205y2z2u+2/93xz3u+8679yz3u+97/125x3u2-43/37x2yu2-39/80xy2u2+12139y3u2+4417x2zu2+4294xyzu2+11/58y2zu2-56/75xz2u2+8338yz2u2+113/118x2u3-10190xyu3-37/16y2u3+40xzu3+74/23yzu3+8554xu4+115/22yu4-39/79x4v+61/72x3yv+8048x2y2v-9201xy3v+16/121y4v+113/93x3zv+109/75x2yzv+12700xy2zv-10607y3zv+50/11x2z2v+1223xyz2v-103/79y2z2v-123/58xz3v+31/26yz3v-15/122z4v+122/25x3uv-99/17x2yuv+1723xy2uv-38/121y3uv+11016x2zuv-25/102xyzuv-14970y2zuv-61/6xz2uv-14981yz2uv+15900z3uv+3268x2u2v-75/19xyu2v-1436y2u2v-1764xzu2v-57/41yzu2v+12741z2u2v-14615xu3v+119/61yu3v-115/119zu3v+10501u4v-8502x3v2-51/76x2yv2-6281xy2v2+17/49y3v2-106/7x2zv2+63/101xyzv2-27/95y2zv2-1606xz2v2+9245yz2v2+1912z3v2+11155x2uv2+223xyuv2-13/18y2uv2+110/43xzuv2+76/81yzuv2-6291z2uv2+1400xu2v2-95/23yu2v2-9701zu2v2+106/105u3v2+72/47x2v3-13118xyv3+14409y2v3+37/86xzv3+44/69yzv3-325z2v3+113/71xuv3+16/81yuv3+6/19zuv3-119/39u2v3-89/9xv4+72/53yv4+112/55zv4-8587uv4-6604v5,-x3y2+48/11x2y3+8356xy4-35/121y5-12750x3yz+100/111x2y2z+45/74xy3z+99/74y4z-6/7x3z2-47/67x2yz2+11465xy2z2-11865y3z2+7776x2z3+124/45xyz3-98/115y2z3+117/85xz4-59/120yz4-8748z5+61/102x3yu+12668x2y2u-6653xy3u-41/54y4u+13408x3zu-2185x2yzu-1240xy2zu+1161y3zu+44/27x2z2u-11164xyz2u-13388y2z2u-107/13xz3u+90/71yz3u+4204z4u+97/125x2yu2+58/11xy2u2-46/73y3u2+55/48x2zu2+121/31xyzu2+126/61y2zu2-55/69xz2u2+5988yz2u2+3755z3u2+113/118xyu3-115/4y2u3+3390xzu3-5762yzu3+30/61z2u3+8554yu4-14317zu4+99/116x3yv-113/119x2y2v+50/23xy3v-37/79y4v-8668x3zv+14049x2yzv+111/35xy2zv+61/28y3zv-10171x2z2v+68/21xyz2v+2023y2z2v-9/109xz3v+8520yz3v-2683z4v-13547x3uv+28/65x2yuv-5988xy2uv+61/111y3uv+12314x2zuv+29/44xyzuv+6141y2zuv+11280xz2uv+79/22yz2uv-38/111z3uv+19/51x2u2v+5093xyu2v-10291y2u2v-5009xzu2v-111/49yzu2v+3813z2u2v-61/37xu3v+15914yu3v-3218zu3v-12915u4v-118/101x3v2-7/57x2yv2+13128xy2v2+11606y3v2+42/101x2zv2-54/17xyzv2-43/49y2zv2-119/110xz2v2+9742yz2v2-43/4z3v2-55/8x2uv2-29/88xyuv2+12042y2uv2+101/37xzuv2-57/62yzuv2+106/97z2uv2+38/83xu2v2+8152yu2v2-5492zu2v2-47/79u3v2+15112x2v3+69/44xyv3-6/71y2v3+113/54xzv3-13210yzv3-707z2v3-119/8xuv3+3845yuv3-19/20zuv3+4852u2v3+15761xv4-12372yv4+74/69zv4-2100uv4-12833v5,-x3yz+48/11x2y2z+8356xy3z-35/121y4z-31/66x3z2+54/83x2yz2+61/18xy2z2-11526y3z2-7372x2z3+91/60xyz3+95/97y2z3+45/71xz4-71/115yz4-25/54z5+15/52x3yu+6039x2y2u+74/99xy3u-17/40y4u+29/50x3zu-7775x2yzu+6368xy2zu+14170y3zu+52/41x2z2u+7003xyz2u-5787y2z2u-101/37xz3u-23/28yz3u-20/63z4u+41/77x3u2+8650x2yu2-15922xy2u2-16/83y3u2+7278x2zu2+31/30xyzu2-2/107y2zu2+35/122xz2u2+85/58yz2u2-757z3u2+2/101x2u3+86/17xyu3+95/59y2u3+123/22xzu3-6869yzu3-9311z2u3-105/97xu4+5699yu4+15925zu4+13528u5-154x3yv+4187x2y2v+56/107xy3v-15932y4v-5137x3zv-37/56x2yzv+9401xy2zv+92/123y3zv-79/97x2z2v+9201xyz2v+19/53y2z2v+107/20xz3v+17/77yz3v-15306z4v+3215x3uv-79/117x2yuv-9/76xy2uv-6352y3uv+93/13x2zuv-65/89xyzuv-115/4y2zuv-34/57xz2uv+39/107yz2uv+31/9z3uv+107/48x2u2v+2632xyu2v+29/96y2u2v-125/89xzu2v+29/113yzu2v+3940z2u2v-116/111xu3v+6145yu3v-105/62zu3v+101/17u4v-9281x3v2-49/107x2yv2-12154xy2v2+4/19y3v2-114/71x2zv2-15/118xyzv2+4372y2zv2+45/121xz2v2+46/111yz2v2+6614z3v2+17x2uv2+10806xyuv2-10617y2uv2-25/111xzuv2-116/27yzuv2-7/58z2uv2-686xu2v2+3/13yu2v2-17/49zu2v2-40/107u3v2+47/90x2v3-83/43xyv3-6326y2v3+49/64xzv3+113/76yzv3-122/73z2v3+10232xuv3-116/109yuv3-1990zuv3+70/51u2v3-118/19xv4-27/55yv4+21/19zv4-23/57uv4-11721v5,-3399x4y+1849x3y2-3/29x2y3+28/87xy4+10/29y5-9788x4z-49/73x3yz+13829x2y2z+118/73xy3z+13129y4z-618x3z2+92/13x2yz2+101/117xy2z2-162y3z2+24/5x2z3-29/74xyz3+2687y2z3-74/39xz4+2/57yz4+68/73x4u-13787x3yu-11659x2y2u+14729xy3u+92/53y4u+15/71x3zu-62/15x2yzu+21/85xy2zu+4938y3zu-120/37x2z2u-77/102xyz2u-4785y2z2u-83/70xz3u-12128yz3u-13592z4u-123/20x3u2+2607x2yu2+40/19xy2u2+6361y3u2-3091x2zu2+89/113xyzu2+149y2zu2-2890xz2u2-8374yz2u2+11886z3u2-49/43x2u3-9854xyu3-6943y2u3+10743xzu3-122/45yzu3-13902z2u3-103/19xu4-48/59yu4+27/86zu4+46/35u5-117/17x4v-15/7x3yv+8409x2y2v-83/28xy3v+86/35y4v+37/45x3zv+4/3x2yzv+35/38xy2zv+4015y3zv-49/111x2z2v-1260xyz2v-25/33y2z2v+116/19xz3v+93/8yz3v+5755z4v-25/89x3uv-11669x2yuv-64/107xy2uv+2993y3uv+7767x2zuv-17/95xyzuv-103/80y2zuv-14576xz2uv+80/47yz2uv+25/107z3uv+103/2x2u2v+125/117xyu2v-2/89y2u2v-5298xzu2v-50/27yzu2v-71/53z2u2v+2652xu3v+15761yu3v+2124zu3v+11/82u4v+100/63x3v2+4180x2yv2+11/39xy2v2-1221y3v2+108/125x2zv2+97/126xyzv2-7698y2zv2+13984xz2v2+1342yz2v2-84/121z3v2-26/73x2uv2-14/15xyuv2-22/37y2uv2-71/82xzuv2+12430yzuv2+103/52z2uv2-13095xu2v2+10114yu2v2-8/73zu2v2-33/97u3v2+83/105x2v3+22/45xyv3-7961y2v3-9654xzv3-54/55yzv3-3/71z2v3-10148xuv3-117/98yuv3+101/102zuv3-606u2v3+97/43xv4-68/21yv4+63/16zv4+42/17uv4+5834v5,-3399x3y2-32/113x2y3+14/99xy4+15001y5-121/115x3yz+4604x2y2z+7/2xy3z+9532y4z-3267x3z2+97/118x2yz2-14238xy2z2-80/21y3z2-12332x2z3-19/69xyz3+116/15y2z3-103/32xz4+15340yz4+10509z5+112/109x3yu-97x2y2u-40/11xy3u+90/29y4u-95/106x3zu-114/67x2yzu+113/48xy2zu+12080y3zu-44x2z2u+18/17xyz2u-4814y2z2u-103/100xz3u-96/61yz3u-205z4u-87/82x3u2-97/108x2yu2+3230xy2u2+104/83y3u2+41/86x2zu2+116/49xyzu2-59/110y2zu2+14/59xz2u2-6962yz2u2-2185z3u2+59/91x2u3+2497xyu3+3/37y2u3-13010xzu3+6/83yzu3-11448z2u3+13/72xu4-69/62yu4-2869zu4+23/73u5-20/43x3yv+5074x2y2v+28/125xy3v-2706y4v+13010x3zv-17/109x2yzv+21/4xy2zv+59/93y3zv-2406x2z2v+117/11xyz2v-14978y2z2v+70/89xz3v-33/7yz3v-13676z4v-13690x3uv+9825x2yuv-117/107xy2uv+12760y3uv-93/98x2zuv-113/64xyzuv+113/103y2zuv-9748xz2uv+11016yz2uv-10729z3uv+90/13x2u2v-13/47xyu2v-11/39y2u2v+20/69xzu2v+5531yzu2v+125/49z2u2v-11025xu3v-9621yu3v+113/109zu3v+4710u4v-107/7x3v2+110/119x2yv2-10025xy2v2-6644y3v2-5041x2zv2+5/96xyzv2+11472y2zv2-5128xz2v2+2927yz2v2+121/18z3v2-125/89x2uv2+12936xyuv2-71/47y2uv2+34/47xzuv2-75/103yzuv2-2654z2uv2-2350xu2v2-7707yu2v2+47/72zu2v2-952u3v2-21/67x2v3+58/37xyv3-8757y2v3+3615xzv3+44/123yzv3-13027z2v3-9/10xuv3+75/43yuv3+115/18zuv3+8071u2v3-26/3xv4-67/65yv4+14186zv4-41/122uv4+33/28v5,-3399x3yz-32/113x2y2z+14/99xy3z+15001y4z-9788x3z2+37/96x2yz2+7743xy2z2+31/55y3z2-618x2z3-8171xyz3+82/109y2z3+24/5xz4+88/85yz4-74/39z5-13165x3yu+3407x2y2u-12509xy3u-23/45y4u-11774x3zu-10/67x2yzu+69/79xy2zu-10/123y3zu-7636x2z2u+83/32xyz2u+51/112y2z2u+19/8xz3u+9309yz3u-44/49z4u+4089x3u2-374x2yu2-919xy2u2+98/107y3u2+2776x2zu2+85/26xyzu2+31/13y2zu2-103/82xz2u2+35/76yz2u2+59/45z3u2+2950x2u3+27/44xyu3+88/71y2u3+7/114xzu3-72/77yzu3+12917z2u3-34/67xu4-85/82yu4-55/84zu4+4690u5+11/42x3yv-19/125x2y2v-8288xy3v+9199y4v-12929x3zv+13357x2yzv-4903xy2zv-584y3zv-10/33x2z2v+59/113xyz2v+103/92y2z2v+101/69xz3v+8708yz3v-8/7z4v+13560x3uv-43/49x2yuv-121/98xy2uv+75/79y3uv-39x2zuv-88/69xyzuv-89/78y2zuv+110/67xz2uv+61/4yz2uv-98/45z3uv+82/7x2u2v-85/41xyu2v+6548y2u2v+9367xzu2v-59/81yzu2v-14408z2u2v+2363xu3v-80/11yu3v-50/17zu3v-14799u4v-53/21x3v2+9437x2yv2-117/80xy2v2+81/85y3v2-8/45x2zv2-6428xyzv2+15126y2zv2+68/89xz2v2+7/122yz2v2+9639z3v2+113/4x2uv2-8678xyuv2-104/45y2uv2-79/90xzuv2+39/101yzuv2-7234z2uv2-28/43xu2v2+1251yu2v2-97/56zu2v2+17/41u3v2+107/24x2v3+2747xyv3+9933y2v3-4199xzv3+53/83yzv3+6364z2v3-5456xuv3+618yuv3-123/55zuv3+2375u2v3+63/76xv4-115/106yv4-8811zv4-31/75uv4+10/109v5,13/89x4y+77/31x3y2+36/83x2y3-11411xy4+6936y5-12223x4z+7400x3yz+33/118x2y2z-12146xy3z+108/79y4z+82/99x3z2-9877x2yz2-79/70xy2z2-19/123y3z2-1491x2z3+7953xyz3-43/126y2z3+60/17xz4+98/57yz4-13317x4u-77/27x3yu-6811x2y2u-69/61xy3u+6144y4u+5404x3zu+121/120x2yzu-91/23xy2zu-71/106y3zu+1435x2z2u-120/13xyz2u-12019y2z2u-68/7xz3u-113/82yz3u+11526z4u-8706x3u2-89/53x2yu2-14804xy2u2+120/107y3u2+71/94x2zu2-1/70xyzu2+1532y2zu2+4470xz2u2+13/60yz2u2-115/102z3u2-82/21x2u3+27/121xyu3-4439y2u3-101/47xzu3-3186yzu3-106/101z2u3-10169xu4+19/58yu4-96/73zu4-7959u5-10526x4v-107/92x3yv+47/6x2y2v-23/43xy3v-69/62y4v+59/65x3zv-28/95x2yzv+5479xy2zv-39/77y3zv+11/69x2z2v-11713xyz2v+43/79y2z2v-15602xz3v+16/73yz3v-13952z4v+61/82x3uv-2219x2yuv-91/106xy2uv+5/37y3uv-148x2zuv+31/51xyzuv+18/101y2zuv+97/68xz2uv-73/32yz2uv+47/2z3uv+2/41x2u2v-13009xyu2v-7/60y2u2v+15779xzu2v+72/7yzu2v-11/73z2u2v-119/44xu3v-9067yu3v+3249zu3v+61/51u4v+12525x3v2-118/9x2yv2-3270xy2v2-4/25y3v2-5075x2zv2+77/40xyzv2-89/65y2zv2+17/58xz2v2-15609yz2v2+95/54z3v2-75/79x2uv2-4907xyuv2+12418y2uv2-57/17xzuv2-8746yzuv2+13/95z2uv2-124/67xu2v2+16/13yu2v2+28/23zu2v2-10847u3v2-645x2v3+106/75xyv3+6/115y2v3-8495xzv3+58/35yzv3-9398z2v3-101/72xuv3-71/20yuv3-124/65zuv3-8971u2v3+27/28xv4+12/29yv4-4276zv4+10858uv4+29/12v5,13/89x3y2+12068x2y3-15543xy4-77/79y5+6626x3yz+64/53x2y2z-6/23xy3z-47/125y4z+14403x3z2-43/78x2yz2-31/115xy2z2+94/59y3z2-118/117x2z3-11229xyz3+2268y2z3-116/85xz4+25/58yz4+3085z5+59/27x3yu+67/82x2y2u+11/6xy3u+103/47y4u-63/80x3zu-81/47x2yzu+7760xy2zu-115/56y3zu-10/17x2z2u+101/5xyz2u+15634y2z2u+1/107xz3u-9282yz3u+43/62z4u+62/55x3u2+100/113x2yu2-9205xy2u2-46/13y3u2+43/96x2zu2+10159xyzu2+692y2zu2+859xz2u2-19/74yz2u2+123/47z3u2-9/20x2u3-11391xyu3-2375y2u3+109/24xzu3-57/53yzu3-925z2u3-82/45xu4+97/34yu4+13/82zu4-108/29u5+63/10x3yv+38/17x2y2v-19/115xy3v+3150y4v+22/69x3zv+26/57x2yzv+110/27xy2zv+87/77y3zv+85/18x2z2v+39/47xyz2v-48/17y2z2v-7/27xz3v-13/100yz3v-11662z4v-17/8x3uv+37/11x2yuv+29/11xy2uv-109/88y3uv-2817x2zuv-61/44xyzuv+10/31y2zuv+10010xz2uv+51/86yz2uv-97/83z3uv-89/96x2u2v+4030xyu2v-58/77y2u2v-114/43xzu2v-37/10yzu2v-2011z2u2v+14483xu3v-109/101yu3v+121/102zu3v-79/92u4v+15113x3v2+10781x2yv2-14259xy2v2-113/48y3v2-7/94x2zv2-17/74xyzv2-5/117y2zv2-59/75xz2v2+13188yz2v2+103/43z3v2+4/125x2uv2-52/59xyuv2+85/92y2uv2-1/46xzuv2-9106yzuv2-83/11z2uv2-23/94xu2v2+6742yu2v2-35/107zu2v2-14596u3v2-117/43x2v3+1026xyv3+90/19y2v3+14671xzv3-101/100yzv3+6962z2v3+61/68xuv3+108/37yuv3-4157zuv3-3974u2v3+15677xv4+8661yv4+8459zv4-16/23uv4-37/119v5,13/89x3yz+12068x2y2z-15543xy3z-77/79y4z-12223x3z2-13941x2yz2+115/84xy2z2+13/98y3z2+82/99x2z3+7751xyz3+122/17y2z3-1491xz4+1327yz4+60/17z5+15363x3yu+9780x2y2u+19/117xy3u-1924y4u-14600x3zu+46/41x2yzu-5466xy2zu-73/12y3zu+10838x2z2u-8302xyz2u-89/113y2z2u+53/69xz3u-9224yz3u+47/33z4u-7399x3u2+89/77x2yu2+9312xy2u2-41/80y3u2-732x2zu2-6781xyzu2-8608y2zu2-9270xz2u2-117/58yz2u2-115/68z3u2-48/31x2u3-9067xyu3+97/107y2u3+73/57xzu3-2719yzu3-110/59z2u3-37/86xu4-15796yu4-61/4zu4-115/72u5+6161x3yv+4134x2y2v+677xy3v-8375y4v+1150x3zv+1551x2yzv+4157xy2zv+112/87y3zv+8171x2z2v+6040xyz2v+15651y2z2v-7/66xz3v-47/61yz3v+77/64z4v+14848x3uv+48/119x2yuv-9534xy2uv-117/95y3uv+5/4x2zuv+122xyzuv+90/31y2zuv-41/26xz2uv+31/30yz2uv-10428z3uv-9896x2u2v-71/21xyu2v-55/38y2u2v-29/22xzu2v-11092yzu2v+39/122z2u2v+93/73xu3v+22/49yu3v-21/106zu3v+56u4v+8565x3v2-1695x2yv2+2/17xy2v2+1/78y3v2-113/71x2zv2-41/100xyzv2+55/14y2zv2+15286xz2v2+17/53yz2v2+126/71z3v2-79/87x2uv2+109/97xyuv2-28/31y2uv2-6533xzuv2+22/5yzuv2-10449z2uv2+10830xu2v2-15516yu2v2+28/57zu2v2-81/22u3v2+4198x2v3+5667xyv3-7133y2v3-8408xzv3+11066yzv3-26/125z2v3-808xuv3+95/54yuv3-64/17zuv3-5267u2v3-15333xv4+42/89yv4+63/85zv4+119/113uv4-2011v5,5583x4y+1725x3y2-8652x2y3-91/25xy4-8495x4z-13731x3yz+9298x2y2z-41/111xy3z-15503y4z-13805x3z2+3962x2yz2-2/63xy2z2+3314y3z2+2522x2z3-10/87xyz3-408y2z3+7/16xz4+69/22yz4-7254z5-59/21x4u+115/7x3yu-1718x2y2u+7851xy3u+2632y4u-82/3x3zu+37/86x2yzu+101/113xy2zu+6747y3zu-109/113x2z2u+7399xyz2u+24/103y2z2u+89/9xz3u-14630yz3u+15066z4u-12561x3u2+113/115x2yu2+87/97xy2u2-126/67y3u2-48/7x2zu2+123/103xyzu2-11/107y2zu2-2747xz2u2+8158yz2u2-3/107z3u2+41/6x2u3+12767xyu3+3873y2u3+74/83xzu3-55/119yzu3-24/83z2u3+55xu4-7/95yu4+57/44zu4+2/101u5-6928x4v-121/57x3yv+111/104x2y2v+946xy3v-29y4v+3057x3zv-14/25x2yzv+43/31xy2zv-105/2y3zv+2336x2z2v+61/77xyz2v-7880y2z2v+5/58xz3v+10593yz3v+7094z4v+63/59x3uv-5/69x2yuv-11/81xy2uv-4157y3uv+73/65x2zuv-1676xyzuv-2376y2zuv-85/63xz2uv-95/2yz2uv-14903z3uv-119/110x2u2v-115/24xyu2v+125/9y2u2v+106/87xzu2v-13/12yzu2v-4/19z2u2v+7838xu3v-43/111yu3v+7/113zu3v-12500u4v+7743x3v2-2023x2yv2-85/83xy2v2+49/41y3v2+20/87x2zv2+3932xyzv2-77/6y2zv2+47/90xz2v2-15580yz2v2+39/4z3v2-61/8x2uv2+2518xyuv2+29/98y2uv2+11057xzuv2-18/107yzuv2+708z2uv2+14720xu2v2-3175yu2v2-113/59zu2v2-14735u3v2+7/69x2v3-4029xyv3+54/91y2v3+12372xzv3+67/2yzv3+8856z2v3-2178xuv3+995yuv3+64/95zuv3+4039u2v3-37/44xv4+23/17yv4-3035zv4-103/124uv4+69/64v5,-5583x3y2-1725x2y3+8652xy4+91/25y5+6201x3yz-73/49x2y2z-3844xy3z+10548y4z-11057x3z2-105/122x2yz2+31/53xy2z2+79/89y3z2-24/101x2z3+107/119xyz3-126y2z3+8164xz4+2/77yz4-51/8z5-14941x3yu-106x2y2u+8695xy3u+125/62y4u+4328x3zu+29/117x2yzu-6249xy2zu-2791y3zu+67/49x2z2u-38/29xyz2u+122/41y2z2u+10603xz3u-3029yz3u+5578z4u+14754x3u2-108/79x2yu2+4408xy2u2-12401y3u2-1426x2zu2-1741xyzu2-83/86y2zu2+79/95xz2u2+122/121yz2u2+81/2z3u2-1172x2u3-41/68xyu3-70/3y2u3+24/107xzu3+120/79yzu3+18/119z2u3-65/122xu4+1018yu4+22/107zu4+15189u5+5/8x3yv-12060x2y2v+3/62xy3v-227y4v+60/41x3zv-123/115x2yzv+110/123xy2zv+12864y3zv-86/121x2z2v-69/94xyz2v+14/79y2z2v+118/45xz3v+10842yz3v-37/58z4v+100/69x3uv-47/65x2yuv-7/67xy2uv-93/100y3uv-6262x2zuv-4/75xyzuv+2082y2zuv-9117xz2uv+12450yz2uv-84/67z3uv+123/26x2u2v-51/89xyu2v+19/74y2u2v-104/77xzu2v+318yzu2v+12402z2u2v+95/8xu3v-81/26yu3v-4486zu3v+3872u4v+72/91x3v2-83/63x2yv2+93/92xy2v2-15924y3v2-53/62x2zv2+6046xyzv2+1408y2zv2+60/107xz2v2-1150yz2v2-126/19z3v2-7429x2uv2+2554xyuv2+3602y2uv2+10738xzuv2-57/64yzuv2+86/69z2uv2+8172xu2v2+91/113yu2v2+92/65zu2v2+118/37u3v2+47/83x2v3+12750xyv3+10851y2v3+4216xzv3+6/101yzv3-108z2v3+2920xuv3-101/102yuv3-157zuv3+7742u2v3-7234xv4-2/111yv4+59/33zv4-93/91uv4+24/19v5,1592x4y+75/121x3y2+40/19x2y3-2651xy4+9934x4z+245x3yz+11665x2y2z+30/41xy3z+1823y4z+89/88x3z2-105/46x2yz2+79/58xy2z2-4191y3z2-76/61x2z3-21/32xyz3-9516y2z3-14896xz4-85/77yz4+51/109z5+61/30x4u-10/101x3yu+11796x2y2u+76/101xy3u+123/88y4u-5932x3zu-11857x2yzu+7128xy2zu-45/79y3zu+119/18x2z2u+9/74xyz2u+7042y2z2u-1114xz3u-11/82yz3u-1466z4u-6/85x3u2+27/106x2yu2+14246xy2u2-6216y3u2+47/6x2zu2-45/59xyzu2+89/41y2zu2+41/80xz2u2-7583yz2u2-75/113z3u2-14808x2u3-10873xyu3-90/67y2u3-11081xzu3-7369yzu3-7131z2u3-1402xu4-15386yu4-108/73zu4-5039u5+120/113x4v+10617x3yv-50/87x2y2v-2395xy3v-20/69y4v-8587x3zv+12960x2yzv-41/50xy2zv-13844y3zv-65/32x2z2v-77/122xyz2v-85/66y2z2v+13/100xz3v-20/51yz3v-13676z4v+76/97x3uv+1046x2yuv-8059xy2uv-117/59y3uv-29/105x2zuv+7287xyzuv-107/119y2zuv-35/118xz2uv+79/86yz2uv-2211z3uv+5448x2u2v+62/35xyu2v-2275y2u2v+29/121xzu2v-1674yzu2v-56/43z2u2v-3377xu3v-43/110yu3v+23/10zu3v-24/61u4v+121/53x3v2-4745x2yv2-57/64xy2v2+9554y3v2-12741x2zv2+10449xyzv2+37/108y2zv2+8621xz2v2-11/57yz2v2+1566z3v2+125/49x2uv2-121/118xyuv2+109/84y2uv2-335xzuv2+10167yzuv2-59/109z2uv2-103/119xu2v2+43/13yu2v2-73/87zu2v2+2037u3v2+13002x2v3+83/48xyv3-10713y2v3+1026xzv3-105/64yzv3-37/6z2v3+14779xuv3-6448yuv3+19/69zuv3-1/110u2v3+10010xv4+79/12yv4+12/19zv4-35/61uv4-11/57v5,-1592x3y2-75/121x2y3-40/19xy4+2651y5+39/121x3yz+122/77x2y2z-114/31xy3z+1544y4z+2/3x3z2-10271x2yz2-8373xy2z2+56/61y3z2+55/48x2z3-116xyz3-25/7y2z3-108/113xz4-34/53yz4+5548z5-122x3yu-9690x2y2u+43/87xy3u-5/19y4u+97/54x3zu-17/19x2yzu+4355xy2zu+12/5y3zu-1/100x2z2u+12754xyz2u+13600y2z2u+17/45xz3u-12091yz3u+5145z4u-63/64x3u2-84/31x2yu2-97/41xy2u2+7/13y3u2-79/62x2zu2-80/103xyzu2-69/14y2zu2+119/4xz2u2-35/87yz2u2-13840z3u2+14101x2u3+7952xyu3-1857y2u3-9861xzu3+3180yzu3+75/107z2u3-250xu4-15134yu4+4717zu4-2/41u5+22/27x3yv-8983x2y2v+10520xy3v-113/2y4v+10/73x3zv-1986x2yzv-110/13xy2zv+1550y3zv+32/111x2z2v-111/35xyz2v+101/98y2z2v+8045xz3v-2/89yz3v+2924z4v-79/11x3uv-15178x2yuv+10874xy2uv+54/11y3uv-8950x2zuv+70/53xyzuv-2403y2zuv-8249xz2uv+6935yz2uv+20/89z3uv+885x2u2v-76/71xyu2v-4/17y2u2v-31/52xzu2v-4/99yzu2v+10333z2u2v-93/104xu3v+82/101yu3v-71/37zu3v+9397u4v-15/112x3v2-6614x2yv2+119/2xy2v2+88/119y3v2+306x2zv2+2790xyzv2+10992y2zv2-115/74xz2v2-14711yz2v2+11612z3v2-1788x2uv2-75/97xyuv2+79/30y2uv2+99/59xzuv2-11439yzuv2-121/113z2uv2+108/37xu2v2+37/36yu2v2-3/65zu2v2-55/42u3v2+13/100x2v3-209xyv3-1272y2v3-117/68xzv3+63/94yzv3+32/59z2v3+1013xuv3-3463yuv3+6946zuv3-37/86u2v3+67/117xv4+85/28yv4-3024zv4-82/9uv4-32/65v5,-35/52x4y-12140x3y2+23/83x2y3+69/5xy4-80/79y5+120/43x4z-11865x3yz-3487x2y2z+53/59xy3z+53/102y4z-14083x3z2-14430x2yz2-2442xy2z2-33/104y3z2-91/38x2z3+4/87xyz3-26/57y2z3+4097xz4-9/122yz4+6364z5+9634x4u-97/95x3yu-46/99x2y2u+3847xy3u+121/106y4u+12765x3zu-5292x2yzu+1607xy2zu-67/121y3zu-12/35x2z2u+4/55xyz2u-17/27y2z2u+91/122xz3u-23/31yz3u+65/49z4u+73/46x3u2-124/27x2yu2-9933xy2u2+46/75y3u2+53/114x2zu2+3503xyzu2-14147y2zu2-11283xz2u2+11889yz2u2+99/104z3u2+3117x2u3+12624xyu3-10060y2u3+2193xzu3-80/47yzu3-77/13z2u3+11/31xu4-47/90yu4+49/48zu4-2/105u5-92/61x4v+7443x3yv+35/76x2y2v+114/67xy3v-73/126y4v+97/107x3zv+9464x2yzv+10869xy2zv+15718y3zv-37/33x2z2v+124/13xyz2v-11/26y2z2v-61/40xz3v+91/100yz3v-18/103z4v+60/29x3uv+21/125x2yuv-11117xy2uv+11748y3uv-16/117x2zuv+18/103xyzuv-1711y2zuv+1872xz2uv-109/123yz2uv-18/113z3uv-26/103x2u2v+14140xyu2v+11065y2u2v+8686xzu2v-5/111yzu2v+30/101z2u2v-10501xu3v-36/113yu3v-73/74zu3v+12753u4v-43/52x3v2-76/15x2yv2-5793xy2v2+18/13y3v2+1/79x2zv2+84/23xyzv2-172y2zv2+86/77xz2v2+15/37yz2v2+11835z3v2-6482x2uv2+94/113xyuv2+10727y2uv2-102/41xzuv2+15914yzuv2-12973z2uv2-9038xu2v2-13107yu2v2+1533zu2v2+12549u3v2-13528x2v3+903xyv3+23/114y2v3-123/64xzv3-81/5yzv3+111/103z2v3+4734xuv3-33/20yuv3-7954zuv3-2478u2v3+15518xv4-6723yv4-14/31zv4-3482uv4+10919v5,-3/94x4y-12936x3y2+2/11x2y3+32/23xy4-15921y5+61/93x4z+82/111x3yz-93/2x2y2z-6659xy3z-97/90y4z+402x3z2-14586x2yz2-121/39xy2z2+68/7y3z2+1212x2z3-2980xyz3+49/52y2z3-72/89xz4+92/47yz4+8478z5+2733x4u-103/89x3yu+1166x2y2u-7/53xy3u-106/23y4u+677x3zu+907x2yzu+7891xy2zu-9014y3zu+76/47x2z2u+49/116xyz2u-49/78y2z2u+12261xz3u+118/105yz3u-126/13z4u-8812x3u2-97/120x2yu2-9534xy2u2+92/5y3u2-54/71x2zu2+94/103xyzu2+2256y2zu2+4182xz2u2-5798yz2u2-31/115z3u2-73/98x2u3+15822xyu3+1004y2u3-578xzu3+9494yzu3-6779z2u3+14506xu4+10/121yu4+58/27zu4-2817u5-19/119x4v+7128x3yv+75/64x2y2v-65/109xy3v+5129y4v-53/55x3zv+54/125x2yzv-3009xy2zv+6144y3zv+15601x2z2v+123/55xyz2v-58/77y2z2v-56/61xz3v+121/10yz3v-103/86z4v-93/25x3uv+94/123x2yuv-25/107xy2uv+14807y3uv+65/7x2zuv+87/44xyzuv+6605y2zuv+23/99xz2uv-413yz2uv-17/15z3uv-79/46x2u2v+15240xyu2v-42/67y2u2v+8932xzu2v-5888yzu2v-4204z2u2v+7002xu3v-36/97yu3v-1634zu3v+61/102u4v-14/33x3v2-6520x2yv2+9004xy2v2-67/36y3v2-7/8x2zv2-24/11xyzv2-9373y2zv2+1556xz2v2-79/74yz2v2-6691z3v2+108x2uv2-76/61xyuv2+220y2uv2-1191xzuv2-4/9yzuv2+4546z2uv2+12205xu2v2+9/22yu2v2+64/93zu2v2-44/125u3v2+292x2v3+41/74xyv3+16/79y2v3-15892xzv3+5733yzv3+6796z2v3-42/55xuv3+71/79yuv3-19/104zuv3-38/15u2v3+6436xv4+28/15yv4+87/55zv4+2270uv4-30/41v5,-117/4x3y+97/122x2y2-3618xy3+6566y4+97/113x3z-12634x2yz+9865xy2z-1764y3z+114/31x2z2+5006xyz2+7/44y2z2-15040xz3+8/125yz3+11134z4-12980x3u-79/41x2yu-79/98xy2u+89/65y3u-1217x2zu+89/87xyzu+83/66y2zu+115/11xz2u+123/107yz2u+10920z3u-86/73x2u2-11/94xyu2-14054y2u2+6752xzu2-123/124yzu2+12129z2u2-13310xu3-52/63yu3+12847zu3-1545u4-11064x3v+11499x2yv-37/64xy2v+50/103y3v+123/94x2zv-126xyzv-111/44y2zv+95/14xz2v+113/83yz2v-77/103z3v+41/64x2uv+91/90xyuv-4932y2uv+103/31xzuv+62/63yzuv+1161z2uv-99/106xu2v-3181yu2v-11741zu2v-33/8u3v-3/118x2v2-9369xyv2+527y2v2-113/39xzv2-88/49yzv2-113/101z2v2+95/68xuv2-5930yuv2-20/43zuv2+7/41u2v2+109/93xv3-107/61yv3-8352zv3-5255uv3+12021v4,-2159x4-94/3x3y-4602x2y2+1609xy3+10721y4+28/9x3z-99/35x2yz+1/110xy2z+113/114y3z-118/75x2z2-103/93xyz2-68/67y2z2+13687xz3-1531yz3+61/107z4+6076x3u+9004x2yu+2211xy2u+110/53y3u+47/102x2zu+8495xyzu-9238y2zu+57/121xz2u-8543yz2u+8/19z3u-13527x2u2-13293xyu2+1138y2u2+26/115xzu2+78/53yzu2-12556z2u2+7299xu3+70/19yu3-14687zu3+13559u4+113/9x3v-85/126x2yv-83/3xy2v-3/46y3v+1814x2zv+28/79xyzv+103/51y2zv+78/31xz2v-14387yz2v+1/88z3v+116/75x2uv-101/59xyuv-70/3y2uv+109/71xzuv+13/88yzuv-147z2uv-113/76xu2v-9661yu2v+13855zu2v-6162u3v-1857x2v2-8208xyv2-4634y2v2-6178xzv2-7352yzv2-8247z2v2-113/15xuv2+99/40yuv2+21/97zuv2+11/37u2v2-6605xv3+8964yv3+35/121zv3+8543uv3-6008v4;TestGRRes(Name, 2, I); kill R, Name, @p; "";
+
+  if( defined(assumeLevel0) ){ assumeLevel = assumeLevel0; } else { kill assumeLevel; } // restore the state of aL
 }
 
 /////////////////////////////////////////////////////////
 
-// ????
+// Q@Woflram?
 proc grzero()
 "presentation of S(0)^1
@@ -584 +479 @@
 "
 {
- module Z = 0;
- return ( grobj(Z,intvec(0)) );
+ return ( grobj(module([0]), intvec(0)) );
 }
 
@@ -593 +487 @@
 "
 {
-  if(p==0){ return ( grzero() ); } // just ERROR ???
+  if(p==0){ ERROR("Sorry, we don't know what is A^0!?!?"); } // grzero!?
 
   ASSUME(0, p > 0);
+  ASSUME(1, grtest(A) );
 
   if(p==1){ return(A); }
@@ -617 +512 @@
 "
 {
+  ASSUME(1, grtest(A) );
+  ASSUME(1, grtest(B) );
+  
   intvec a = attrib(A, "isHomog");
   intvec b = attrib(B, "isHomog");
@@ -622 +520 @@
   int r = nrows(A);
   
-  ASSUME( 0, r == size(a) );
-  
-  module T; T[r] = 0; T = T, module(transpose(B));
-  module AB = module(A), transpose(T);
-  
-  return(grobj(AB, c));
+  module T = align(module(B), r); //  T;  print(T);  nrows(T); // BUG!!!!
+  module S = module(A), T;
+  
+  return(grobj(S, c));
 }
 example
-{ "EXAMPLE:"; echo = 2; int assumeLevel = 5;
-   
+{ "EXAMPLE:"; echo = 2; 
+
+  if( defined(assumeLevel) ){ int assumeLevel0 = assumeLevel; } else { int assumeLevel; export(assumeLevel); }  // store the state of aL
+  assumeLevel = 5;
+  
   ring r=32003,(x,y,z),dp;
 
@@ -681 +580 @@
   grview(T);
 
+  if( defined(assumeLevel0) ){ assumeLevel = assumeLevel0; } else { kill assumeLevel; } // restore the state of aL
 }
 
@@ -689 +589 @@
 {
   intvec a = attrib(M, "isHomog");
-  return (grobj(M, intvec( a + intvec(d:size(a))) ));
+  return (grobj(M, intvec( a - intvec(d:size(a))) ));
+}
+
+
+proc grtest(def N)
+"test if N is a matrix/module with matching isHomog attribute (intvec)"
+{
+  string t = typeof(N);
+  if( (t != "ideal") && (t != "module") && (t != "matrix") ) { return (0); };  // N should be something like a matrix
+
+  // N must be graded!
+  if ( typeof(attrib(N, "isHomog")) != "intvec" ){ return (0); };
+
+  intvec gr = attrib(N, "isHomog"); // grading weights...
+  if ( nrows(N) > size(gr) ) { return (0); };  // wrong number of rows?
+
+//  if( attrib(N, "rank") != size(gr) ){ return (0); } // wrong rank :(
+
+
+//  if( !homog(N) ) { return (0); };  // N should be homogenous
+
+  
+  return (1);
 }
 
 proc grisequal (def A, def B)
 "TODO"
-{
-  return (1==1); // TODO!
+{ 
+  ASSUME(1, grtest(A) ); 
+  ASSUME(1, grtest(B) );
+ 
+  int ra = nrows(A);
+  int rb = nrows(B);
+  
+  intvec wa = attrib(A, "isHomog");
+  intvec wb = attrib(B, "isHomog");
+  
+  if( (ra != rb) || (ncols(A) != ncols(B)) ){ return (0); } // TODO: ???
+  return ( (wa == wb) && 
+           (size(module(matrix(A) - matrix(B))) == 0)    ); // TODO: ???
 }
 
@@ -703 +636 @@
 "
 {
-  module Z; Z[a] = 0; 
-  Z = grobj(Z, intvec(d:a));
+  module Z; Z[a] = [0]; 
+  Z = grobj(Z, -intvec(d:a)); // will set the rank as well
 
   ASSUME(2, grisequal(Z, grpower( grshift(grzero(), d), a ) )); // optional check
@@ -710 +643 @@
 }
 
-proc grobj(module M, intvec w)
+proc grobj(def A, intvec w)
 ""
 {
+  module M = module(A);
+  ASSUME(0, size(w) >= nrows(M) );
+  attrib(M, "rank", size(w));
   attrib(M, "isHomog", w);
-  attrib(M, "rank", size(w));
+  ASSUME(1, grtest(M) );
   return (M);
 }
+
+
+static proc align( def A, int d)
+"analog of align kernel command for older Singular versions
+ this is static since it should not be used by @code{align}-able (newer)
+ Singular releases. 
+ Note that this proc does not care about any attributes (of A)
+"
+{
+  module T; T[d] = 0; 
+  T = T, module(transpose(A));
+  return( module(transpose(T)) );
+}

Singular/LIB/grobcov.lib

@@ -1 @@
-//////////////////////////////////////////////////////////////////////////////
-version="version grobcov.lib 4-0-0-0 Dec_2013 ";
+//
+version="version grobcov.lib 4.0.1.2 Jan_2015 "; // $Id$
 category="General purpose";
 info="
 LIBRARY:  grobcov.lib   Groebner Cover for parametric ideals.
-PURPOSE:  Comprehensive Groebner Systems, Groebner Cover, Canonical Forms,
-          Parametric Polynomial Systems.
-          The library contains Montes-Wibmer's algorithms to compute the
-          canonical Groebner cover of a parametric ideal as described in
-          the paper:
-
-          Montes A., Wibmer M.,
-          \"Groebner Bases for Polynomial Systems with parameters\".
-          Journal of Symbolic Computation 45 (2010) 1391-1425.
-          The locus algorithm and definitions will be
-          published in a forthcoming paper by Abanades, Botana, Montes, Recio:
-          \''An Algebraic Taxonomy for Locus Computation in Dynamic Geometry\''.
-
-          The central routine is grobcov. Given a parametric
-          ideal, grobcov outputs its Canonical Groebner Cover, consisting
-          of a set of pairs of (basis, segment). The basis (after
-          normalization) is the reduced Groebner basis for each point
-          of the segment. The segments are disjoint, locally closed
-          and correspond to constant lpp (leading power product)
-          of the basis, and are represented in canonical prime
-          representation. The segments are disjoint and cover the
-          whole parameter space. The output is canonical, it only
-          depends on the given parametric ideal and the monomial order.
-          This is much more than a simple Comprehensive Groebner System.
-          The algorithm grobcov allows options to solve partially the
-          problem when the whole automatic algorithm does not finish
-          in reasonable time.
-
-          grobcov uses a first algorithm cgsdr that outputs a disjoint
-          reduced Comprehensive Groebner System with constant lpp.
-          For this purpose, in this library, the implemented algorithm is
-          Kapur-Sun-Wang algorithm, because it is the most efficient
-          algorithm known for this purpose.
-
-          D. Kapur, Y. Sun, and D.K. Wang.
-          \"A New Algorithm for Computing Comprehensive Groebner Systems\".
-          Proceedings of ISSAC'2010, ACM Press, (2010), 29-36.
-
-          cgsdr can be called directly if only a disjoint reduced
-          Comprehensive Groebner System (CGS) is required.
-
-AUTHORS:  Antonio Montes , Hans Schoenemann.
-OVERVIEW: see \"Groebner Bases for Polynomial Systems with parameters\"
-          Montes A., Wibmer M.,
+
+          Comprehensive Groebner Systems, Groebner Cover, Canonical Forms,
+          Parametric Polynomial Systems, Dynamic Geometry, Loci, Envelop,
+          Constructible sets.
+          See
+
+          A. Montes A, M. Wibmer,
+          \"Groebner Bases for Polynomial Systems with parameters\",
           Journal of Symbolic Computation 45 (2010) 1391-1425.
           (http://www-ma2.upc.edu/~montes/).
+
+AUTHORS:  Antonio Montes , Hans Schoenemann.
+
+OVERVIEW:
+            In 2010, the library was designed to contain Montes-Wibmer's
+            algorithms for compute the canonical Groebner Cover of a
+            parametric ideal as described in the paper:
+
+            Montes A., Wibmer M.,
+            \"Groebner Bases for Polynomial Systems with parameters\".
+            Journal of Symbolic Computation 45 (2010) 1391-1425.
+
+            The central routine is grobcov. Given a parametric
+            ideal, grobcov outputs its Canonical Groebner Cover, consisting
+            of a set of pairs of (basis, segment). The basis (after
+            normalization) is the reduced Groebner basis for each point
+            of the segment. The segments are disjoint, locally closed
+            and correspond to constant lpp (leading power product)
+            of the basis, and are represented in canonical prime
+            representation. The segments are disjoint and cover the
+            whole parameter space. The output is canonical, it only
+            depends on the given parametric ideal and the monomial order.
+            This is much more than a simple Comprehensive Groebner System.
+            The algorithm grobcov allows options to solve partially the
+            problem when the whole automatic algorithm does not finish
+            in reasonable time.
+
+            grobcov uses a first algorithm cgsdr that outputs a disjoint
+            reduced Comprehensive Groebner System with constant lpp.
+            For this purpose, in this library, the implemented algorithm is
+            Kapur-Sun-Wang algorithm, because it is the most efficient
+            algorithm known for this purpose.
+
+            D. Kapur, Y. Sun, and D.K. Wang.
+            \"A New Algorithm for Computing Comprehensive Groebner Systems\".
+            Proceedings of ISSAC'2010, ACM Press, (2010), 29-36.
+
+            The library has evolved to include new applications of the
+            Groebner Cover, and new theoretical developments have been done.
+
+            A new set of routines (locus, locusdg, locusto) has been included to
+            compute loci of points. The routines are used in the Dynamic
+            Geometry software Geogebra. They are described in:
+
+            Abanades, Botana, Montes, Recio:
+            \''An Algebraic Taxonomy for Locus Computation in Dynamic Geometry\''.
+            Computer-Aided Design 56 (2014) 22-33.
+
+            Recently also routines for computing the envelop of a family
+            of curves (enverlop, envelopdg), to be used in Dynamic Geometry,
+            has been included and will be described in a forthcoming paper:
+
+             Abanades, Botana, Montes, Recio:
+             \''Envelops in Dynamic Geometry using the Groebner cover\''.
+
+            The actual version also includes two routines (AddCons and AddconsP)
+            for computing the canonical form of a constructible set, given as a
+            union of locally closed sets. They are used in the new version for the
+            computation of loci and envelops. It determines the canonical locally closed
+            level sets of a constructuble. They will be described in a forthcoming paper:
+
+            A. Montes, J.M. Brunat,
+            \"Canonical representations of constructible sets\".
+
+            This version was finished on 31/01/2015
 
 NOTATIONS: All given and determined polynomials and ideals are in the
@@ -62 +91 @@
 @*         If you want to use some internal routine you must
 @*         create before the above rings by calling setglobalrings();
-@*         because most of the internal routines use these rings.
+@*         because some of the internal routines use these rings.
 @*         The call to the basic routines grobcov, cgsdr will
 @*         kill these rings.
 
 PROCEDURES:
-
 grobcov(F);  Is the basic routine giving the canonical
-                   Groebner cover of the parametric ideal F.
-                   This routine accepts many options, that
-                   allow to obtain results even when the canonical
-                   computation does not finish in reasonable time.
-
-cgsdr(F);      Is the procedure for obtaining a first disjoint,
-                   reduced Comprehensive Groebner System that
-                   is used in grobcov, that can also be used
-                   independently if only the CGS is required.
-                   It is a more efficient routine than buildtree
-                   (the own routine that is no more used). It uses
-                   now KSW algorithm.
-
-setglobalrings();  Generates the global rings @R, @P and @PR that are
-                   respectively the rings Q[a][x], Q[a], Q[x,a].
-                   It is called inside each of the fundamental routines
-                   of the library: grobcov, cgsdr, locus, locusdg and killed before
-                   the output.
-                   So, if the user want to use some other internal routine,
-                   then setglobalrings() is to be called before, as most of them use
-                   these rings.
-
-pdivi(f,F);     Performs a pseudodivision of a parametric polynomial
-                   by a parametric ideal.
-                   Can be used without calling presiouly setglobalrings(),
+             Groebner Cover of the parametric ideal F.
+             This routine accepts many options, that
+             allow to obtain results even when the canonical
+             computation does not finish in reasonable time.
+
+cgsdr(F); Is the procedure for obtaining a first disjoint,
+             reduced Comprehensive Groebner System that
+             is used in grobcov, that can also be used
+             independently if only the CGS is required.
+             It is a more efficient routine than buildtree
+             (the own routine of 2010 that is no more used).
+             Now, KSW algorithm is used.
+
+setglobalrings();  Generates the global rings @R, @P and @PR that
+              are respectively the rings Q[a][x], Q[a], Q[x,a].  (a=parameters,
+             x=variables) It is called inside each of the fundamental
+             routines of the library: grobcov, cgsdr, locus, locusdg and
+             killed before the output.
+             If the user want to use some other internal routine,
+             then setglobalrings() is to be called before, as most of
+             them use these rings.
+
+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.
-                   Can be used without calling presiouly setglobalrings(),
-
-Prep(N,M);   Computes the P-representation of V(N) \ V(M).
-                   Can be used without calling previously setglobalrings().
-
-extend(GC); When the grobcov of an ideal has been computed
-                   with the default option ('ext',0) and the explicit
-                   option ('rep',2) (which is not the default), then
-                   one can call extend (GC) (and options) to obtain the
-                   full representation of the bases. With the default
-                   option ('ext',0) only the generic representation of
-                   the bases are computed, and one can obtain the full
-                   representation using extend.
-                   Can be used without calling presiouly setglobalrings(),
-
-locus(G):     Special routine for determining the locus of points
-                   of  objects. Given a parametric ideal J with
-                   parameters (a_1,..a_m) and variables (x_1,..,xn),
-                   representing the system determining
-                   the locus of points (a_1,..,a_m)) who verify certain
-                   properties, computing the grobcov G of
-                   J and applying to it locus, determines the different
-                   classes of locus components. They can be
-                   Normal, Special, Accumulation point, Degenerate.
-                   The output are the components given in P-canonical form
-                   of two constructible sets: Normal, and Non-Normal
-                   The Normal Set has Normal and Special components
-                   The Non-Normal set has Accumulation and Degenerate components.
-                   The description of the algorithm and definitions will be
-                   given in a forthcoming paper by Abanades, Botana, Montes, Recio:
-                   ''An Algebraic Taxonomy for Locus Computation in Dynamic Geometry''.
-                   Can be used without calling presiouly setglobalrings(),
-
-locusdg(G):  Is a special routine for computing the locus in dinamical geometry.
-                    It detects automatically a possible point that is to be avoided by the mover,
-                    whose coordinates must be the last two coordinates in the definition of the ring.
-                    If such a point is detected, then it eliminates the segments of the grobcov
-                    depending on the point that is to be avoided.
-                    Then it calls locus.
-                    Can be used without calling presiouly setglobalrings(),
-
-locusto(L):  Transforms the output of locus to a string that
-                  can be reed from different computational systems.
-                  Can be used without calling presiouly setglobalrings(),
-
-addcons(L): Given a disjoint set of locally closed subsets in P-representation,
-                  it returns the canonical P-representation of the constructible set
-                  formed by the union of them.
-                  Can be used without calling presiouly setglobalrings(),
+             ( E is the null ideal and N the non-null ideal )
+             over the parameters.
+
+Crep(N,M);  Computes the canonical C-representation of V(N) \ V(M).
+
+Prep(N,M);  Computes the canonical P-representation of V(N) \ V(M).
+
+PtoCrep(L)  Starting from the canonical Prep of a locally closed set
+             computes its Crep.
+
+extend(GC);  When the grobcov of an ideal has been computed
+             with the default option ('ext',0) and the explicit
+             option ('rep',2) (which is not the default), then
+             one can call extend (GC) (and options) to obtain the
+             full representation of the bases. With the default
+             option ('ext',0) only the generic representation of
+             the bases are computed, and one can obtain the full
+             representation using extend.
+
+locus(G);    Special routine for determining the geometrical locus of points
+             verifying given conditions. Given a parametric ideal J with
+             parameters (x,y) and variables (x_1,..,xn), representing the
+             system determining the locus of points (x,y) who verify certain
+             properties, one can apply locus to the output of  grobcov(J),
+             locus determines the different classes of locus components.
+              described in the paper:
+
+             \"An Algebraic Taxonomy for Locus Computation in Dynamic Geometry\",
+             M. Abanades, F. Botana, A. Montes, T. Recio,
+             Computer-Aided Design 56 (2014) 22-33.
+
+             The components can be 'Normal', 'Special', 'Accumulation', 'Degenerate'.
+             The output are the components is given in P-canonical form
+             It also detects automatically a possible point that is to be
+             avoided by the mover, whose coordinates must be the last two
+             coordinates in the definition of the ring. If such a point is detected,
+             then it eliminates the segments of the grobcov depending on the
+             point that is to be avoided.
+
+locusdg(G);  Is a special routine that determines the  'Relevant' components
+             of the locus in dynamic geometry. It is to be called to the output of locus
+             and selects from it the useful components.
+
+envelop(F,C); Special routine for determining the envelop of a family of curves
+             F in Q[x,y][x_1,..xn] depending on a ideal of constraints C in Q[x_1,..,xn].
+             It detemines the different components as well as its type:
+             'Normal', 'Special', 'Accumulation', 'Degenerate'. And
+             it also classifies the Special components, determining the
+             zero dimensional antiimage of the component and verifying if
+             the component is a special curve of the family or not.
+             It calls internally first grobcov and then locus with special options
+             to obtain the complete result.
+             The taxonomy that it provides, as well as the algorithms involved
+             will be described in a forthcoming paper:
+
+             Abanades, Botana, Montes, Recio:
+             \''Envelops in Dynamic Geometry using the Gr\"obner cover\''.
+
+
+envelopdg(ev); Is a special routine to determine the 'Relevant' components
+             of the envelop of a family of curves to be used in Dynamic Geometry.
+             It must be called to the output of envelop(F,C).
+
+locusto(L); Transforms the output of locus, locusdg, envelop and  envelopdg
+             into a string that can be reed from different computational systems.
+
+AddCons(L); Uses the routine AddConsP. Given a set of locally closed sets as
+             difference of of varieties (does not need to be in C-representation)
+             it returns the canonical P-representation of the constructible set
+             formed by the union of them. The result is formed by a set of embedded,
+             disjoint, locally closed sets (levels).
+
+AddConsP(L);  Given a set of locally closed P-components, it returns the
+             canonical P-representation of the constructible set
+             formed by the union of them, consisting of a set of embedded,
+             disjoint locally closed sets (levels).
+             The routines AddCons and AddConsP and the canonical structure
+             of the constructible sets will be described in a forthcoming paper.
+
+             A. Montes, J.M. Brunat,
+             \"Canonical representations of constructible sets\".
 
 SEE ALSO: compregb_lib
@@ -157 +214 @@
 
 // Library grobcov.lib
-// (Groebner cover):
-// Release 1: (public)
+// (Groebner Cover):
+// Release 0: (public)
 // Initial data: 21-1-2008
+// Uses buildtree for cgsdr
 // Final data: 3-7-2008
 // Release 2: (private)
 // Initial data: 6-9-2009
+// Last version using buildtree for cgsdr
 // Final data: 25-10-2011
-// Release 3:
+// Release B: (private)
 // Initial data: 1-7-2012
+// Uses both buildtree and KSW for cgsdr
 // Final data: 4-9-2012
-// Release 4: (this release, public)
+// Release G: (public)
 // Initial data: 4-9-2012
+// Uses KSW algorithm for cgsdr
 // Final data: 21-11-2013
-// basering Q[a][x];
-
-// ************ Begin of buildtree ******************************
-
-proc setglobalrings()
-"USAGE:   setglobalrings();
-          No arguments
-RETURN:   After its call the rings @R=Q[a][x], @P=Q[a], @RP=Q[x,a] are
-          defined as global variables.
-NOTE: It is called internally by the fundamental routines of the
-          library grobcov, cgsdr, extend, pdivi, pnormalf, locus, locusdg, locusto,
-          and killed before the output.
-          The user does not need to call it, except when it is interested
-          in using some internal routine of the library that
-          uses these rings.
-          The basering R, must be of the form Q[a][x], a=parameters,
-          x=variables, and should be defined previously.
-KEYWORDS: ring, rings
-EXAMPLE:  setglobalrings; shows an example"
-{
-  if (defined(@P))
-  {
-    kill @P; kill @R; kill @RP;
-  }
+// Release K: (public)
+// Updated locus: 8-1-2015
+// Updated AddConsP and AddCons: 8-1-2015
+// Reformed many routines, examples and helps: 8-1-2015
+// New routines for computing the envelop of a family of curves: 22-1-2015
+// Final data: 22-1-2015
+
+//*************Auxiliary routines**************
+
+// elimintfromideal: elimine the constant numbers from an ideal
+//        (designed for W, nonnull conditions)
+// Input: ideal J
+// Output:ideal K with the elements of J that are non constants, in the
+//        ring K[x1,..,xm]
+static proc elimintfromideal(ideal J)
+{
+  int i;
+  int j=0;
+  ideal K;
+  if (size(J)==0){return(ideal(0));}
+  for (i=1;i<=ncols(J);i++){if (size(variables(J[i])) !=0){j++; K[j]=J[i];}}
+  return(K);
+}
+
+// delfromideal: deletes the i-th polynomial from the ideal F
+//    Works in any kind of ideal
+static proc delfromideal(ideal F, int i)
+{
+  int j;
+  ideal G;
+  if (size(F)<i){ERROR("delfromideal was called with incorrect arguments");}
+  if (size(F)<=1){return(ideal(0));}
+  if (i==0){return(F)};
+  for (j=1;j<=ncols(F);j++)
+  {
+    if (j!=i){G[ncols(G)+1]=F[j];}
+  }
+  return(G);
+}
+
+// delidfromid: deletes the polynomials in J that are in I
+// Input: ideals I, J
+// Output: the ideal J without the polynomials in I
+//   Works in any kind of ideal
+static proc delidfromid(ideal I, ideal J)
+{
+  int i; list r;
+  ideal JJ=J;
+  for (i=1;i<=size(I);i++)
+  {
+    r=memberpos(I[i],JJ);
+    if (r[1])
+    {
+      JJ=delfromideal(JJ,r[2]);
+    }
+  }
+  return(JJ);
+}
+
+// eliminates the ith element from a list or an intvec
+static proc elimfromlist(l, int i)
+{
+  if(typeof(l)=="list"){list L;}
+  if (typeof(l)=="intvec"){intvec L;}
+  if (typeof(l)=="ideal"){ideal L;}
+  int j;
+  if((size(l)==0) or (size(l)==1 and i!=1)){return(l);}
+  if (size(l)==1 and i==1){return(L);}
+  // L=l[1];
+  if(i==1)
+  {
+    for(j=2;j<=size(l);j++)
+    {
+      L[j-1]=l[j];
+    }
+  }
+  else
+  {
+    for(j=1;j<=i-1;j++)
+    {
+      L[j]=l[j];
+    }
+    for(j=i+1;j<=size(l);j++)
+    {
+      L[j-1]=l[j];
+    }
+  }
+  return(L);
+}
+
+// eliminates repeated elements form an ideal or matrix or module or intmat or bigintmat
+static proc elimrepeated(F)
+{
+  int i;
+  int nt;
+  if (typeof(F)=="ideal"){nt=ncols(F);}
+  else{nt=size(F);}
+
+  def FF=F;
+  FF=F[1];
+  for (i=2;i<=nt;i++)
+  {
+    if (not(memberpos(F[i],FF)[1]))
+    {
+      FF[size(FF)+1]=F[i];
+    }
+  }
+  return(FF);
+}
+
+
+// equalideals
+// Input: ideals F and G;
+// Output: 1 if they are identical (the same polynomials in the same order)
+//         0 else
+static proc equalideals(ideal F, ideal G)
+{
+  int i=1; int t=1;
+  if (size(F)!=size(G)){return(0);}
+  while ((i<=size(F)) and (t))
+  {
+    if (F[i]!=G[i]){t=0;}
+    i++;
+  }
+  return(t);
+}
+
+// returns 1 if the two lists of ideals are equal and 0 if not
+static proc equallistideals(list L, list M)
+{
+  int t; int i;
+  if (size(L)!=size(M)){return(0);}
+  else
+  {
+    t=1;
+    if (size(L)>0)
+    {
+      i=1;
+      while ((t) and (i<=size(L)))
+      {
+        if (equalideals(L[i],M[i])==0){t=0;}
+        i++;
+      }
+    }
+    return(t);
+  }
+}
+
+// idcontains
+// Input: ideal p, ideal q
+// Output: 1 if p contains q,  0 otherwise
+// If the routine is to be called from the top, a previous call to
+// setglobalrings() is needed.
+static proc idcontains(ideal p, ideal q)
+{
+  int t; int i;
+  t=1; i=1;
+  def P=p; def Q=q;
+  attrib(P,"isSB",1);
+  poly r;
+  while ((t) and (i<=size(Q)))
+  {
+    r=reduce(Q[i],P);
+    if (r!=0){t=0;}
+    i++;
+  }
+  return(t);
+}
+
+
+// selectminideals
+//   given a list of ideals returns the list of integers corresponding
+//   to the minimal ideals in the list
+// Input: L (list of ideals)
+// Output: the list of integers corresponding to the minimal ideals in L
+//   Works in Q[u_1,..,u_m]
+static proc selectminideals(list L)
+{
+  list P; int i; int j; int t;
+  if(size(L)==0){return(L)};
+  if(size(L)==1){P[1]=1; return(P);}
+  for (i=1;i<=size(L);i++)
+  {
+    t=1;
+    j=1;
+    while ((t) and (j<=size(L)))
+    {
+      if (i!=j)
+      {
+        if(idcontains(L[i],L[j])==1)
+        {
+          t=0;
+        }
+      }
+      j++;
+    }
+    if (t){P[size(P)+1]=i;}
+  }
+  return(P);
+}
+
+
+// Auxiliary routine
+// elimconstfac: eliminate the factors in the polynom f that are in Q[a]
+// Input:
+//   poly f:
+//   list L: of components of the segment
+// Output:
+//   poly f2  where the factors of f in Q[a] that are non-null on any component
+//   have been dropped from f
+static proc elimconstfac(poly f)
+{
+  int cond; int i; int j; int t;
+  if (f==0){return(f);}
   def RR=basering;
-  def @R=basering;  // must be of the form K[a][x], a=parameters, x=variables
-  list Rx=ringlist(RR);
-  def @P=ring(Rx[1]);
-  list Lx;
-  Lx[1]=0;
-  Lx[2]=Rx[2]+Rx[1][2];
-  Lx[3]=Rx[1][3];
-  Lx[4]=Rx[1][4];
-  Rx[1]=0;
-  def D=ring(Rx);
-  def @RP=D+@P;
-  export(@R);      // global ring K[a][x]
-  export(@P);      // global ring K[a]
-  export(@RP);     // global ring K[x,a] with product order
+  setring(@R);
+  def ff=imap(RR,f);
+  def l=factorize(ff,0);
+  poly f1=1;
+  for(i=2;i<=size(l[1]);i++)
+  {
+      f1=f1*(l[1][i])^(l[2][i]);
+  }
   setring(RR);
-}
-example
-{ "EXAMPLE:"; echo = 2;
-  ring R=(0,a,b),(x,y,z),dp;
-  setglobalrings();
-  Grobcov::@R;
-  Grobcov::@P;
-  Grobcov::@RP;
-ringlist(Grobcov::@R);
-ringlist(Grobcov::@P);
-ringlist(Grobcov::@RP);
-}
-
-//*************Auxiliary routines**************
-
-// cld : clears denominators of an ideal and normalizes to content 1
-//        can be used in @R or @P or @RP
-// input:
-//        ideal J (J can be also poly), but the output is an ideal;
-// output:
-//        ideal Jc (the new form of ideal J without denominators and
-//        normalized to content 1)
-proc cld(ideal J)
-{
-  if (size(J)==0){return(ideal(0));}
-  int te=0;
-  def RR=basering;
-  if(not(defined(@RP)))
-  {
-    te=1;
-    setglobalrings();
-  }
-  setring(@RP);
-  def Ja=imap(RR,J);
-  ideal Jb;
-  if (size(Ja)==0){setring(RR); return(ideal(0));}
-  int i;
-  def j=0;
-  for (i=1;i<=ncols(Ja);i++){if (size(Ja[i])!=0){j++; Jb[j]=cleardenom(Ja[i]);}}
-  setring(RR);
-  def Jc=imap(@RP,Jb);
-  if(te){kill @R; kill @RP; kill @P;}
-  return(Jc);
-}
-
-proc memberpos(f,J)
+  def f2=imap(@R,f1);
+  return(f2);
+};
+
+static proc memberpos(f,J)
 //"USAGE:  memberpos(f,J);
 //         (f,J) expected (polynomial,ideal)
@@ -264 +468 @@
 //               or       (ideal,list(ideal))
 //               or       (list(intvec),  list(list(intvec))).
-//         The ring can be @R or @P or @RP or any other.
+//         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.
@@ -379 +583 @@
 //}
 
-proc subset(J,K)
+// Auxiliary routine
+// pos
+// Input:  intvec p of zeros and ones
+// Output: intvec W of the positions where p has ones.
+static proc pos(intvec p)
+{
+  int i;
+  intvec W; int j=1;
+  for (i=1; i<=size(p); i++)
+  {
+    if (p[i]==1){W[j]=i; j++;}
+  }
+  return(W);
+}
+
+// Input:
+//  A,B: lists of ideals
+// Output:
+//   1 if both lists of ideals are equal, or 0 if not
+static proc equallistsofideals(list A, list B)
+{
+ int i;
+ int tes=0;
+ if (size(A)!=size(B)){return(tes);}
+ tes=1; i=1;
+ while(tes==1 and i<=size(A))
+ {
+   if (equalideals(A[i],B[i])==0){tes=0; return(tes);}
+   i++;
+ }
+ return(tes);
+}
+
+// Input:
+//  A,B:  lists of P-rep, i.e. of the form [p_i,[p_{i1},..,p_{ij_i}]]
+// Output:
+//   1 if both lists of P-reps are equal, or 0 if not
+static proc equallistsA(list A, list B)
+{
+  int tes=0;
+  if(equalideals(A[1],B[1])==0){return(tes);}
+  tes=equallistsofideals(A[2],B[2]);
+  return(tes);
+}
+
+// Input:
+//  A,B:  lists lists of of P-rep, i.e. of the form [[p_1,[p_{11},..,p_{1j_1}]] .. [p_i,[p_{i1},..,p_{ij_i}]]
+// Output:
+//   1 if both lists of lists of P-rep are equal, or 0 if not
+static proc equallistsAall(list A,list B)
+{
+ int i; int tes;
+ if(size(A)!=size(B)){return(tes);}
+ tes=1; i=1;
+ while(tes and i<=size(A))
+ {
+   tes=equallistsA(A[i],B[i]);
+   i++;
+ }
+ return(tes);
+}
+
+// idint: ideal intersection
+//        in the ring @P.
+//        it works in an extended ring
+// input: two ideals in the ring @P
+// output the intersection of both (is not a GB)
+static proc idint(ideal I, ideal J)
+{
+  def RR=basering;
+  ring T=0,t,lp;
+  def K=T+RR;
+  setring(K);
+  def Ia=imap(RR,I);
+  def Ja=imap(RR,J);
+  ideal IJ;
+  int i;
+  for(i=1;i<=size(Ia);i++){IJ[i]=t*Ia[i];}
+  for(i=1;i<=size(Ja);i++){IJ[size(Ia)+i]=(1-t)*Ja[i];}
+  ideal eIJ=eliminate(IJ,t);
+  setring(RR);
+  return(imap(K,eIJ));
+}
+
+//purpose ideal intersection called in @R and computed in @P
+static proc idintR(ideal N, ideal M)
+{
+  def RR=basering;
+  setring(@P);
+  def Np=imap(RR,N);
+  def Mp=imap(RR,M);
+  def Jp=idint(Np,Mp);
+  setring(RR);
+  return(imap(@P,Jp));
+}
+
+// Auxiliary routine
+// comb: the list of combinations of elements (1,..n) of order p
+static proc comb(int n, int p)
+{
+  list L; list L0;
+  intvec c; intvec d;
+  int i; int j; int last;
+  if ((n<0) or (n<p))
+  {
+    return(L);
+  }
+  if (p==1)
+  {
+    for (i=1;i<=n;i++)
+    {
+      c=i;
+      L[size(L)+1]=c;
+    }
+    return(L);
+  }
+  else
+  {
+    L0=comb(n,p-1);
+    for (i=1;i<=size(L0);i++)
+    {
+      c=L0[i]; d=c;
+      last=c[size(c)];
+      for (j=last+1;j<=n;j++)
+      {
+        d[size(c)+1]=j;
+        L[size(L)+1]=d;
+      }
+    }
+    return(L);
+  }
+}
+
+// Auxiliary routine
+// combrep
+// Input: V=(n_1,..,n_i)
+// Output: L=(v_1,..,v_p) where p=prod_j=1^i (n_j)
+//    is the list of all intvec v_j=(v_j1,..,v_ji) where 1<=v_jk<=n_i
+static proc combrep(intvec V)
+{
+  list L; list LL;
+  int i; int j; int k;  intvec W;
+  if (size(V)==1)
+  {
+    for (i=1;i<=V[1];i++)
+    {
+      L[i]=intvec(i);
+    }
+    return(L);
+  }
+  for (i=1;i<size(V);i++)
+  {
+    W[i]=V[i];
+  }
+  LL=combrep(W);
+  for (i=1;i<=size(LL);i++)
+  {
+    W=LL[i];
+    for (j=1;j<=V[size(V)];j++)
+    {
+      W[size(V)]=j;
+      L[size(L)+1]=W;
+    }
+  }
+  return(L);
+}
+
+static proc subset(J,K)
 //"USAGE:   subset(J,K);
 //          (J,K)  expected (ideal,ideal)
@@ -405 +776 @@
 //}
 
-// elimintfromideal: elimine the constant numbers from an ideal
-//        (designed for W, nonnull conditions)
-// input: ideal J
-// output:ideal K with the elements of J that are non constants, in the
-//        ring @P
-proc elimintfromideal(ideal J)
-{
+proc setglobalrings()
+"USAGE:   setglobalrings();
+          No arguments
+RETURN:   After its call the rings @R=Q[a][x], @P=Q[a], @RP=Q[x,a] are
+          defined as global variables.  (a=parameters, x=variables)
+NOTE: It is called internally by many basic routines of the
+          library grobcov, cgsdr, extend, pdivi, pnormalf, locus, locusdg,
+          envelop, envelopdg, and killed before the output.
+          The user does not need to call it, except when it is interested
+          in using some internal routine of the library that
+          uses these rings.
+          The basering R, must be of the form Q[a][x], (a=parameters,
+          x=variables), and should be defined previously.
+KEYWORDS: ring, rings
+EXAMPLE:  setglobalrings; shows an example"
+{
+  if (defined(@P))
+  {
+    kill @P; kill @R; kill @RP;
+  }
+  def RR=basering;
+  def @R=basering;  // must be of the form Q[a][x], (a=parameters, x=variables)
+  def Rx=ringlist(RR);
+  def @P=ring(Rx[1]);
+  list Lx;
+  Lx[1]=0;
+  Lx[2]=Rx[2]+Rx[1][2];
+  Lx[3]=Rx[1][3];
+  Lx[4]=Rx[1][4];
+  Rx[1]=0;
+  def D=ring(Rx);
+  def @RP=D+@P;
+  export(@R);      // global ring Q[a][x]
+  export(@P);      // global ring Q[a]
+  export(@RP);     // global ring K[x,a] with product order
+  setring(RR);
+};
+example
+{
+  "EXAMPLE:"; echo = 2;
+  ring R=(0,a,b),(x,y,z),dp;
+  setglobalrings();
+  R;
+  Grobcov::@R;
+  Grobcov::@P;
+  Grobcov::@RP;
+  ringlist(Grobcov::@R);
+  ringlist(Grobcov::@P);
+  ringlist(Grobcov::@RP);
+}
+
+// cld : clears denominators of an ideal and normalizes to content 1
+//        can be used in @R or @P or @RP
+// input:
+//        ideal J (J can be also poly), but the output is an ideal;
+// output:
+//        ideal Jc (the new form of ideal J without denominators and
+//        normalized to content 1)
+static proc cld(ideal J)
+{
+  if (size(J)==0){return(ideal(0));}
+  int te=0;
+  def RR=basering;
+  if(not(defined(@RP)))
+  {
+    te=1;
+    setglobalrings();
+  }
+  setring(@RP);
+  def Ja=imap(RR,J);
+  ideal Jb;
+  if (size(Ja)==0){setring(RR); return(ideal(0));}
   int i;
-  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);
-}
+  def j=0;
+  for (i=1;i<=ncols(Ja);i++){if (size(Ja[i])!=0){j++; Jb[j]=cleardenom(Ja[i]);}}
+  setring(RR);
+  def Jc=imap(@RP,Jb);
+  if(te){kill @R; kill @RP; kill @P;}
+  return(Jc);
+};
 
 // simpqcoeffs : simplifies a quotient of two polynomials
 // input: two coeficients (or terms), that are considered as a quotient
 // output: the two coeficients reduced without common factors
-proc simpqcoeffs(poly n,poly m)
-{
-  number nc=content(n);
-  number mc=content(m);
-  number gc=gcd(nc,mc);
+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);
@@ -439 +876 @@
 //   list (poly r, ideal q, poly mu)
 proc pdivi(poly f,ideal F)
-"USAGE:   pdivi(f,F);
-          poly f: the polynomial to be divided
-          ideal F: the divisor ideal
-RETURN:   A list (poly r, ideal q, poly m). r is the remainder of the
+"USAGE: pdivi(f,F);
+          poly f: the polynomialin Q[a][x] to be divided
+          ideal F: the divisor ideal in Q[a][x].
+RETURN: A list (poly r, ideal q, poly m). r is the remainder of the
           pseudodivision, q is the set of quotients, and m is the
           coefficient factor by which f is to be multiplied.
-NOTE:     pseudodivision of a poly f by an ideal F in @R. Returns a
+NOTE: pseudodivision of a poly f by an ideal F in Q[a][x]. Returns a
           list (r,q,m) such that m*f=r+sum(q.G), and no lpp of a divisor
           divides a pp of r.
@@ -451 +888 @@
 EXAMPLE:  pdivi; shows an example"
 {
+  F=simplify(F,2);
   int i;
   int j;
-//  string("T_f=",f);
   poly v=1;
   for(i=1;i<=nvars(basering);i++){v=v*var(i);}
-  int te=0;
-  def R=basering;
-  if (defined(@P)==1){te=1;}
-  else{setglobalrings();}
   poly r=0;
   poly mu=1;
   def p=f;
   ideal q;
-  for (i=1; i<=size(F); i++){q[i]=0;};
+  for (i=1; i<=ncols(F); i++){q[i]=0;};
   ideal lpf;
   ideal lcf;
-  for (i=1;i<=size(F);i++){lpf[i]=leadmonom(F[i]);}
-  for (i=1;i<=size(F);i++){lcf[i]=leadcoef(F[i]);}
+  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;
@@ -492 +925 @@
         rho=qlc[1]*qlm;
         p=nu*p-rho*F[i];
-//        string("i=",i,"  coef(p,v)=",coef(p,v));
         r=nu*r;
         for (j=1;j<=size(F);j++){q[j]=nu*q[j];}
@@ -504 +936 @@
       r=r+lcp*lpp;
       p=p-lcp*lpp;
-//      string("pasa al rem p=",p);
     }
   }
   list res=r,q,mu;
-  if(te==0){kill @P; kill @R; kill @RP;}
   return(res);
 }
 example
-{ "EXAMPLE:"; echo = 2;
+{
+  "EXAMPLE:"; echo = 2;
   ring R=(0,a,b,c),(x,y),dp;
-  "Divisor=";
   poly f=(ab-ac)*xy+(ab)*x+(5c);
-  "Dividends=";
+  // Divisor=";
+  f;
   ideal F=ax+b,cy+a;
-  "(Remainder, quotients, factor)=";
+  // Dividends=";
+  F;
   def r=pdivi(f,F);
+  // (Remainder, quotients, factor)=";
   r;
-  "Verifying the division: r[3]*f-(r[2][1]*F[1]+r[2][2]*F[2])-r[1] =";
-  r[3]*f-(r[2][1]*F[1]+r[2][2]*F[2])-r[1];
+  // Verifying the division:
+  r[3]*f-(r[2][1]*F[1]+r[2][2]*F[2]+r[1]);
 }
 
@@ -534 +967 @@
 //                if red then S reduces by Buchberger 1st criterion
 //                (not used)
-proc pspol(poly f,poly g)
+static proc pspol(poly f,poly g)
 {
   def lcf=leadcoef(f);
@@ -559 +992 @@
 //         of ideal J (non repeated). Integer factors are ignored,
 //         even 0 is ignored. It can be called from ideal @R.
-proc facvar(ideal J)
+static proc facvar(ideal J)
 //"USAGE:   facvar(J);
 //          J: an ideal in the parameters
@@ -605 +1038 @@
 //   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)
+//        (E,N) represents V(E) \ V(N),
+//        Ered eliminates the non-null factors of f in V(E) \ V(N)
 // output:
 //   poly f2  where the non-null conditions have been dropped from f
-proc Ered(poly f,ideal E, ideal N)
+static proc Ered(poly f,ideal E, ideal N)
 {
   def RR=basering;
@@ -663 +1096 @@
 }
 
-// pnormalf: reduces a polynomial f wrt a V(E)\V(N)
+// pnormalf: reduces a polynomial f wrt a V(E) \ V(N)
 //           dividing by E and eliminating factors in N.
 //           called in the ring @R,
@@ -671 +1104 @@
 //         ideal E  (depends only on the parameters)
 //         ideal N  (depends only on the parameters)
-//                  (E,N) represents V(E)\V(N)
+//                  (E,N) represents V(E) \ V(N)
 //         optional:
-// output: poly f2 reduced wrt to V(E)\V(N)
+// output: poly f2 reduced wrt to V(E) \ V(N)
 proc pnormalf(poly f, ideal E, ideal N)
-"USAGE:   pnormalf(f,E,N);
-          f: the polynomial to be reduced modulo V(E)\V(N)
-          of a segment in the parameters.
-          E: the null conditions ideal
-          N: the non-null conditions
-RETURN:   a reduced polynomial g of f, whose coefficients are reduced
+"USAGE: pnormalf(f,E,N);
+          f: the polynomial in Q[a][x]  (a=parameters, x=variables) to be
+          reduced modulo V(E) \ V(N) of a segment in Q[a].
+          E: the null conditions ideal in Q[a]
+          N: the non-null conditions in Q[a]
+RETURN: a reduced polynomial g of f, whose coefficients are reduced
           modulo E and having no factor in N.
 NOTE: Should be called from ring Q[a][x].
-          Ideals E and N must be given by polynomials
-          in the parameters.
+          Ideals E and N must be given by polynomials in Q[a].
 KEYWORDS: division, pdivi, reduce
-EXAMPLE:  pnormalf; shows an example"
+EXAMPLE: pnormalf; shows an example"
 {
     def RR=basering;
@@ -705 +1137 @@
     if(te==0){kill @R; kill @RP; kill @P;}
     return(f2)
-}
+};
 example
-{ "EXAMPLE:"; echo = 2;
+{
+  "EXAMPLE:"; echo = 2;
   ring R=(0,a,b,c),(x,y),dp;
+  short=0;
   poly f=(b^2-1)*x^3*y+(c^2-1)*x*y^2+(c^2*b-b)*x+(a-bc)*y;
   ideal E=(c-1);
   ideal N=a-b;
+
   pnormalf(f,E,N);
-}
-
-// idint: ideal intersection
-//        in the ring @P.
-//        it works in an extended ring
-// input: two ideals in the ring @P
-// output the intersection of both (is not a GB)
-proc idint(ideal I, ideal J)
-{
-  def RR=basering;
-  ring T=0,t,lp;
-  def K=T+RR;
-  setring(K);
-  def Ia=imap(RR,I);
-  def Ja=imap(RR,J);
-  ideal IJ;
-  int i;
-  for(i=1;i<=size(Ia);i++){IJ[i]=t*Ia[i];}
-  for(i=1;i<=size(Ja);i++){IJ[size(Ia)+i]=(1-t)*Ja[i];}
-  ideal eIJ=eliminate(IJ,t);
-  setring(RR);
-  return(imap(K,eIJ));
 }
 
@@ -740 +1153 @@
 // input:  two polynomials f,g in the ring @R
 // output: 0 if f<g,  1 if f>=g
-proc lesspol(poly f, poly g)
+static proc lesspol(poly f, poly g)
 {
   if (leadmonom(f)<leadmonom(g)){return(1);}
@@ -752 +1165 @@
     }
   }
-}
-
-// delfromideal: deletes the i-th polynomial from the ideal F
-proc delfromideal(ideal F, int i)
-{
-  int j;
-  ideal G;
-  if (size(F)<i){ERROR("delfromideal was called incorrect arguments");}
-  if (size(F)<=1){return(ideal(0));}
-  if (i==0){return(F)};
-  for (j=1;j<=size(F);j++)
-  {
-    if (j!=i){G[size(G)+1]=F[j];}
-  }
-  return(G);
-}
-
-// delidfromid: deletes the polynomials in J that are in I
-// input: ideals I,J
-// output: the ideal J without the polynomials in I
-proc delidfromid(ideal I, ideal J)
-{
-  int i; list r;
-  ideal JJ=J;
-  for (i=1;i<=size(I);i++)
-  {
-    r=memberpos(I[i],JJ);
-    if (r[1])
-    {
-      JJ=delfromideal(JJ,r[2]);
-    }
-  }
-  return(JJ);
-}
+};
 
 // sortideal: sorts the polynomials in an ideal by lm in ascending order
-proc sortideal(ideal Fi)
+static proc sortideal(ideal Fi)
 {
   def RR=basering;
@@ -802 +1182 @@
     j=1;
     p=H[1];
-    for (i=1;i<=size(H);i++)
+    for (i=1;i<=ncols(H);i++)
     {
       if(lesspol(H[i],p)){j=i;p=H[j];}
     }
-    G[size(G)+1]=p;
+    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 is ascending order and
+// order the polynomials in ascending order and
 // eliminates the polynomials whose lpp are divisible by some
 // smaller one
-proc mingb(ideal F)
+static proc mingb(ideal F)
 {
   int t; int i; int j;
@@ -841 +1223 @@
 // redgbn: given a minimal basis (gb reducing) it
 // reduces each polynomial wrt to V(E) \ V(N)
-proc redgbn(ideal F, ideal E, ideal N)
+static proc redgbn(ideal F, ideal E, ideal N)
 {
   int te=0;
@@ -859 +1241 @@
 }
 
-// eliminates repeated elements form an ideal or matrix or module or intmat or bigintmat
-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);
-}
-
-proc elimrepeatedvl(F)
-{
-  int i;
-  def FF=F;
-  FF=F[1];
-  for (i=2;i<=size(F);i++)
-  {
-    if (not(memberpos(F[i],FF)[1]))
-    {
-      FF[size(FF)+1]=F[i];
-    }
-  }
-  return(FF);
-}
-
-
-// equalideals
-// input: 2 ideals F and G;
-// output: 1 if they are identical (the same polynomials in the same order)
-//         0 else
-proc equalideals(ideal F, ideal G)
-{
-  int i=1; int t=1;
-  if (size(F)!=size(G)){return(0);}
-  while ((i<=size(F)) and (t))
-  {
-    if (F[i]!=G[i]){t=0;}
-    i++;
-  }
-  return(t);
-}
-
-// delintvec
-// input: intvec V
-//        int i
-// output:
-//        intvec W (equal to V but the coordinate i is deleted
-proc delintvec(intvec V, int i)
-{
-  int j;
-  intvec W;
-  for (j=1;j<i;j++){W[j]=V[j];}
-  for (j=i+1;j<=size(V);j++){W[j-1]=V[j];}
-  return(W);
-}
-
 //**************Begin homogenizing************************
 
@@ -931 +1247 @@
 // input:  poly f
 // output  1 if f is homogeneous, 0 if not
-proc ishomog(f)
+static proc ishomog(f)
 {
   int i; poly r; int d; int dr;
@@ -947 +1263 @@
 // postredgb: given a minimal basis (gb reducing) it
 // reduces each polynomial wrt to the others
-proc postredgb(ideal F)
+static proc postredgb(ideal F)
 {
   int te=0;
@@ -964 +1280 @@
 }
 
-//purpose ideal intersection called in @R and computed in @P
-proc idintR(ideal N, ideal M)
-{
-  def RR=basering;
-  setring(@P);
-  def Np=imap(RR,N);
-  def Mp=imap(RR,M);
-  def Jp=idint(Np,Mp);
-  setring(RR);
-  return(imap(@P,Jp));
-}
 
 //purpose reduced Groebner basis called in @R and computed in @P
-proc gbR(ideal N)
+static proc gbR(ideal N)
 {
   def RR=basering;
@@ -997 +1302 @@
 //        poly f:
 // Output: Na = N:<f>
-proc incquotient(ideal N, poly f)
+static proc incquotient(ideal N, poly f)
 {
   poly g; int i;
@@ -1044 +1349 @@
 }
 
-// eliminates the ith element from a list or an intvec
-proc elimfromlist(def 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=2;j<=size(l);j++)
-    {
-      L[j-1]=l[j];
-    }
-  }
-  else
-  {
-    for(j=1;j<=i-1;j++)
-    {
-      L[j]=l[j];
-    }
-    for(j=i+1;j<=size(l);j++)
-    {
-      L[j-1]=l[j];
-    }
-  }
-  return(L);
-}
-
 // 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
-proc idbefid(ideal a, ideal b)
+static proc idbefid(ideal a, ideal b)
 {
   poly fa; poly fb; poly la; poly lb;
@@ -1111 +1385 @@
 
 // sort a list of ideals using idbefid
-proc sortlistideals(list L)
+static proc sortlistideals(list L)
 {
   int i; int j; int n;
@@ -1134 +1408 @@
 }
 
-// returns 1 if the two lists of ideals are equal and 0 if not
-proc equallistideals(list L, list M)
-{
-  int t; int i;
-  if (size(L)!=size(M)){return(0);}
-  else
-  {
-    t=1;
-    if (size(L)>0)
-    {
-      i=1;
-      while ((t) and (i<=size(L)))
-      {
-        if (equalideals(L[i],M[i])==0){t=0;}
-        i++;
-      }
-    }
-    return(t);
-  }
+// Crep
+// Computes the C-representation of V(N) \ V(M).
+// input:
+//    ideal N (null ideal) (not necessarily radical nor maximal)
+//    ideal M (hole ideal) (not necessarily containing N)
+// output:
+//    the list (a,b) of the canonical ideals
+//    the Crep of V(N) \ V(M)
+// Assumed to be called in the ring @R
+// Works on the ring @P
+proc Crep(ideal N, ideal M)
+ "USAGE:  Crep(N,M);
+ Input: ideal N (null ideal) (not necessarily radical nor maximal)
+           ideal M (hole ideal) (not necessarily containing N)
+ RETURN:  The canonical C-representation of the locally closed set.
+           [ P,Q ], a pair of radical ideals with P included in Q,
+           representing the set V(P) \ V(Q) = V(N) \ V(M)
+ NOTE:   Operates in a ring R=Q[a] (a=parameters)
+ KEYWORDS: locally closed set, canoncial form
+ EXAMPLE:  Crep; shows an example"
+{
+  list l;
+  ideal Np=std(N);
+  if (equalideals(Np,ideal(1)))
+  {
+    l=ideal(1),ideal(1);
+    return(l);
+  }
+  int i;
+  list L;
+  ideal Q=Np+M;
+  ideal P=ideal(1);
+  L=minGTZ(Np);
+  for(i=1;i<=size(L);i++)
+  {
+    if(idcontains(L[i],Q)==0)
+    {
+      P=intersect(P,L[i]);
+    }
+  }
+  P=std(P);
+  Q=std(radical(Q+P));
+  list T=P,Q;
+  return(T);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  if(defined(Grobcov::@P)){kill Grobcov::@R; kill Grobcov::@P; kill Grobcov::@RP;}
+  ring R=0,(x,y,z),lp;
+  short=0;
+  ideal E=x*(x^2+y^2+z^2-25);
+  ideal N=x*(x-3),y-4;
+  def Cr=Crep(E,N);
+  Cr;
+  def L=Prep(E,N);
+  L;
+  def Cr1=PtoCrep(L);
+  Cr1;
 }
 
@@ -1162 +1476 @@
 // output:
 //    the ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r)));
-//    the Prep of V(N)\V(M)
-// Assumed to work in the ring @P of the parameters
-// Can be used without calling previously setglobalrings().
+//    the Prep of V(N) \ V(M)
+// Assumed to be called in the ring @R
+// Works on the ring @P
 proc Prep(ideal N, ideal M)
+ "USAGE:  Prep(N,M);
+ Input: ideal N (null ideal) (not necessarily radical nor maximal)
+           ideal M (hole ideal) (not necessarily containing N)
+ RETURN: The canonical P-representation of the locally closed set V(N) \ V(M)
+           Output: [ Comp_1, .. , Comp_s ] where
+              Comp_i=[p_i,[p_i1,..,p_is_i]]
+ NOTE:   Operates in a ring R=Q[a]  (a=parameters)
+ KEYWORDS: Locally closed set, Canoncial form
+ EXAMPLE:  Prep; shows an example"
 {
   int te;
@@ -1172 +1495 @@
     return(list(list(ideal(1),list(ideal(1)))));
   }
-  def RR=basering;
-  if(defined(@P)){te=1;}
-  else{te=0; setglobalrings();}
-  setring(@P);
-  ideal Np=imap(RR,N);
-  ideal Mp=imap(RR,M);
   int i; int j; list L0;
-
-  list Ni=minGTZ(Np);
+  list Ni=minGTZ(N);
   list prep;
   for(j=1;j<=size(Ni);j++)
@@ -1190 +1506 @@
   for (i=1;i<=size(Ni);i++)
   {
-    Mij=minGTZ(Ni[i]+Mp);
+    Mij=minGTZ(Ni[i]+M);
     for(j=1;j<=size(Mij);j++)
     {
@@ -1203 +1519 @@
   }
   if (size(prep)==0){prep=list(list(ideal(1),list(ideal(1))));}
-  setring(RR);
-  def pprep=imap(@P,prep);
-  if(te==0){kill @P; kill @R; kill @RP;}
-  return(pprep);
+  def Lout=CompleteA(prep,prep);
+  return(Lout);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  if(defined(Grobcov::@P)){kill Grobcov::@R; kill Grobcov::@P; kill Grobcov::@RP;}
+  short=0;
+  ring R=0,(x,y,z),lp;
+  ideal E=x*(x^2+y^2+z^2-25);
+  ideal N=x*(x-3),y-4;
+  Prep(E,N);
 }
 
@@ -1213 +1537 @@
 // input:
 //    list ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r)));
-//         the P-representation of V(N)\V(M)
+//         the P-representation of V(N) \ V(M)
 // output:
 //    list (ideal ida, ideal idb)
-//    the C-representaion of V(N)\V(M) = V(ida)\V(idb)
-// Assumed to work in the ring @P of the parameters
-// If the routine is to be called from the top, a previous call to
-// setglobalrings() is needed.
+//    the C-representaion of V(N) \ V(M) = V(ida) \ V(idb)
+// Assumed to be called in the ring @R
+// Works on the ring @P
 proc PtoCrep(list L)
+"USAGE: PtoCrep(L)
+ Input L: list  [ Comp_1, .. , Comp_s ] where
+          Comp_i=[p_i,[p_i1,..,p_is_i] ], is
+          the P-representation of a locally closed set V(N) \ V(M)
+ RETURN:  The canonical C-representation of the locally closed set
+          [ P,Q ], a pair of radical ideals with P included in Q,
+          representing the set V(P) \ V(Q)
+ NOTE: Operates in a ring R=Q[a]  (a=parameters)
+ KEYWORDS: locally closed set, canoncial form
+ EXAMPLE:  PtoCrep; shows an example"
+{
+  int te;
+  def RR=basering;
+  if(defined(@P)){te=1; setring(@P); list Lp=imap(RR,L);}
+  else {  te=0; def Lp=L;}
+  def La=PtoCrep0(Lp);
+  if(te==1) {setring(RR); def LL=imap(@P,La);}
+  if(te==0){def LL=La;}
+  return(LL);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  if(defined(Grobcov::@P)){kill Grobcov::@R; kill Grobcov::@P; kill Grobcov::@RP;}
+  short=0;
+  ring R=0,(x,y,z),lp;
+
+  // (P,Q) represents a locally closed set
+  ideal P=x^3+x*y^2+x*z^2-25*x;
+  ideal Q=y-4,x*z,x^2-3*x;
+
+  // Now compute the P-representation=
+  def L=Prep(P,Q);
+  L;
+  // Now compute the C-representation=
+  def J=PtoCrep(L);
+  J;
+  // Now come back recomputing the P-represetation of the C-representation=
+  Prep(J[1],J[2]);
+}
+
+// PtoCrep0
+// Computes the C-representation from the P-representation.
+// input:
+//    list ((p_1,(p_11,p_1k_1)),..,(p_r,(p_r1,p_rk_r)));
+//         the P-representation of V(N) \ V(M)
+// output:
+//    list (ideal ida, ideal idb)
+//    the C-representation of V(N) \ V(M) = V(ida) \ V(idb)
+// Works in a ring Q[u_1,..,u_m] and is called on it.
+static proc PtoCrep0(list L)
 {
   int te=0;
-  def RR=basering;
-  if(defined(@P)){te=1;}
-  else{te=0; setglobalrings();}
-  setring(@P);
-  def Lp=imap(RR,L);
+  def Lp=L;
   int i; int j;
   ideal ida=ideal(1); ideal idb=ideal(1); list Lb; ideal N;
@@ -1241 +1611 @@
     }
   }
+  //idb=radical(idb);
   def La=list(ida,idb);
-  setring(RR);
-  def LL=imap(@P,La);
-  if(te==0){kill @P; kill @R; kill @RP;}
-  return(LL);
+  return(La);
 }
 
@@ -1254 +1622 @@
 //      and dehomogenize the result.
 proc cgsdr(ideal F, list #)
-"USAGE: cgsdr(F); To compute a disjoint, reduced CGS.
+"USAGE: cgsdr(F);
+          F: ideal in Q[a][x] (a=parameters, x=variables) to be discussed.
+          To compute a disjoint, reduced Comprehensive Groebner System (CGS).
           cgsdr is the starting point of the fundamental routine grobcov.
-          Inside grobcov it is used only with options 'can' set to 0,1 and
+          Inside grobcov it is used with options 'can' set to 0,1 and
           not with options ('can',2).
           It is to be used if only a disjoint reduced CGS is required.
-          F: ideal in Q[a][x] (parameters and variables) to be discussed.
 
           Options: To modify the default options, pairs of arguments
@@ -1275 +1644 @@
             \"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).
+                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
@@ -1285 +1654 @@
                 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.
+                is set to (\"out\",0) because it is not compatible.
           One can give none or whatever of these options.
           With the default options (\"can\",2,\"out\",1), only the
-          Kapur-Sun-Wang algorithm is computed. This is very effectif
-          but is only the starting point for the computation.
+          Kapur-Sun-Wang algorithm is computed. This is very efficient
+          but is only the starting point for the computation of grobcov.
           When grobcov is computed, the call to cgsdr inside uses
           specific options that are more expensive ("can",0-1,"out",0).
-RETURN:   Returns a list T describing a reduced and disjoint
+RETURN: Returns a list T describing a reduced and disjoint
           Comprehensive Groebner System (CGS),
           With option (\"out\",0)
-           the segments are grouped by
-           leading power products (lpp) of the reduced Groebner
-           basis and given in P-representation.
-           The returned list is of the form:
+          the segments are grouped by
+          leading power products (lpp) of the reduced Groebner
+          basis and given in P-representation.
+          The returned list is of the form:
            (
              (lpp, (num,basis,segment),...,(num,basis,segment),lpp),
@@ -1304 +1673 @@
              (lpp, (num,basis,segment),...,(num,basis,segment),lpp)
            )
-           The bases are the reduced Groebner bases (after normalization)
-           for each point of the corresponding segment.
-
-           The third element of each lpp segment is the lpp of the
-           used ideal in the CGS as a string:
+          The bases are the reduced Groebner bases (after normalization)
+          for each point of the corresponding segment.
+
+          The third element of each lpp segment is the lpp of the
+          used ideal in the CGS as a string:
             with option (\"can\",0) the homogenized basis is used
             with option (\"can\",1) the homogenized ideal is used
@@ -1314 +1683 @@
 
           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:
+          only KSW is applied and segments are given as
+          difference of varieties and are not grouped
+          The returned list is of the form:
            (
              (E,N,B),..(E,N,B)
            )
-           E is the null variety
-           N is the nonnull variety
-           segment = V(E)\V(N)
-           B is the reduced Groebner basis
-
-NOTE:     The basering R, must be of the form Q[a][x], a=parameters,
-          x=variables, and should be defined previously, and the ideal
+          E is the null variety
+          N is the nonnull variety
+          segment = V(E) \ V(N)
+          B is the reduced Groebner basis
+
+NOTE:  The basering R, must be of the form Q[a][x], (a=parameters,
+          x=variables), and should be defined previously, and the ideal
           defined on R.
 KEYWORDS: CGS, disjoint, reduced, Comprehensive Groebner System
@@ -1337 +1706 @@
   // 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;
-  def E=ideal(0);
-  def N=ideal(1);
-  int comment=0;
   int start=timer;
   list L=#;
@@ -1404 +1773 @@
     {
       // COMPUTING THE HOMOGOENIZED IDEAL
-      if(comment>0){string("Homogenizing the whole ideal: option can=1");}
+      if(comment>=1){string("Homogenizing the whole ideal: option can=1");}
       list RRL=ringlist(RR);
       RRL[3][1][1]="dp";
@@ -1440 +1809 @@
       BH[i]=homog(BH[i],@t);
     }
-    if (comment>=1){string("Homogenized system = "); BH;}
+    if (comment>=2){string("Homogenized system = "); BH;}
     def GSH=KSW(BH,LH);
     setglobalrings();
@@ -1488 +1857 @@
           GS[i][3]=postredgb(mingb(GS[i][3]));
           if (typeof(GS[i][7])=="ideal")
-          {
-            GS[i][7]=postredgb(mingb(GS[i][7]));
-          }
+          { GS[i][7]=postredgb(mingb(GS[i][7]));}
         }
       }
@@ -1499 +1866 @@
 }
 example
-{ "EXAMPLE:"; echo = 2;
-  "Casas conjecture for degree 4";
+{
+  "EXAMPLE:"; echo = 2;
+  // Casas conjecture for degree 4:
   ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp;
+  short=0;
   ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0),
           x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1),
@@ -1514 +1883 @@
 //            lpp segments and improve the output
 // output: grouped segments by lpp obtained in cgsdr
-proc grsegments(list T)
+static proc grsegments(list T)
 {
   int i;
@@ -1545 +1914 @@
 }
 
-// idcontains
-// input: ideal p, ideal q
-// output: 1 if p contains q,  0 otherwise
-// If the routine is to be called from the top, a previous call to
-// setglobalrings() is needed.
-proc idcontains(ideal p, ideal q)
-{
-  int t; int i;
-  t=1; i=1;
-  def RR=basering;
-  setring(@P);
-  def P=imap(RR,p);
-  def Q=imap(RR,q);
-  attrib(P,"isSB",1);
-  poly r;
-  while ((t) and (i<=size(Q)))
-  {
-    r=reduce(Q[i],P);
-    if (r!=0){t=0;}
-    i++;
-  }
-  setring(RR);
-  return(t);
-}
-
-// selectminindeals
-//   given a list of ideals returns the list of integers corresponding
-//   to the minimal ideals in the list
-// input: L (list of ideals)
-// output: the list of integers corresponding to the minimal ideals in L
-// If the routine is to be called from the top, a previous call to
-// setglobalrings() is needed.
-proc selectminideals(list L)
-{
-  if (size(L)==0){return(L)};
-  def RR=basering;
-  setring(@P);
-  def Lp=imap(RR,L);
-  int i; int j; int t; intvec notsel;
-  list P;
-  for (i=1;i<=size(Lp);i++)
-  {
-    if(memberpos(i,notsel)[1])
-    {
-      i++;
-      if(i>size(Lp)){break;}
-    }
-    t=1;
-    j=1;
-    while ((t) and (j<=size(Lp)))
-    {
-      if (i==j){j++;}
-      if ((j<=size(Lp)) and (memberpos(j,notsel)[1]==0))
-      {
-
-        if (idcontains(Lp[i],Lp[j]))
-        {
-          notsel[size(notsel)+1]=i;
-          t=0;
-        }
-      }
-      j++;
-    }
-    if (t){P[size(P)+1]=i;}
-  }
-  setring(RR);
-  return(P);
-}
-
 // LCUnion
-// Given a list of the P-representations of disjoint locally closed segments
+// 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
@@ -1623 +1923 @@
 // output: P-representation of the union
 //       ((P_j,(P_j1,...,P_jk_j | j=1..t)))
-// If the routine is to be called from the top, a previous call to
-// setglobalrings() is needed.
-// Auxiliary routine called by grobcov and addcons
-// A previous call to setglobarings() is needed if it is to be used at the top.
-proc LCUnion(list LL)
+static proc LCUnion(list LL)
 {
   def RR=basering;
   setring(@P);
-  int care;
   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;  intvec Qi;
-  if(not(defined(@Q2))){list @Q2;}
-  else{kill @Q2; list @Q2;}
-  exportto(Top,@Q2);
+  list L0; list P0; list P; list Q0; list Q;
   for (i=1;i<=size(L);i++)
   {
@@ -1647 +1939 @@
   }
   Q0=selectminideals(P0);
-  //"T_3Q0="; Q0;
-  kill P; list P;
   for (i=1;i<=size(Q0);i++)
   {
-    Qi=L0[Q0[i]];
-    Q[size(Q)+1]=Qi;
-      P[size(P)+1]=L[Qi[1]][Qi[2]];
-  }
-  @Q2=Q;
-  if(size(P)==0)
-  {
-    setring(RR);
-    list TT;
-    return(TT);
+    Q[i]=L0[Q0[i]];
+    P[i]=L[Q[i][1]][Q[i][2]];
   }
   // P is the list of the maximal components of the union
@@ -1676 +1958 @@
         {
           C[size(C)+1]=L[i][j];
-          C[size(C)][3]=intvec(i,j);
         }
       }
     }
-    if(size(P[k])>=3){T[size(T)+1]=list(Q[k],P[k][1],addpart(H,C),P[k][3]);}
-    else{T[size(T)+1]=list(Q[k],P[k][1],addpart(H,C));}
-  }
-  @Q2=sortpairs(@Q2);
+    T[size(T)+1]=list(Q[k],P[k][1],addpart(H,C));
+  }
   setring(RR);
   def TT=imap(@P,T);
@@ -1699 +1978 @@
 // posi=(i,j) position of the component
 //        the list of segments to be added to the holes
-proc addpart(list H, list C)
+static proc addpart(list H, list C)
 {
   list Q; int i; int j; int k; int l; int t; int t1;
@@ -1717 +1996 @@
         if (equalideals(q,C[j][1]))
         {
-          @Q2[size(@Q2)+1]=C[j][3];
+          // \\ @Q2[size(@Q2)+1]=C[j][3];
           t=0;
           for (k=1;k<=size(C[j][2]);k++)
           {
             t1=1;
-            //kill addq;
-            //list addq;
             l=1;
             while((t1) and (l<=size(Q)))
@@ -1739 +2016 @@
             {
               addq[size(addq)+1]=C[j][2][k];
-              @Q2[size(@Q2)+1]=C[j][3];
+              // \\ @Q2[size(@Q2)+1]=C[j][3];
             }
           }
@@ -1756 +2033 @@
         list addq;
       }
-      //print("Q="); Q; print("notQ="); notQ;
     }
     i++;
@@ -1775 +2051 @@
 // that part.
 // Works on @P ring.
-proc addpartfine(list H, list C0)
-{
+static proc addpartfine(list H, list C0)
+{
+  //"T_H="; H;
   int i; int j; int k; int te; intvec notQ; int l; list sel; int used;
   intvec jtesC;
   if ((size(H)==1) and (equalideals(H[1],ideal(1)))){return(H);}
   if (size(C0)==0){return(H);}
-  def RR=basering;
-  setring(@P);
   list newQ; list nQ; list Q; list nQ1; list Q0;
-  def Q1=imap(RR,H);
-  //Q1=sortlistideals(Q1);
-  def C=imap(RR,C0);
+  def Q1=H;
+  //Q1=sortlistideals(Q1,idbefid);
+  def C=C0;
   while(equallistideals(Q0,Q1)==0)
   {
@@ -1848 +2123 @@
       }
     }
+    //"T_Q1_0="; Q1;
     sel=selectminideals(Q1);
     kill nQ1; list nQ1;
@@ -1856 +2132 @@
     Q1=nQ1;
   }
-  setring(RR);
+  if(size(Q1)==0){Q1=ideal(1),ideal(1);}
+  //"T_Q1_1="; Q1;
   //if(used>0){string("addpartfine was ", used, " times used");}
-  return(imap(@P,Q1));
-}
+  return(Q1);
+}
+
 
 // Auxiliary routine for grobcov: ideal F is assumed to be homogeneous
@@ -1869 +2147 @@
 //      where a Prep is ( (p1,(p11,..,p1k_1)),..,(pj,(pj1,..,p1k_j)) )
 //            a Crep is ( ida, idb )
-proc gcover(ideal F,list #)
+static proc gcover(ideal F,list #)
 {
   int i; int j; int k; ideal lpp; list GPi2; list pairspP; ideal B; int ti;
@@ -2017 +2295 @@
 // grobcov
 // input:
-//    ideal F: a parametric ideal in Q[a][x], where a are the parameters
-//             and x the variables
+//    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
@@ -2039 +2316 @@
 proc grobcov(ideal F,list #)
 "USAGE:   grobcov(F); This is the fundamental routine of the
-          library. It computes the Groebner cover of a parametric ideal
-          (see (*) Montes A., Wibmer M., Groebner Bases for Polynomial
-          Systems with parameters. JSC 45 (2010) 1391-1425.)
-          The Groebner cover of a parametric ideal consist of a set of
+          library. It computes the Groebner cover of a parametric ideal. See
+                Montes A., Wibmer M.,
+               \"Groebner Bases for Polynomial Systems with parameters\".
+               JSC 45 (2010) 1391-1425.)
+          The Groebner Cover of a parametric ideal consist of a set of
           pairs(S_i,B_i), where the S_i are disjoint locally closed
           segments of the parameter space, and the B_i are the reduced
@@ -2048 +2326 @@
 
           The ideal F must be defined on a parametric ring Q[a][x].
+          (a=parameters, x=variables)
           Options: To modify the default options, pair of arguments
           -option name, value- of valid options must be added to the call.
@@ -2053 +2332 @@
           Options:
             \"null\",ideal E: The default is (\"null\",ideal(0)).
-            \"nonnull\",ideal N: The default (\"nonnull\",ideal(1)).
+            \"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).
+                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
+                Groebner Cover, so that the result is the canonical
+                Groebner Cover. Setting (\"can\",0) only homogenizes the
                 basis so the result is not exactly canonical, but the
                 computation is shorter.
             \"ext\",0-1: The default is (\"ext\",0). With the default
-                (\"ext\",0), only the generic representation is computed
-                (single polynomials, but not specializing to non-zero at
-                each point of the segment. With option (\"ext\",1) the
+                (\"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.
+                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)
@@ -2087 +2367 @@
           of the segment.
           The lpph corresponds to the lpp of the homogenized ideal
-          and is different for each segment. It is given as a string.
-
-          Basis: to each element of lpp corresponds an I-regular function given
-          in full representation (by option (\"ext\",1)) or in
+          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
@@ -2110 +2392 @@
           The P-representation of a segment is of the form
           ((p_1,(p_11,..,p_1k1)),..,(p_r,(p_r1,..,p_rkr))
-          representing the segment U_i (V(p_i) \ U_j (V(p_ij))),
+          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
+          (E,N) representing V(E) \ V(N), and the ideals E and N are
           radical and N contains E.
 
-NOTE: The basering R, must be of the form Q[a][x], a=parameters,
-          x=variables, and should be defined previously. The ideal must
+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
@@ -2244 +2526 @@
 }
 example
-{ "EXAMPLE:"; echo = 2;
-  "Casas conjecture for degree 4";
+{
+  "EXAMPLE:"; echo = 2;
+  // Casas conjecture for degree 4:
   ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp;
+  short=0;
   ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0),
-          x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1),
-          x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0),
-          x2^2+(2*a3)*x2+(a2),
-          x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0),
-          x3+(a3);
+            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);
 }
@@ -2262 +2546 @@
 // Can be called from the top
 proc extend(list GC, list #);
-"USAGE:   extend(GC); When the grobcov of an ideal has been computed
-          with the default option (\"ext\",0) and the explicit option
-          (\"rep\",2) (which is not the default), then one can call
-          extend (GC) (and options) to obtain the full representation
-          of the bases. With the default option (\"ext\",0) only the
-          generic representation of the bases are computed, and one can
-          obtain the full representation using extend.
-            \"rep\",0-1-2: The default is (\"rep\",0) and then the segments
-                are given in canonical P-representation. Option (\"rep\",1)
-                represents the segments in canonical C-representation,
-                and option (\"rep\",2) gives both representations.
-            \"comment\",0-1: The default is (\"comment\",0). Setting
-                \"comment\" higher will provide information about the
-                time used in the computation.
-          One can give none or whatever of these options.
-RETURN:   The list
-          (
-           (lpp_1,basis_1,segment_1,lpph_1),
-           ...
-           (lpp_s,basis_s,segment_s,lpph_s)
-          )
-
-          The lpp are constant over a segment and correspond to the
-          set of lpp of the reduced Groebner basis for each point
-          of the segment.
-          The lpph corresponds to the lpp of the homogenized ideal
-          and is different for each segment. It is given as a string.
-
-          Basis: to each element of lpp corresponds an I-regular function given
-          in full representation. The
-          I-regular function is the corresponding element of the reduced
-          Groebner basis for each point of the segment with the given lpp.
-          For each point in the segment, the polynomial or the set of
-          polynomials representing it, if they do not specialize to 0,
-          then after normalization, specializes to the corresponding
-          element of the reduced Groebner basis. In the full representation
-          at least one of the polynomials representing the I-regular
-          function specializes to non-zero.
-
-          With the default option (\"rep\",0) the segments are given
-          in P-representation.
-          With option (\"rep\",1) the segments are given
-          in C-representation.
-          With option (\"rep\",2) both representations of the segments are
-          given.
-
-          The P-representation of a segment is of the form
-          ((p_1,(p_11,..,p_1k1)),..,(p_r,(p_r1,..,p_rkr))
-          representing the segment U_i (V(p_i) \ U_j (V(p_ij))),
-          where the p's are prime ideals.
-
-          The C-representation of a segment is of the form
-          (E,N) representing V(E)\V(N), and the ideals E and N are
-          radical and N contains E.
-
-NOTE: The basering R, must be of the form Q[a][x], a=parameters,
-          x=variables, and should be defined previously. The ideal must
+"USAGE: extend(GC); The default option of grobcov provides
+          the bases in generic representation (the I-regular functions
+          of the bases ara given by a single polynomial. It can specialize
+          to zero for some points of the segments, but in general, it
+          is sufficient for many pouposes. Nevertheless the I-regular
+          functions allow a full representation given bey a set of
+          polynomials specializing to the value of the function (after normalization)
+          or to zero, but at least one of the polynomials specializes to non-zero.
+          The full representation can be obtained by computing the
+          grobcov with option \"ext\",1. The default option is \"ext\",0.
+          With option \"ext\",1 the computation can be much more
+          time consuming, even if the result can be simpler.
+          Alternatively, one can compute the full representation of the
+          bases after computing grobcov with the defaoult option \"ext\",0
+          and the option \"rep\",2, that outputs both the Prep and the Crep
+          of the segments and then call \"extend\" to the output.
+
+RETURN: When calling extend(grobcov(S,\"rep\",2)) the result is of the form
+              (
+                 (lpp_1,basis_1,segment_1,lpph_1),
+                  ...
+                 (lpp_s,basis_s,segment_s,lpph_s)
+              )
+          where each function of the basis can be given by an ideal
+          of representants.
+
+NOTE: The basering R, must be of the form Q[a][x], (a=parameters,
+          x=variables), and should be defined previously. The ideal must
           be defined on R.
 KEYWORDS: Groebner cover, parametric ideal, canonical, discussion of
@@ -2336 +2591 @@
   }
   if(m==1){"Warning! grobcov has already extended bases"; return(S);}
-  if(size(GC[1])!=5){"Warning! extend make sense only when grobcov has been called with options 'rep',2,'ext',0"; " "; return();}
+  if(size(GC[1])!=5){"Warning! extend make sense only when grobcov has been called with options  'rep',2,'ext',0"; " "; return();}
   int repop=0;
   int start3=timer;
@@ -2434 +2689 @@
 example
 {
-  ring R=(0,a0,b0,c0,a1,b1,c1,a2,b2,c2),(x), dp;
+  "EXAMPLE:"; echo = 2;
+  ring R=(0,a0,b0,c0,a1,b1,c1),(x), dp;
   short=0;
   ideal S=a0*x^2+b0*x+c0,
-          a1*x^2+b1*x+c1,
-          a2*x^2+b2*x+c2;
-  "System S="; S;
-
-  def GCS=grobcov(S,"rep",2,"comment",1);
-  "grobcov(S,'rep',2,'comment',1)="; GCS;
-  def FGC=extend(GCS,"rep",0,"comment",1);
-  "Full representation="; FGC;
+            a1*x^2+b1*x+c1;
+  def GCS=grobcov(S,"rep",2);
+  GCS;
+  def FGC=extend(GCS,"rep",0);
+  // Full representation=
+  FGC;
 }
 
@@ -2450 +2704 @@
 // nonzerodivisor
 // input:
-//    poly g in K[a],
+//    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
-proc nonzerodivisor(poly gr, list Pr)
+static proc nonzerodivisor(poly gr, list Pr)
 {
   def RR=basering;
@@ -2492 +2746 @@
 // input:
 //   int i:
-//   list LPr: (p1,..,pr) of prime components of an ideal in K[a]
+//   list LPr: (p1,..,pr) of prime components of an ideal in Q[a]
 // output:
 //   list (fr,fnr) of two polynomials that are equal on V(pi)
 //       and fr=0 on V(P) \ V(pi), and fnr is nonzero on V(pj) for all j.
-proc deltai(int i, list LPr)
+static proc deltai(int i, list LPr)
 {
   def RR=basering;
@@ -2502 +2756 @@
   def LP=imap(RR,LPr);
   int j; poly p;
-  ideal F=ideal(1);
+  def F=ideal(1);
   poly f;
   poly fn;
@@ -2539 +2793 @@
 // output:
 //    poly P on an open and dense set of V(p_1 int ... p_r)
-proc combine(list L, ideal F)
+static proc combine(list L, ideal F)
 {
   // ATTENTION REVISE AND USE Pci and F
@@ -2553 +2807 @@
 }
 
-// Auxiliary routine
-// elimconstfac: eliminate the factors in the polynom f that are in K[a]
-// input:
-//   poly f:
-//   list L: of components of the segment
-// output:
-//   poly f2  where the factors of f in K[a] that are non-null on any component
-//   have been dropped from f
-proc elimconstfac(poly f)
-{
-  int cond; int i; int j; int t;
-  if (f==0){return(f);}
-  def RR=basering;
-  setring(@R);
-  def ff=imap(RR,f);
-  list l=factorize(ff,0);
-  poly f1=1;
-  for(i=2;i<=size(l[1]);i++)
-  {
-      f1=f1*(l[1][i])^(l[2][i]);
-  }
-  setring(RR);
-  def f2=imap(@R,f1);
-  return(f2);
-}
 
 //Auxiliary routine
@@ -2587 +2816 @@
 // output:
 //   t:  with value 1 if f reduces modulo P, 0 if not.
-proc nullin(poly f,ideal P)
+static proc nullin(poly f,ideal P)
 {
   int t;
@@ -2605 +2834 @@
 // Input: A polynomial f
 // Output: The list of leading terms
-proc monoms(poly f)
+static proc monoms(poly f)
 {
   list L;
@@ -2628 +2857 @@
 //   poly f: a generic polynomial in the basis
 //   ideal idp: such that ideal(S)=idp
-//   ideal idq: such that S=V(idp)\V(idq)
+//   ideal idq: such that S=V(idp) \ V(idq)
 ////   NW the list of ((N1,W1),..,(Ns,Ws)) of red-rep of the grouped
 ////      segments in the lpp-segment  NO MORE USED
 // output:
-proc extend0(poly f, ideal idp, ideal idq)
+static proc extend0(poly f, ideal idp, ideal idq)
 {
   matrix CC; poly Q; list NewMonoms;
@@ -2692 +2921 @@
 //        each intvec v=(i_1,..,is) corresponds to a polynomial in the sheaf
 //        that will be built from it in extend procedure.
-proc findindexpolys(list coefs)
+static proc findindexpolys(list coefs)
 {
   int i; int j; intvec numdens;
@@ -2733 +2962 @@
 
 // Auxiliary routine
-// extendcoef: given Q,P in K[a] where P/Q specializes on an open and dense subset
+// 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,
-proc extendcoef(poly P, poly Q, ideal idp, ideal idq)
+static proc extendcoef(poly P, poly Q, ideal idp, ideal idq)
 {
   def RR=basering;
@@ -2765 +2994 @@
 // 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)
+//      (we assume that one of them is non-null on V(N) \ V(M))
+//   ideal N, ideal M: ideals representing the locally closed set V(N) \ V(M)
 // assume to work in @P
-proc selectregularfun(matrix CC, ideal NN, ideal MM)
+static proc selectregularfun(matrix CC, ideal NN, ideal MM)
 {
   int numcombused;
@@ -2827 +3056 @@
 // output:
 //   object T with index c
-proc searchinlist(intvec c,list L)
+static proc searchinlist(intvec c,list L)
 {
   int i; list T;
@@ -2841 +3070 @@
 }
 
-// Auxiliary routine
-// comb: the list of combinations of elements (1,..n) of order p
-proc comb(int n, int p)
-{
-  list L; list L0;
-  intvec c; intvec d;
-  int i; int j; int last;
-  if ((n<0) or (n<p))
-  {
-    return(L);
-  }
-  if (p==1)
-  {
-    for (i=1;i<=n;i++)
-    {
-      c=i;
-      L[size(L)+1]=c;
-    }
-    return(L);
-  }
-  else
-  {
-    L0=comb(n,p-1);
-    for (i=1;i<=size(L0);i++)
-    {
-      c=L0[i]; d=c;
-      last=c[size(c)];
-      for (j=last+1;j<=n;j++)
-      {
-        d[size(c)+1]=j;
-        L[size(L)+1]=d;
-      }
-    }
-    return(L);
-  }
-}
 
 // Auxiliary routine
@@ -2895 +3088 @@
 // The selection is done to obtian the minimal number of elements
 //    of the sheaf that specializes to non-null everywhere.
-proc selectminsheaves(list L)
+static proc selectminsheaves(list L)
 {
   list C=allsheaves(L);
@@ -2909 +3102 @@
 //   list LL of the subsets of C that cover all the subsegments
 //   (the union of the corresponding L(C) has all 1).
-proc smsheaves(list C, list L)
+static proc smsheaves(list C, list L)
 {
   int i; int i0; intvec W;
@@ -2950 +3143 @@
 //    LL is the list of all combrep
 //    LLS is the list of intvec of the corresponding elements of LL
-proc allsheaves(list L)
+static proc allsheaves(list L)
 {
   intvec V; list LL; intvec W; int r; intvec U;
@@ -2986 +3179 @@
 // Output:
 //   int nor: the nuber of 1 of v in the positions given by pos.
-proc numones(intvec v, intvec pos)
+static proc numones(intvec v, intvec pos)
 {
   int i; int n;
@@ -2994 +3187 @@
   }
   return(n);
-}
-
-// Auxiliary routine
-// pos
-// Input:  intvec p of zeros and ones
-// Output: intvec W of the positions where p has ones.
-proc pos(intvec p)
-{
-  int i;
-  intvec W; int j=1;
-  for (i=1; i<=size(p); i++)
-  {
-    if (p[i]==1){W[j]=i; j++;}
-  }
-  return(W);
 }
 
@@ -3019 +3197 @@
 //   intvec pp: of zeroes and ones, where a 0 stays in pp[i] if either
 //   already p[i]==0 or c[i]==1.
-proc actualize(intvec p, intvec c)
+static proc actualize(intvec p, intvec c)
 {
   int i; intvec pp=p;
@@ -3029 +3207 @@
 }
 
+
 // 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
-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);
-}
-
-// Auxiliary routine
-proc reducemodN(poly f,ideal E)
+static proc reducemodN(poly f,ideal E)
 {
   def RR=basering;
@@ -3081 +3226 @@
 // Auxiliary routine
 // intersp: computes the intersection of the ideals in S in @P
-proc intersp(list S)
+static proc intersp(list S)
 {
   def RR=basering;
@@ -3094 +3239 @@
 // Auxiliary routine
 // radicalmember
-proc radicalmember(poly f,ideal ida)
+static proc radicalmember(poly f,ideal ida)
 {
   int te;
@@ -3115 +3260 @@
 }
 
-// Auxiliary routine
-// NonNull: returns 1 if the poly f is nonnull on V(E)\V(N), 0 otherwise.
-proc NonNull(poly f, ideal E, ideal N)
-{
-  int te=1; int i;
-  def RR=basering;
-  setring(@P);
-  def fp=imap(RR,f);
-  def Ep=imap(RR,E);
-  def Np=imap(RR,N);
-  ideal H;
-  ideal Ef=Ep+fp;
-  for (i=1;i<=size(Np);i++)
-  {
-    te=radicalmember(Np[i],Ef);
-    if (te==0){break;}
-  }
-  setring(RR);
-  return(te);
-}
+// // Auxiliary routine
+// // NonNull: returns 1 if the poly f is nonnull on V(E) \ V(N), 0 otherwise.
+// // Input:
+// //   f: polynomial
+// //   E: null ideal
+// //   N: nonnull ideal
+// // Output:
+// //   1 if f is nonnul in V(P) \ V(Q),
+// //   0 if it has zeroes inside.
+// //  Works in @P
+// proc NonNull(poly f, ideal E, ideal N)
+// {
+//   int te=1; int i;
+//   def RR=basering;
+//   setring(@P);
+//   def fp=imap(RR,f);
+//   def Ep=imap(RR,E);
+//   def Np=imap(RR,N);
+//   ideal H;
+//   ideal Ef=Ep+fp;
+//   for (i=1;i<=size(Np);i++)
+//   {
+//     te=radicalmember(Np[i],Ef);
+//     if (te==0){break;}
+//   }
+//   setring(RR);
+//   return(te);
+// }
 
 // Auxiliary routine
@@ -3149 +3302 @@
 //            points where the q's are null on S.
 //            The elements of caout are of the form (p,q,prep);
-proc selectextendcoef(matrix CC, ideal ida, ideal idb)
+static proc selectextendcoef(matrix CC, ideal ida, ideal idb)
 {
   def RR=basering;
@@ -3196 +3349 @@
 
 // Auxiliary routine
-// input:
+// 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
-proc plusP(ideal E1,ideal E2)
+// Output:
+//   ideal Ep=E1+E2; computed in @P
+static proc plusP(ideal E1,ideal E2)
 {
   def RR=basering;
@@ -3221 +3375 @@
 //      All the vi without zeroes are in outcomb, and those with zeroes are
 //         combined to form new intvec with the rest
-proc reform(list combpolys, intvec numdens)
+static proc reform(list combpolys, intvec numdens)
 {
   list combp0; list combp1; int i; int j; int k; int l; list rest; intvec notfree;
@@ -3296 +3450 @@
 }
 
-// Auxiliary routine
-// nonnullCrep
-proc nonnullCrep(poly f0,ideal ida0,ideal idb0)
-{
-  int i;
-  def RR=basering;
-  setring(@P);
-  def f=imap(RR,f0);
-  def ida=imap(RR,ida0);
-  def idb=imap(RR,idb0);
-  def idaf=ida+f;
-  int te=1;
-  for(i=1;i<=size(idb);i++)
-  {
-    if(radicalmember(idb[i],idaf)==0)
-    {
-      te=0; break;
-    }
-  }
-  setring(RR);
-  return(te);
-}
 
 // Auxiliary routine
@@ -3326 +3458 @@
 //             L=(p1,..,ps);  F0=(f1,..,fs);
 //             F0[i] \in intersect_{j#i} p_i
-proc precombint(list L)
+static proc precombint(list L)
 {
   int i; int j; int tes;
@@ -3369 +3501 @@
 // Auxiliary routine
 // minAssGTZ eliminating denominators
-proc minGTZ(ideal N);
+static proc minGTZ(ideal N);
 {
   int i; int j;
@@ -3389 +3521 @@
 // Input:
 //   ideal E: of null conditions
-//   ideal N: of non-null conditions representing V(E)\V(N)
+//   ideal N: of non-null conditions representing V(E) \ V(N)
 // Output:
-//   1 if V(E) \V(N) = empty
+//   1 if V(E) \ V(N) = empty
 //   0 if not
-proc inconsistent(ideal E, ideal N)
+static proc inconsistent(ideal E, ideal N)
 {
   int j;
@@ -3423 +3555 @@
 // Auxiliary routine
 // MDBasis: Minimal Dickson Basis
-proc MDBasis(ideal G)
+static proc MDBasis(ideal G)
 {
   int i; int j; int te=1;
@@ -3449 +3581 @@
 // Auxiliary routine
 // primepartZ
-proc primepartZ(poly f);
-{
-  def R=basering;
+static proc primepartZ(poly f);
+{
   def cp=content(f);
   def fp=f/cp;
@@ -3458 +3589 @@
 
 // LCMLC
-proc LCMLC(ideal H)
+static proc LCMLC(ideal H)
 {
   int i;
@@ -3498 +3629 @@
 //                    but without canonical description of the segments.
 //     ((B,E,N),..,(B,E,N))
-proc KSW(ideal F, list #)
+static proc KSW(ideal F, list #)
 {
   setglobalrings();
@@ -3570 +3701 @@
 // Auxiliary routine
 // sqf
-proc sqf(poly f)
+static proc sqf(poly f)
 {
   def RR=basering;
@@ -3589 +3720 @@
 //   The ideal must be defined on C[parameters][variables]
 // Output:
-proc KSW0(ideal F, ideal E, ideal N, int comment)
+static proc KSW0(ideal F, ideal E, ideal N, int comment)
 {
   def R=basering;
@@ -3612 +3743 @@
   E0=imap(@P,E1);
   N0=imap(@P,N1);
-//   E0=elimrepeated(E0);
-//   N0=elimrepeated(N0);
   if (inconsistent(E0,N0)==1)
   {
@@ -3711 +3840 @@
 // Auxiliary routine
 // KSWtocgsdr
-proc KSWtocgsdr(list L)
+static proc KSWtocgsdr(list L)
 {
   int i; list CG; ideal B; ideal lpp; int j; list NKrep;
@@ -3737 +3866 @@
 //    the ((p_1,(p_11,..,p_1k_1)),..,(p_r,(p_r1,..,p_rk_r)));
 //    the Prep of V(N) \ V(W)
-proc KtoPrep(ideal N, ideal W)
+static proc KtoPrep(ideal N, ideal W)
 {
   int i; int j;
@@ -3781 +3910 @@
 // input:  the list of vertices of KSW
 // output: the same terminal vertices grouped by lpp
-proc groupKSWsegments(list T)
-{
-  //"T_T="; T;
+static proc groupKSWsegments(list T)
+{
   int i; int j;
   list L;
@@ -3814 +3942 @@
 
 // indepparameters
-// Auxiliary routine to detect "Special" components of the locus
+// Auxiliary routine to detect 'Special' components of the locus
 // Input: ideal B
 // Output:
 //   1 if the solutions of the ideal do not depend on the parameters
 //   0 if they depend
-proc indepparameters(ideal B)
+static proc indepparameters(ideal B)
 {
   def R=basering;
@@ -3837 +3965 @@
 // if the dimension in @P of an ideal in the parameters has dimension 0 then it returns 0
 // else it retuns 1
-proc dimP0(ideal N)
+static proc dimP0(ideal N)
 {
   def R=basering;
@@ -3850 +3978 @@
 }
 
-// Auxiliary routine.
-// input:    ideals E and F (assumed in ring @P
-// returns: 1 if ideal E is contained in ideal F (assumed F is std basis)
-//              0 if not
-proc containedideal(ideal E, ideal F)
-{
-  int i; int t; poly f;
-  if(equalideals(F,ideal(0)))
-  {
-    if(equalideals(E,ideal(0))==0){return(0);}
-    else(return(1));
-  }
-  t=1; i=1;
-  while((t==1) and (i<=size(E)))
-  {
-    attrib(F,"isSB",1);
-    f=reduce(E[i],F);
-    if(f!=0){t=0;}
-    i++;
-  }
-  return(t);
-}
-
-// addcons: given a set of locally closed sets given in P-representation,
-//       (particularly a subset of components of a selection of segments
-//       of the Grobner cover), it builds the canonical P-representation
-//       of the corresponding constructible set, of its union, including levels it they are.
-// input: a list L of disjoint segments (for exmple a selection of segments of the Groebner cover)
-//       of the form
-//       ( ( (p1_1,(p1_11,..,p1_1k_1).. (p1_s,(p1_s1,..,p1_sk_s)),..,
-//         ( (pn_1,(pn_11,..,pn_1j_1).. (pn_s,(pn_s1,..,pn_sj_s))   )
-// output: The canonical P-representation of the  constructible set resulting by
-//       adding the given components.
-//       The first element of each component can be ignored. It is only for internal use.
-//       The fifth element of each component is the level of the constructible set
-//       of the component.
-//       It is no need a previous call to setglobalrings()
-proc addcons(list L)
-"USAGE:   addcons(  ( ( (p1_1,(p1_11,..,p1_1k_1).. (p1_s,(p1_s1,..,p1_sk_s)),..,
-  ( (pn_1,(pn_11,..,pn_1j_1).. (pn_s,(pn_s1,..,pn_sj_s))   )   )
-  a list L of locally closed sets in P-representation
-RETURN: the canonical P-representation of the constructible set of the union.
-NOTE:   It is called internally by the routines locus, locusdg,
-
-KEYWORDS: P-representation, constructible set
-EXAMPLE:  adccons; shows an example"
-{
-  int te;
-  if(defined(@P)){te=1;}
-  else{setglobalrings();}
-  list notadded; list P0;
-  int i; int j; int k; int l;
-  intvec v; intvec v1; intvec v0;
-  ideal J; list K; list L1;
-  for(i=1;i<=size(L);i++)
-  {
-    if(equalideals(L[i][1][1],ideal(1))==0)
-    {L1[size(L1)+1]=L[i];}
-  }
-  L=L1;
-  for(i=1;i<=size(L);i++)
-  {
-    for(j=1;j<=size(L[i]);j++)
-    {
-        v=i,j;
-        notadded[size(notadded)+1]=v;
-    }
-  }
-  //"T_notadded="; notadded;
-  int level=1;
-  list P=L;
-  list LL;
-  list Ln;
-  //int cont=0;
-  while(size(notadded)>0)
-  {
-    kill notadded; list notadded;
-    for(i=1;i<=size(P);i++)
-    {
-      for(j=1;j<=size(P[i]);j++)
-      {
-          v=i,j;
-          notadded[size(notadded)+1]=v;
-      }
-    }
-    //"T_notadded="; notadded;
-     //cont++;
-     P0=P;
-     Ln=LCUnion(P);
-     //"T_Ln="; Ln;
-     notadded=minuselements(notadded,@Q2);
-     //"T_@Q2="; @Q2;
-     //"T_notadded="; notadded;
-     for(i=1;i<=size(Ln);i++)
-     {
-       //Ln[i][4]=  ; JA HAURIA DE VENIR DE LCUnion
-       Ln[i][5]=level;
-       LL[size(LL)+1]=Ln[i];
-     }
-     i=0;j=0;
-     v=0,0;
-     v0=0,0;
-     kill P; list P;
-     for(l=1;l<=size(notadded);l++)
-     {
-       v1=notadded[l];
-       if(v1[1]<>v0[1]){i++;j=1; P[i]=K;} // list P[i];
-       else
-       {
-         if(v1[2]<>v0[2])
-         {
-           j=j+1;
-         }
-       }
-       v=i,j;
-       //"T_v1="; v1;
-       P[i][j]=K; P[i][j][1]=J; P[i][j][2]=K;
-       P[i][j][1]=P0[v1[1]][v1[2]][1];
-       //"T_P0[v1[1]][v1[2]][2]="; P0[v1[1]][v1[2]][2];
-       //"T_size(P0[v1[1]][v1[2]][2])="; size(P0[v1[1]][v1[2]][2]);
-       for(k=1;k<=size(P0[v1[1]][v1[2]][2]);k++)
-       {
-         P[i][j][2][k]=P0[v1[1]][v1[2]][2][k];
-         //string("T_P[",i,"][",j,"][2][",k,"]="); P[i][j][2][k];
-       }
-       v0=v1;
-       //"T_P_1="; P;
-     }
-     //"T_P="; P;
-     level++;
-     //if(cont>3){break;}
-     //kill notadded; list notadded;
-  }
-  if(defined(@Q2)){kill @Q2;}
-  if(te==0){kill @P; kill @R; kill @RP;}
-  //"T_LL_sortida addcons="; LL; "Fi sortida";
-  return(LL);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-   ring R=(0,a,b),(x1,y1,x2,y2,p,r),dp;
-   ideal S93=(a+1)^2+b^2-(p+1)^2,
-                    (a-1)^2+b^2-(1-r)^2,
-                    a*y1-b*x1-y1+b,
-                    a*y2-b*x2+y2-b,
-                    -2*y1+b*x1-(a+p)*y1+b,
-                    2*y2+b*x2-(a+r)*y2-b,
-                    (x1+1)^2+y1^2-(x2-1)^2-y2^2;
-  short=0;
-  "System S93="; S93; " ";
-  def GC93=grobcov(S93);
-  "grobcov(S93)="; GC93; " ";
-  int i;
-  list H;
-  for(i=1;i<=size(GC93);i++){H[i]=GC93[i][3];}
-  "H="; H;
-  "addcons(H)="; addcons(H);
-}
-
 // Takes a list of intvec and sorts it and eliminates repeated elements.
-// Auxiliary routine used in addcons
-proc sortpairs(L)
+// Auxiliary routine
+static proc sortpairs(L)
 {
   def L1=sort(L);
@@ -4020 +3988 @@
 
 // Eliminates the pairs of L1 that are also in L2.
-// Auxiliary routine used in addcons
-proc minuselements(list L1,list L2)
+// Auxiliary routine
+static proc minuselements(list L1,list L2)
 {
   int i;
@@ -4032 +4000 @@
 }
 
+// NorSing
+// Input:
+//   ideal B: the basis of a component of the grobcov
+//   ideal E: the top of the component (assumed to be of dimension > 0 (single equation)
+//   ideal N: the holes of the component
+// Output:
+//   int d: the dimension of B on the variables reduced by the holes.
+//     if d>0    then the component is 'Normal'
+//     if d==0 then the component is 'Singular'
+static proc NorSing(ideal B, ideal E, ideal N, list #)
+{
+  int i; int j; int Fenv=0; int env; int dd;
+  list DD=#;
+  def RR=basering;
+  int moverdim=2;
+  int version=0;
+  int nv=nvars(RR);
+  if(nv<4){version=1;}
+  int d;
+  poly F;
+  for(i=1;i<=(size(DD) div 2);i++)
+  {
+    if(DD[2*i-1]=="movdim"){moverdim=DD[2*i];}
+    if(DD[2*i-1]=="version"){version=DD[2*i];}
+    if(DD[2*i-1]=="family"){F=DD[2*i];}
+  }
+  if(F!=0){Fenv=1;}
+  //"T_B="; B; "E="; E; "N="; N;
+  list L0;
+  if(dimP0(E)==0){L0=2,"Normal";} // 2 es fals pero ha de ser >0 encara que sigui 0
+  else
+  {
+    if(version==0)
+    {
+      //"T_B="; B;  // Computing std(B+E,plex(x,y,x1,..xn)) one can detect if there is a first part
+      // independent of parameters giving the variables with dimension 0
+     dd=indepparameters(B);
+      if (dd==1){d=0; L0=d,string(B),determineF(B,F,E);}
+      else{d=1; L0=2,"Normal";}
+    }
+    else
+    {
+      def RH=ringlist(RR);
+      //"T_RH="; RH;
+      def H=RH;
+      H[1]=0;
+      H[2]=RH[1][2]+RH[2];
+      int n=size(H[2]);
+      intvec ll;
+      for(i=1;i<=n;i++)
+      {
+        ll[i]=1;
+      }
+      H[3][1][1]="lp";
+      H[3][1][2]=ll;
+      def RRH=ring(H);
+      setring(RRH);
+      ideal BH=imap(RR,B);
+      ideal EH=imap(RR,E);
+      ideal NH=imap(RR,N);
+      if(Fenv==1){poly FH=imap(RR,F);}
+      for(i=1;i<=size(EH);i++){BH[size(BH)+1]=EH[i];}
+      BH=std(BH);  // MOLT COSTOS!!!
+      ideal G;
+      ideal r; poly q;
+      for(i=1;i<=size(BH);i++)
+      {
+        r=factorize(BH[i],1);
+        q=1;
+        for(j=1;j<=size(r);j++)
+        {
+          if((pdivi(r[j],NH)[1] != 0) or (equalideals(ideal(NH),ideal(1))))
+          {
+            q=q*r[j];
+          }
+        }
+        if(q!=1){G[size(G)+1]=q;}
+      }
+      setring RR;
+      def GG=imap(RRH,G);
+      ideal GGG;
+      if(defined(L0)){kill L0; list L0;}
+      for(i=1;i<=size(GG);i++)
+      {
+        if(indepparameters(GG[i])){GGG[size(GGG)+1]=GG[i];}
+      }
+      GGG=std(GGG);
+      ideal GLM;
+      for(i=1;i<=size(GGG);i++)
+      {
+        GLM[i]=leadmonom(GGG[i]);
+      }
+      attrib(GLM,"IsSB",1);
+      d=dim(std(GLM));
+      string antiim=string(GGG);
+      L0=d,antiim;
+      if(d==0)
+      {
+        //" ";string("Antiimage of Special component = ", GGG);
+         if(Fenv==1)
+        {
+          L0[3]=determineF(GGG,F,E);
+        }
+      }
+      else
+      {
+        L0[2]="Normal";
+      }
+    }
+  }
+  return(L0);
+}
+
+static proc determineF(ideal A,poly F,ideal E)
+{
+  int env; int i;
+  def RR=basering;
+  def RH=ringlist(RR);
+  def H=RH;
+  H[1]=0;
+  H[2]=RH[1][2]+RH[2];
+  int n=size(H[2]);
+  intvec ll;
+  for(i=1;i<=n;i++)
+  {
+    ll[i]=1;
+  }
+  H[3][1][1]="lp";
+  H[3][1][2]=ll;
+  def RRH=ring(H);
+
+        //" ";string("Antiimage of Special component = ", GGG);
+
+   setring(RRH);
+   list LL;
+   def AA=imap(RR,A);
+   def FH=imap(RR,F);
+   def EH=imap(RR,E);
+   ideal M=std(AA+FH);
+   def rh=reduce(EH,M);
+   if(rh==0){env=1;} else{env=0;}
+   setring RR;
+          //L0[3]=env;
+    return(env);
+}
+
+// DimPar
+// Auxilliary routine to NorSing determining the dimension of a parametric ideal
+// Does not use @P and define it directly because is assumes that
+//                             variables and parameters have been inverted
+static proc DimPar(ideal E)
+{
+  def RRH=basering;
+  def RHx=ringlist(RRH);
+  def Prin=ring(RHx[1]);
+  setring(Prin);
+  def E2=std(imap(RRH,E));
+  def d=dim(E2);
+  setring RRH;
+  return(d);
+}
+
 // locus0(G): Private routine used by locus (the public routine), that
 //                builds the diferent components.
 // input:      The output G of the grobcov (in generic representation, which is the default option)
+//       Options: The algorithm allows the following options ar pair of arguments:
+//                "movdim", d  : by default movdim is 2 but it can be set to other values
+//                "version", v   :  There are two versions of the algorithm. ('version',1) is
+//                 a full algorithm that always distinguishes correctly between 'Normal'
+//                 and 'Special' components, whereas ('version',0) can decalre a component
+//                 as 'Normal' being really 'Special', but is more effective. By default ('version',1)
+//                 is used when the number of variables is less than 4 and 0 if not.
+//                 The user can force to use one or other version, but it is not recommended.
+//                 "system", ideal F: if the initial systrem is passed as an argument. This is actually not used.
+//                 "comments", c: by default it is 0, but it can be set to 1.
+//                 Usually locus problems have mover coordinates, variables and tracer coordinates.
+//                 The mover coordinates are to be placed as the last variables, and by default,
+//                 its number is 2. If one consider locus problems in higer dimensions, the number of
+//                 mover coordinates (placed as the last variables) is to be given as an option.
 // output:
 //         list, the canonical P-representation of the Normal and Non-Normal locus:
@@ -4046 +4190 @@
 //              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.
-proc locus0(list GG, list #)
-{
-  int t1=1; int t2=1;
+static proc locus0(list GG, list #)
+{
+  int te=0;
+  int t1=1; int t2=1; int i;
   def R=basering;
+  //if(defined(@P)==1){te=1; kill @P; kill @R; kill @RP; }
+  //setglobalrings();
+  // Options
+  list DD=#;
   int moverdim=nvars(R);
+  int version=0;
+  int nv=nvars(R);
+  if(nv<4){version=1;}
+  int comment=0;
+  ideal Fm;
+  poly F;
+  for(i=1;i<=(size(DD) div 2);i++)
+  {
+    if(DD[2*i-1]=="movdim"){moverdim=DD[2*i];}
+    if(DD[2*i-1]=="version"){version=DD[2*i];}
+    if(DD[2*i-1]=="system"){Fm=DD[2*i];}
+    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
+    if(DD[2*i-1]=="family"){F=DD[2*i];}
+  }
   list HHH;
   if (GG[1][1][1]==1 and GG[1][2][1]==1 and GG[1][3][1][1][1]==0 and GG[1][3][1][2][1]==1){return(HHH);}
-  list Lop=#;
-  if(size(Lop)>0){moverdim=Lop[1];}
-  setglobalrings();
-  list G1; list G2;
+   list G1; list G2;
   def G=GG;
   list Q1; list Q2;
-  int i; int d; int j; int k;
+  int d; int j; int k;
   t1=1;
   for(i=1;i<=size(G);i++)
@@ -4115 +4275 @@
   }
   else{list P1;}
-  ideal BB;
+  ideal BB; int dd; list NS;
   for(i=1;i<=size(P1);i++)
   {
-    if (indepparameters(P1[i][3])==1){P1[i][3]="Special";}
-    else{P1[i][3]="Normal";}
+       NS=NorSing(P1[i][3],P1[i][1],P1[i][2][1],DD);
+       dd=NS[1];
+       if(dd==0){P1[i][3]=NS;} //"Special";
+       else{P1[i][3]="Normal";}
   }
   list P2;
@@ -4133 +4295 @@
   for(i=1;i<=size(P2);i++){Q2[i]=l; Q2[i][1]=P2[i];} P2=Q2;
 
-  setring @P;
+  setring(@P);
   ideal J;
   if(t1==1)
   {
     def C1=imap(R,P1);
-    def L1=addcons(C1);
+    def L1=AddLocus(C1);
    }
   else{list C1; list L1; kill P1; list P1;}
@@ -4144 +4306 @@
   {
     def C2=imap(R,P2);
-    def L2=addcons(C2);
-  }
+    def L2=AddLocus(C2);
+ }
   else{list L2; list C2; kill P2; list P2;}
     for(i=1;i<=size(L2);i++)
@@ -4169 +4331 @@
   def L=imap(@P,LN);
   for(i=1;i<=size(L);i++){if(size(L[i][2])==0){L[i][2]=ideal(1);}}
-  kill @R; kill @RP; kill @P;
+  //if(te==0){kill @R; kill @RP; kill @P;}
   list LL;
   for(i=1;i<=size(L);i++)
   {
-    LL[i]=l;
-    LL[i]=elimfromlist(L[i],1);
-  }
-  return(LL);
-}
-
-// locus0dg(G): Private routine used by locusdg (the public routine), that
-//                builds the diferent components used in Dynamical Geometry
-// input:      The output G of the grobcov (in generic representation, which is the default option)
-// output:
-//         list, the canonical P-representation of the Relevant components of the locus in Dynamical Geometry,
-//              i.e. the Normal and Accumulation components.
-//              This routine is compemented by locusdg 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)
-//              whwre type_i is always "Relevant".
-//         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.
-proc locus0dg(list GG, list #)
-{
-  int t1=1; int t2=1; int ojo;
-  def R=basering;
-  int moverdim=nvars(R);
-  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 Lop=#;
-  if(size(Lop)>0){moverdim=Lop[1];}
-  setglobalrings();
-  list G1; list G2;
-  def G=GG;
-  list Q1; list Q2;
-  int i; int d; int j; int k;
-  t1=1;
-  for(i=1;i<=size(G);i++)
-  {
-    attrib(G[i][1],"IsSB",1);
-    d=locusdim(G[i][2],moverdim);
-    if(d==0){G1[size(G1)+1]=G[i];}
-    else
-    {
-      if(d>0){G2[size(G2)+1]=G[i];}
-    }
-  }
-  if(size(G1)==0){t1=0;}
-  if(size(G2)==0){t2=0;}
-  setring(@RP);
-  if(t1)
-  {
-    list G1RP=imap(R,G1);
-  }
-  else {list G1RP;}
-  list P1RP;
-  ideal B;
-  for(i=1;i<=size(G1RP);i++)
-  {
-     kill B;
-     ideal B;
-     for(k=1;k<=size(G1RP[i][3]);k++)
-     {
-       attrib(G1RP[i][3][k][1],"IsSB",1);
-       G1RP[i][3][k][1]=std(G1RP[i][3][k][1]);
-       for(j=1;j<=size(G1RP[i][2]);j++)
-       {
-         B[j]=reduce(G1RP[i][2][j],G1RP[i][3][k][1]);
-       }
-       P1RP[size(P1RP)+1]=list(G1RP[i][3][k][1],G1RP[i][3][k][2],B);
-     }
-  }
-  setring(R);
-  ideal h;
-  if(t1)
-  {
-    def P1=imap(@RP,P1RP);
-    for(i=1;i<=size(P1);i++)
-    {
-      for(j=1;j<=size(P1[i][3]);j++)
-      {
-        h=factorize(P1[i][3][j],1);
-        P1[i][3][j]=h[1];
-        for(k=2;k<=size(h);k++)
-        {
-          P1[i][3][j]=P1[i][3][j]*h[k];
-        }
-      }
-    }
-  }
-  else{list P1;}
-  ideal BB;
-  for(i=1;i<=size(P1);i++)
-  {
-    if (indepparameters(P1[i][3])==1){P1[i][3]="Special";}
-    else{P1[i][3]="Relevant";}
-  }
-  list P2;
-  for(i=1;i<=size(G2);i++)
-  {
-    for(k=1;k<=size(G2[i][3]);k++)
-    {
-      P2[size(P2)+1]=list(G2[i][3][k][1],G2[i][3][k][2]);
-    }
-  }
-  list l;
-  for(i=1;i<=size(P1);i++){Q1[i]=l; Q1[i][1]=P1[i];} P1=Q1;
-  for(i=1;i<=size(P2);i++){Q2[i]=l; Q2[i][1]=P2[i];} P2=Q2;
-
-  setring @P;
-  ideal J;
-  if(t1==1)
-  {
-    def C1=imap(R,P1);
-    def L1=addcons(C1);
-   }
-  else{list C1; list L1; kill P1; list P1;}
-  if(t2==1)
-  {
-    def C2=imap(R,P2);
-    def L2=addcons(C2);
-  }
-  else{list L2; list C2; kill P2; list P2;}
-    for(i=1;i<=size(L2);i++)
-    {
-      J=std(L2[i][2]);
-      d=dim(J); // AQUI
-      if(d==0)
-      {
-        L2[i][4]=string("Relevant",L2[i][4]);
-      }
-      else{L2[i][4]=string("Degenerate",L2[i][4]);}
-    }
-  list LN;
-  list ll; list l0;
-  if(t1==1)
-  {
-    for(i=1;i<=size(L1);i++)
-    {
-      if(L1[i][4]=="Relevant")
-      {
-       LN[size(LN)+1]=l0;
-       LN[size(LN)][1]=L1[i];
-     }
-    }
-  }
-  if(t2==1)
-  {
-    for(i=1;i<=size(L2);i++)
-    {
-      if(L2[i][4]=="Relevant")
-      {
-        LN[size(LN)+1]=l0;
-        LN[size(LN)][1]=L2[i];
-      }
-    }
-  }
-  list LNN;
-  kill ll; list ll; list lll; list llll;
-  for(i=1;i<=size(LN);i++)
-  {
-      LNN[size(LNN)+1]=l0;
-      ll=LN[i][1]; lll[1]=ll[2]; lll[2]=ll[3]; lll[3]=ll[4]; LNN[size(LNN)][1]=lll;
-  }
-  kill LN; list LN;
-  LN=addcons(LNN);
-  if(size(LN)==0){ojo=1;}
-  setring(R);
-  if(ojo==1){list L;}
-  else
-  {
-    def L=imap(@P,LN);
-    for(i=1;i<=size(L);i++){if(size(L[i][2])==0){L[i][2]=ideal(1);}}
-  }
-  kill @R; kill @RP; kill @P;
-  list LL;
-  for(i=1;i<=size(L);i++)
-  {
-    LL[i]=l;
-    LL[i]=elimfromlist(L[i],1);
-  }
-  return(LL);
-}
-
+    if(typeof(L[i][4])=="list") {L[i][4][1]="Special";}
+    l[1]=L[i][2];
+    l[2]=L[i][3];
+    l[3]=L[i][4];
+    l[4]=L[i][5];
+    L[i]=l;
+ }
+ return(L);
+}
 
 //  locus(G):  Special routine for determining the locus of points
-//                 of  geometrical constructions. Given a parametric ideal J with
-//                 parameters (a_1,..a_m) and variables (x_1,..,xn),
-//                 representing the system determining
-//                 the locus of points (a_1,..,a_m) who verify certain
-//                 properties, computing the grobcov G of
-//                 J and applying to it locus, determines the different
-//                 classes of locus components. The components can be
-//                 Normal, Special, Accumulation, Degenerate.
-//                 The output are the components given in P-canonical form
-//                 of at most 4 constructible sets: Normal, Special, Accumulation,
-//                 Degenerate.
-//                 The description of the algorithm and definitions is
-//                 given in a forthcoming paper by Abanades, Botana, Montes, Recio.
-//                 Usually locus problems have mover coordinates, variables and tracer coordinates.
-//                 The mover coordinates are to be placed as the last variables, and by default,
-//                 its number is 2. If one consider locus problems in higer dimensions, the number of
-//                 mover coordinates (placed as the last variables) is to be given as an option.
-//
+//                 of  geometrical constructions.
 //  input:      The output G of the grobcov (in generic representation, which is the default option for grobcov)
 //  output:
@@ -4397 +4369 @@
 //               gives the depth of the component of the constructible set.
 proc locus(list GG, list #)
-  "USAGE:   locus(G);
-                 The input must be the grobcov  of a parametrical ideal
-  RETURN:  The  locus.
-                 The output are the components of four constructible subsets of the locus
-                 of the parametrical system: Normal , Special, Accumulation and Degenerate.
-                 These components are
-                 given as a list of  (pi,(pi1,..pis_i),type_i,level_i) varying i,
-                 where the p's are prime ideals, the type can be: Normal, Special,
-                 Accumulation, Degenerate.
-  NOTE:      It can only be called after computing the grobcov of the
-                 parametrical ideal in generic representation ('ext',0),
-                 which is the default.
-                 The basering R, must be of the form Q[a_1,..,a_m][x_1,..,x_n].
-  KEYWORDS: geometrical locus, locus, loci.
-  EXAMPLE:  locus; shows an example"
-{
+"USAGE: locus(G)
+        The input must be the grobcov  of a parametrical ideal in Q[a][x],
+        (a=parameters, x=variables). In fact a must be the tracer coordinates
+        and x the mover coordinates and other auxiliary ones.
+        Special routine for determining the locus of points
+        of  geometrical constructions. Given a parametric ideal J
+        representing the system determining the locus of points (a)
+        who verify certain properties, the call to locus on the output of grobcov( J )
+        determines the different classes of locus components, following
+        the taxonomy defined in
+        Abanades, Botana, Montes, Recio:
+        \"An Algebraic Taxonomy for Locus Computation in Dynamic Geometry\".
+        Computer-Aided Design 56 (2014) 22-33.
+        The components can be Normal, Special, Accumulation, Degenerate.
+OPTIONS: The algorithm allows the following options as pair of arguments:
+         \"movdim\", d  : by default movdim is 2 but it can be set to other values,
+         and represents the number of mever variables. they should be given as
+         the last variables of the ring.
+         \"version\", v   :  There are two versions of the algorithm. (\"version\",1) is
+         a full algorithm that always distinguishes correctly between 'Normal'
+         and 'Special' components, whereas \("version\",0) can decalre a component
+         as 'Normal' being really 'Special', but is more effective. By default (\"version\",1)
+         is used when the number of variables is less than 4 and 0 if not.
+         The user can force to use one or other version, but it is not recommended.
+         \"system\", ideal F: if the initial system is passed as an argument. This is actually not used.
+         \"comments\", c: by default it is 0, but it can be set to 1.
+         Usually locus problems have mover coordinates, variables and tracer coordinates.
+         The mover coordinates are to be placed as the last variables, and by default,
+         its number is 2. If one consider locus problems in higer dimensions, the number of
+         mover coordinates (placed as the last variables) is to be given as an option.
+RETURN: The locus. The output is a list of the components ( C_1,.. C_n ) where
+        C_i =  (p_i,(p_i1,..p_is_i), type_i,l evel_i )   and type_i can be
+        'Normal', 'Special', Accumulation', 'Degenerate'. The 'Special' components
+        return more information, namely the antiimage of the component, that
+        is 0-dimensional for these kind of components.
+        Normal components: for each point in the component,
+        the number of solutions in the variables is finite, and
+        the solutions depend on the point in the component.
+        Special components:  for each point in the component,
+        the number of solutions in the variables is finite. The
+        antiimage of the component is 0-dimensional.
+        Accumlation points: are 0-dimensional components whose
+        antiimage is not zero-dimansional.
+        Degenerate components: are components of dimension greater than 0
+        whose antiimage is not-zero-diemansional.
+        The components are given in canonical P-representation.
+        The levels of a class of locus are 1,
+        because they represent locally closed. sets.
+NOTE: It can only be called after computing the grobcov of the
+      parametrical ideal in generic representation ('ext',0),
+      which is the default.
+      The basering R, must be of the form Q[a_1,..,a_m][x_1,..,x_n].
+KEYWORDS: geometrical locus, locus, loci.
+EXAMPLE: locus; shows an example"
+{
+  int tes=0; int i;
   def R=basering;
-  int dd=2; int comment=0;
+  if(defined(@P)==1){tes=1; kill @P; kill @R; kill @RP;}
+  setglobalrings();
+  // Options
   list DD=#;
-  if (size(DD)>1){comment=1;}
-  if(size(DD)>0){dd=DD[1];}
-  int i; int j; int k;
+  int moverdim=nvars(R);
+  int version=0;
+  int nv=nvars(R);
+  if(nv<4){version=1;}
+  int comment=0;
+  ideal Fm;
+  for(i=1;i<=(size(DD) div 2);i++)
+  {
+    if(DD[2*i-1]=="movdim"){moverdim=DD[2*i];}
+    if(DD[2*i-1]=="version"){version=DD[2*i];}
+    if(DD[2*i-1]=="system"){Fm=DD[2*i];}
+    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
+    if(DD[2*i-1]=="family"){poly F=DD[2*i];}
+  }
+  int j; int k;
   def B0=GG[1][2];
-  if (equalideals(B0,ideal(1)))
-  {
-    return(locus0(GG));
+  def H0=GG[1][3][1][1];
+  if (equalideals(B0,ideal(1)) or (equalideals(H0,ideal(0))==0))
+  {
+    return(locus0(GG,DD));
   }
   else
@@ -4449 +4476 @@
     }
     N=px,py;
-    setglobalrings();
     te=indepparameters(N);
     if(te)
@@ -4469 +4495 @@
         if(te==0){nGP[size(nGP)+1]=GP[j];}
       }
-   }
-  }
+     }
+    else{"Warning! Problem with more than one mover. Not able to solve it."; list L; return(L);}
+  }
+  def LL=locus0(nGP,DD);
   kill @RP; kill @P; kill @R;
-  return(locus0(nGP,dd));
+  return(LL);
 }
 example
-{"EXAMPLE:"; echo = 2;
-  ring R=(0,a,b),(x4,x3,x2,x1),dp;
-  ideal S=(x1-3)^2+(x2-1)^2-9,
-             (4-x2)*(x3-3)+(x1-3)*(x4-1),
-             (3-x1)*(x3-x1)+(4-x2)*(x4-x2),
-             (4-x4)*a+(x3-3)*b+3*x4-4*x3,
-             (a-x1)^2+(b-x2)^2-(x1-x3)^2-(x2-x4)^2;
-  short=0;
-  locus(grobcov(S)); " ";
-  kill R;
+{
+  "EXAMPLE:"; echo = 2;
   ring R=(0,a,b),(x,y),dp;
   short=0;
-  "Concoid";
+
+  // Concoid
   ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
-  "System="; S96; " ";
+  // System S96=
+  S96;
   locus(grobcov(S96));
-  kill R; ring R=(0,x,y),(x1,x2),dp;
+  kill R;
+  // ********************************************
+  ring R=(0,a,b),(x4,x3,x2,x1),dp;
+  short=0;
+  ideal S=(x1-3)^2+(x2-1)^2-9,
+               (4-x2)*(x3-3)+(x1-3)*(x4-1),
+               (3-x1)*(x3-x1)+(4-x2)*(x4-x2),
+               (4-x4)*a+(x3-3)*b+3*x4-4*x3,
+               (a-x1)^2+(b-x2)^2-(x1-x3)^2-(x2-x4)^2;
+  // System S=
+  S;
+  locus(grobcov(S));
+  kill R;
+  "********************************************";
+
+  ring R=(0,x,y),(x1,x2),dp;
+  short=0;
   ideal S=-(x - 5)*(x1 - 1) - (x2 - 2)*(y - 2),
-              (x1 - 5)^2 + (x2 - 2)^2 - 4,
-              -2*(x - 5)*(x2 - 2) + 2*(x1 - 5)*(y - 2);
- "System="; S;
-  locus(grobcov(S)); " ";
-  ideal S1=-(x - x1)*(x1 - 4) + (x2 - y)*(x2 - 4),
-                (x1 - 4)^2 + (x2 - 2)^2 - 4,
-                -2*(x1 - 4)*(2*x2 - y - 4) - 2*(x - 2*x1 + 4)*(x2 - 2);
- "System="; S1;
-  locus(grobcov(S1));
+                (x1 - 5)^2 + (x2 - 2)^2 - 4,
+                -2*(x - 5)*(x2 - 2) + 2*(x1 - 5)*(y - 2);
+  // System S=
+  S;
+  locus(grobcov(S));
 }
 
@@ -4523 +4557 @@
 //                 mover coordinates (placed as the last variables) is to be given as an option.
 //
-//  input:      The output G of the grobcov (in generic representation, which is the default option for grobcov)
+//  input:      The output of locus(G);
 //  output:
 //          list, the canonical P-representation of the Normal and Non-Normal locus:
@@ -4544 +4578 @@
 //               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 GG, list #)
-  "USAGE:   locus(G);
-                 The input must be the grobcov  of a parametrical ideal
-  RETURN:  The  locus.
-                 The output are the components of two constructible subsets of the locus
-                 of the parametrical system.: Normal and Non-normal.
-                 These components are
-                 given as a list of  (pi,(pi1,..pis_i),type_i,level_i) varying i,
-                 where the p's are prime ideals, the type can be: Normal, Special,
-                 Accumulation, Degenerate.
-  NOTE:      It can only be called after computing the grobcov of the
-                 parametrical ideal in generic representation ('ext',0),
-                 which is the default.
-                 The basering R, must be of the form Q[a_1,..,a_m][x_1,..,x_n].
-  KEYWORDS: geometrical locus, locus, loci.
-  EXAMPLE:  locusdg; shows an example"
-{
-  def R=basering;
-  int dd=2; int comment=0;
-  list DD=#;
-  if (size(DD)>1){comment=1;}
-  if(size(DD)>0){dd=DD[1];}
-  int i; int j; int k;
-  def B0=GG[1][2];
-  if (equalideals(B0,ideal(1)))
-  {
-    return(locus0dg(GG));
-  }
-  else
-  {
-    int n=nvars(R);
-    ideal vmov=var(n-1),var(n);
-    ideal N;
-    intvec xw; intvec yw;
-    for(i=1;i<=n-1;i++){xw[i]=0;}
-    xw[n]=1;
-    for(i=1;i<=n;i++){yw[i]=0;}
-    yw[n-1]=1;
-    poly px; poly py;
-    int te=1;
-    i=1;
-    while( te and i<=size(B0))
-    {
-      if((deg(B0[i],xw)==1) and (deg(B0[i])==1)){px=B0[i]; te=0;}
-      i++;
-    }
-    i=1; te=1;
-    while( te and i<=size(B0))
-    {
-      if((deg(B0[i],yw)==1) and (deg(B0[i])==1)){py=B0[i]; te=0;}
-      i++;
-    }
-    N=px,py;
-    setglobalrings();
-    te=indepparameters(N);
-    if(te)
-    {
-      string("locus detected that the mover must avoid point (",N,") in order to obtain the correct locus");
-      // eliminates segments of GG where N is contained in the basis
-      list nGP;
-      def GP=GG;
-      ideal BP;
-      for(j=1;j<=size(GP);j++)
-      {
-        te=1; k=1;
-        BP=GP[j][2];
-        while((te==1) and (k<=size(N)))
-        {
-          if(pdivi(N[k],BP)[1]!=0){te=0;}
-          k++;
-        }
-        if(te==0){nGP[size(nGP)+1]=GP[j];}
-      }
-   }
-  }
-  //"T_nGP="; nGP;
-  kill @RP; kill @P; kill @R;
-  return(locus0dg(nGP,dd));
+proc locusdg(list L)
+"USAGE: locusdg(L)   The call must be locusdg(locus(grobcov(S))).
+RETURN: The output is the list of the 'Relevant' components of the locus
+           in Dynamic Geometry:(C1,..,C:m), where
+                 C_i= ( p_i,(p_i1,..p_is_i), 'Relevant', level_i )
+           The 'Relevant' components are the 'Normal' and 'Accumulation' components
+           of the locus. (See help for locus).
+
+NOTE: It can only be called after computing the locus.
+           Calling sequence:    locusdg(locus(grobcov(S)));
+KEYWORDS: geometrical locus, locus, loci, dynamic geometry
+EXAMPLE: locusdg; shows an example"
+{
+  list LL;
+  int i;
+  for(i=1;i<=size(L);i++)
+  {
+    if(typeof(L[i][3])=="string")
+    {
+      if((L[i][3]=="Normal") or (L[i][3]=="Accumulation")){L[i][3]="Relevant"; LL[size(LL)+1]=L[i];}
+    }
+  }
+  return(LL);
 }
 example
-{"EXAMPLE:"; echo = 2;
+{
+  "EXAMPLE:"; echo = 2;
+  ring R=(0,a,b),(x,y),dp;
+  short=0;
+
+  // Concoid
+  ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
+  // System S96=
+  S96;
+  locus(grobcov(S96));
+  locusdg(locus(grobcov(S96)));
+  kill R;
+  "********************************************";
   ring R=(0,a,b),(x4,x3,x2,x1),dp;
   ideal S=(x1-3)^2+(x2-1)^2-9,
-             (4-x2)*(x3-3)+(x1-3)*(x4-1),
-             (3-x1)*(x3-x1)+(4-x2)*(x4-x2),
-             (4-x4)*a+(x3-3)*b+3*x4-4*x3,
-             (a-x1)^2+(b-x2)^2-(x1-x3)^2-(x2-x4)^2;
+               (4-x2)*(x3-3)+(x1-3)*(x4-1),
+               (3-x1)*(x3-x1)+(4-x2)*(x4-x2),
+               (4-x4)*a+(x3-3)*b+3*x4-4*x3,
+               (a-x1)^2+(b-x2)^2-(x1-x3)^2-(x2-x4)^2;
   short=0;
-  locus(grobcov(S)); " ";
+  locus(grobcov(S));
+  locusdg(locus(grobcov(S)));
   kill R;
-  ring R=(0,a,b),(x,y),dp;
+  "********************************************";
+
+  ring R=(0,x,y),(x1,x2),dp;
   short=0;
-  "Concoid";
-  ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
-  "System="; S96; " ";
-  locusdg(grobcov(S96));
-  kill R; ring R=(0,x,y),(x1,x2),dp;
   ideal S=-(x - 5)*(x1 - 1) - (x2 - 2)*(y - 2),
-              (x1 - 5)^2 + (x2 - 2)^2 - 4,
-              -2*(x - 5)*(x2 - 2) + 2*(x1 - 5)*(y - 2);
- "System="; S;
-  locusdg(grobcov(S)); " ";
-  ideal S1=-(x - x1)*(x1 - 4) + (x2 - y)*(x2 - 4),
-                (x1 - 4)^2 + (x2 - 2)^2 - 4,
-                -2*(x1 - 4)*(2*x2 - y - 4) - 2*(x - 2*x1 + 4)*(x2 - 2);
- "System="; S1;
-  locusdg(grobcov(S1));
+                (x1 - 5)^2 + (x2 - 2)^2 - 4,
+                -2*(x - 5)*(x2 - 2) + 2*(x1 - 5)*(y - 2);
+  locus(grobcov(S));
+  locusdg(locus(grobcov(S)));
 }
 
@@ -4660 +4644 @@
 //     string s: The output of locus converted to a string readable by other programs
 proc locusto(list L)
-"USAGE:   locusto(L);
-          The argument must be the output of locus of a parametrical ideal
+"USAGE: locusto(L);
+          The argument must be the output of locus or locusdg or
+          envelop or envelopdg.
           It transforms the output into a string in standard form
           readable in many languages (Geogebra).
 RETURN: The locus in string standard form
-NOTE:     It can only be called after computing the locus(grobcov(F)) of the
-              parametrical ideal.
-              The basering R, must be of the form Q[a,b,..][x,y,..].
+NOTE: It can only be called after computing either
+              - locus(grobcov(F))                -> locusto( locus(grobcov(F)) )
+              - locusdg(locus(grobcov(F)))  -> locusto( locusdg(locus(grobcov(F))) )
+              - envelop(F,C)                       -> locusto( envelop(F,C) )
+              -envelopdg(envelop(F,C))      -> locusto( envelopdg(envelop(F,C)) )
 KEYWORDS: geometrical locus, locus, loci.
 EXAMPLE:  locusto; shows an example"
@@ -4699 +4686 @@
     if(size(L[i])>=3)
     {
-      s[size(s)]=",";
-      s=string(s,string(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)
@@ -4714 +4713 @@
 }
 example
-{"EXAMPLE:"; echo = 2;
-  ring R=(0,a,b),(x,y),dp;
+{
+  "EXAMPLE:"; echo = 2;
+  ring R=(0,x,y),(x1,y1),dp;
   short=0;
-  ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
-  "System="; S96; " ";
-  locusto(locus(grobcov(S96)));
-}
-
-// Auxiliary routione
+  ideal S=x1^2+y1^2-4,(y-2)*x1-x*y1+2*x,(x-x1)^2+(y-y1)^2-1;
+  locusto(locus(grobcov(S)));
+  locusto(locusdg(locus(grobcov(S))));
+  kill R;
+  "********************************************";
+
+  // 1. Take a fixed line l: x1-y1=0  and consider
+  //    the family F of a lines parallel to l passing through the mover point M
+  // 2. Consider a circle x1^2+x2^2-25, and a mover point M(x1,x2) on it.
+  // 3. Compute the envelop of the family of lines.
+
+  ring R=(0,x,y),(x1,y1),lp;
+  poly F=(y-y1)-(x-x1);
+  ideal C=x1^2+y1^2-25;
+  short=0;
+
+  // Curves Family F=
+  F;
+  // Conditions C=
+  C;
+
+  locusto(envelop(F,C));
+  locusto(envelopdg(envelop(F,C)));
+  kill R;
+  "********************************************";
+
+  // Steiner Deltoid
+  // 1. Consider the circle x1^2+y1^2-1=0, and a mover point M(x1,y1) on it.
+  // 2. Consider the triangle A(0,1), B(-1,0), C(1,0).
+  // 3. Consider lines passing through M perpendicular to two sides of ABC triangle.
+  // 4. Obtain the envelop of the lines above.
+
+  ring R=(0,x,y),(x1,y1,x2,y2),lp;
+  ideal C=(x1)^2+(y1)^2-1,
+               x2+y2-1,
+               x2-y2-x1+y1;
+  matrix M[3][3]=x,y,1,x2,y2,1,x1,0,1;
+  poly F=det(M);
+
+  short=0;
+
+  // Curves Family F=
+  F;
+  // Conditions C=
+  C;
+
+  locusto(envelop(F,C));
+  locusto(envelopdg(envelop(F,C)));
+}
+
+// Auxiliary routine
 // locusdim
 // input:
 //    B:         ideal, a basis of a segment of the grobcov
-//    dgdim: int, the dimension of the mover (for locusdg)
+//    dgdim: int, the dimension of the mover (for locus)
 //                by default dgdim is equal to the number of variables
-proc locusdim(ideal B, list #)
+static proc locusdim(ideal B, list #)
 {
   def R=basering;
@@ -4752 +4797 @@
     }
   }
-  //"T_B0="; B0;
-  //"T_v="; v;
   def RR=ringlist(R);
   def RR0=RR;
   RR0[2]=v;
   def R0=ring(RR0);
-  //ringlist(R0);
   setring(R0);
   def B0r=imap(R,B0);
@@ -4767 +4809 @@
 }
 
-// locusdgto: Transforms the output of locus to a string that
-//      can be read by different computational systems.
-// input:
-//     list L: The output of locus
-// output:
-//     string s: The output of locus converted to a string readable by other programs
-//                   Outputs only the relevant dynamical geometry components.
-//                   Without unnecessary parenthesis
-proc locusdgto(list LL)
+// // locusdgto: Transforms the output of locusdg to a string that
+// //      can be read by different computational systems.
+// // input:
+// //     list L: The output of locus
+// // output:
+// //     string s: The output of locus converted to a string readable by other programs
+// //                   Outputs only the relevant dynamical geometry components.
+// //                   Without unnecessary parenthesis
+// proc locusdgto(list LL)
+// "USAGE: locusdgto(L);
+//           The argument must be the output of locusdg of a parametrical ideal
+//           It transforms the output into a string in standard form
+//           readable in many languages (Geogebra).
+// RETURN: The locusdg in string standard form
+// NOTE: It can only be called after computing the locusdg(grobcov(F)) of the
+//           parametrical ideal.
+//           The basering R, must be of the form Q[a,b,..][x,y,..].
+// KEYWORDS: geometrical locus, locus, loci.
+// EXAMPLE:  locusdgto; shows an example"
+// {
+//   int i; int j; int k; int short0=short; int ojo;
+//   int te=0;
+//   short=0;
+//   if(size(LL)==0){ojo=1; list L;}
+//   else
+//   {
+//     def L=LL;
+//   }
+//   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++)
+//   {
+//     if(L[i][3]=="Relevant")
+//     {
+//       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,"]");
+//       s[size(s)]="]";
+//       s=string(s,",");
+//     }
+//   }
+//   if(s=="["){s="[]";}
+//   else{s[size(s)]="]";}
+//   short=short0;
+//   return(s);
+// }
+// example
+// {"EXAMPLE:"; echo = 2;
+//   ring R=(0,a,b),(x,y),dp;
+//   short=0;
+//   ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
+//   "System="; S96; " ";
+//   "locusdgto(locusdg(grobcov(S96)))=";
+//   locusdgto(locusdg(grobcov(S96)));
+// }
+
+static proc norspec(ideal F)
 {
   def RR=basering;
-  int i; int j; int k; int short0=short; int ojo;
-  int te=0;
-  if(defined(@P)){te=1;}
-  else{setglobalrings();}
-  setring @P;
+  def Rx=ringlist(RR);
+
+   def Rx=ringlist(RR);
+  def @P=ring(Rx[1]);
+  list Lx;
+  Lx[1]=0;
+  Lx[2]=Rx[2]+Rx[1][2];
+  Lx[3]=Rx[1][3];
+  Lx[4]=Rx[1][4];
+  Rx[1]=0;
+  def D=ring(Rx);
+  def @RP=D+@P;
+  exportto(Top,@R);      // global ring Q[a][x]
+  exportto(Top,@P);      // global ring Q[a]
+  exportto(Top,@RP);     // global ring K[x,a] with product order
+  setring(RR);
+
+}
+
+// envelop
+// Input: n arguments
+//   poly F: the polynomial defining the family of curves in ring R=0,(x,y),(x1,..,xn),lp;
+//   ideal C=g1,..,g_{n-1}:  the set of constraints
+// Output: the components of the envolvent;
+proc envelop(poly F, ideal C, list #)
+"USAGE: envelop(F,C);
+          The first argument F must be the family of curves of which
+          on want to compute the envelop.
+          The second argument C must be the ideal of conditions
+          over the variables, and should contain as polynomials
+          as the number of variables -1.
+RETURN: The components of the envelop with its taxonomy:
+           The taxonomy distinguishes 'Normal',
+          'Special', 'Accumulation', 'Degenerate' components.
+          In the case of 'Special' components, it also
+          outputs the antiimage of the component
+          and an integer (0-1). If the integer is 0
+          the component is not a curve of the family and is
+          not considered as 'Relevant' by the envelopdg routine
+          applied to it, but is considered as 'Relevant' if the integer is 1.
+NOTE: grobcov is called internally.
+          The basering R, must be of the form Q[a][x] (a=parameters, x=variables).
+KEYWORDS: geometrical locus, locus, loci, envelop
+EXAMPLE:  envelop; shows an example"
+{
+  int tes=0; int i;   int j;
+  def R=basering;
+  if(defined(@P)==1){tes=1; kill @P; kill @R; kill @RP;}
+  setglobalrings();
+  // Options
+  list DD=#;
+  int moverdim=nvars(R);
+  int version=0;
+  int nv=nvars(R);
+  if(nv<4){version=1;}
+  int comment=0;
+  ideal Fm;
+  for(i=1;i<=(size(DD) div 2);i++)
+  {
+    if(DD[2*i-1]=="movdim"){moverdim=DD[2*i];}
+    if(DD[2*i-1]=="version"){version=DD[2*i];}
+    if(DD[2*i-1]=="system"){Fm=DD[2*i];}
+    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
+  }
+  int n=nvars(R);
+  list v;
+  for(i=1;i<=n;i++){v[size(v)+1]=var(i);}
+  def MF=jacob(F);
+  def TMF=transpose(MF);
+  def Mg=MF;
+  def TMg=TMF;
+  for(i=1;i<=n-1;i++)
+  {
+    Mg=jacob(C[i]);
+    TMg=transpose(Mg);
+    TMF=concat(TMF,TMg);
+  }
+  poly J=det(TMF);
+  ideal S=ideal(F)+C+ideal(J);
+  DD[size(DD)+1]="family";
+  DD[size(DD)+1]=F;
+  def G=grobcov(S,DD);
+   def L=locus(G, DD);
+  return(L);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  // Steiner Deltoid
+  // 1. Consider the circle x1^2+y1^2-1=0, and a mover point M(x1,y1) on it.
+  // 2. Consider the triangle A(0,1), B(-1,0), C(1,0).
+  // 3. Consider lines passing through M perpendicular to two sides of ABC triangle.
+  // 4. Obtain the envelop of the lines above.
+
+  ring R=(0,x,y),(x1,y1,x2,y2),lp;
+  ideal C=(x1)^2+(y1)^2-1,
+               x2+y2-1,
+               x2-y2-x1+y1;
+  matrix M[3][3]=x,y,1,x2,y2,1,x1,0,1;
+  poly F=det(M);
+
   short=0;
-  if(size(LL)==0){ojo=1; list L;}
-  else
-  {
-    def L=imap(RR,LL);
-  }
-  string s="["; string sf="]"; string st=s+sf;
-  if(size(L)==0){return(st);}
-  ideal p;
-  ideal q;
+
+  // Curves Family F=
+  F;
+  // Conditions C=
+  C;
+
+  def Env=envelop(F,C);
+  Env;
+}
+
+// envelopdg
+// Input: list L: the output of envelop(poly F, ideal C, list #)
+// Output: the relevant components of the envolvent in dynamic geometry;
+proc envelopdg(list L)
+"USAGE: envelopdg(L);
+          The input list L must be the output of the call to
+          the routine 'envolop' of the family of curves
+RETURN: The relevant components of the envelop in Dynamic Geometry.
+           'Normal' and 'Accumulation' components are always considered
+           'Relevant'. 'Special' components of the envelop outputs
+           three objects in its characterization: 'Special', the antiimage ideal,
+           and the integer 0 or 1, that indicates that the given component is
+           formed (1) or is not formed (0) by curves of the family. Only if yes,
+           'envelopdg' considers the component as 'Relevant' .
+NOTE: It must be called to the output of the 'envelop' routine.
+          The basering R, must be of the form Q[a,b,..][x,y,..].
+KEYWORDS: geometrical locus, locus, loci, envelop.
+EXAMPLE:  envelop; shows an example"
+{
+  list LL;
+  list Li;
+  int i;
   for(i=1;i<=size(L);i++)
   {
-    if(L[i][3]=="Relevant")
-    {
-      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++)
+    if(typeof(L[i][3])=="string")
+    {
+      if((L[i][3]=="Normal") or (L[i][3]=="Accumulation")){Li=L[i]; Li[3]="Relevant"; LL[size(LL)+1]=Li;}
+    }
+    else
+    {
+      if(typeof(L[i][3])=="list")
+      {
+        if(L[i][3][3]==1)
         {
-          s=string(s,L[i][2][j][k],",");
+          Li=L[i]; Li[3]="Relevant"; LL[size(LL)+1]=Li;
         }
-        s[size(s)]="]";
-        s=string(s,",");
-      }
-      s[size(s)]="]";
-      s=string(s,"]");
-      s[size(s)]="]";
-      s=string(s,",");
-    }
-  }
-  if(s=="["){s="[]";}
-  else{s[size(s)]="]";}
-  setring(RR);
-  short=short0;
-  if(te==0){kill @P; kill @R; kill @RP;}
-  return(s);
+      }
+    }
+  }
+  return(LL);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+
+  // 1. Take a fixed line l: x1-y1=0  and consider
+  //    the family F of a lines parallel to l passing through the mover point M
+  // 2. Consider a circle x1^2+x2^2-25, and a mover point M(x1,x2) on it.
+  // 3. Compute the envelop of the family of lines.
+
+  ring R=(0,x,y),(x1,y1),lp;
+  short=0;
+  poly F=(y-y1)-(x-x1);
+  ideal C=x1^2+y1^2-25;
+  short=0;
+
+  // Curves Family F=
+  F;
+  // Conditions C=
+  C;
+
+  envelop(F,C);
+  envelopdg(envelop(F,C));
 }
 
 //********************* End locus ****************************
+
+//********************* Begin AddCons **********************
+
+// Input: L1,L2: lists of components with common top
+// Output L: list of the union of L1 and L2.
+static proc Add2ComWithCommonTop(list L1, list L2)
+{
+  int i; int j; ideal pij; list L; list Lp; list PR; int k;
+  for(i=1;i<=size(L1[2]);i++)
+  {
+    for(j=1;j<=size(L2[2]);j++)
+    {
+      pij=std(L1[2][i]+L2[2][j]);
+      PR=minGTZ(pij);
+      for(k=1;k<=size(PR);k++)
+      {
+        Lp[size(Lp)+1]=PR[k];
+      }
+    }
+  }
+  for(i=1; i<=size(Lp);i++)
+  {
+    for(j=i+1;j<=size(Lp);j++)
+    {
+      if(idcontains(Lp[i],Lp[j])) {Lp=delete(Lp,j);}
+    }
+    for(j=1;j<i;j++)
+    {
+      if(idcontains(Lp[i],Lp[j])){Lp=delete(Lp,j); j=j-1; i=i-1;}
+    }
+  }
+  L[1]=L1[1];
+  L[2]=Lp;
+  return(L);
+}
+
+// Input: L list od P-rep of components to be added. L[i]=[p_i,[p_{i1},...p_{ir_i}]]
+// Output: lists A,B,L
+//       where no top in the lists are repeated
+//       A: list of components with higher tops
+//       B: list of components with lower tops
+//       L1: A,B
+static proc SepareAB(list L)
+{
+  int i;  int j; list Ln=L;
+  for(i=1;i<=size(Ln);i++)
+  {
+    for(j=i+1;j<=size(Ln);j++)
+    {
+      if (equalideals(Ln[j][1],Ln[i][1]))
+      {
+        Ln[i]=Add2ComWithCommonTop(Ln[i],Ln[j]);
+        Ln=delete(Ln,j);
+        j=j-1;
+      }
+    }
+  }
+  list A; list B; int clas; list T; list L1;
+  for(i=1;i<=size(Ln);i++)
+  {
+    j=1;
+    clas=0;
+    while((clas==0) and  (j<=size(Ln)))
+    {
+      if(j!=i)
+      {
+        if(idcontains(Ln[i][1],Ln[j][1]) ) {B[size(B)+1]=Ln[i]; clas=1;}
+      }
+      j++;
+    }
+    if(clas==0) {A[size(A)+1]=Ln[i];}
+  }
+  L1=A; for(j=1;j<=size(B);j++){L1[size(L1)+1]=B[j];}
+  T[1]=A; T[2]=B; T[3]=L1;
+  return(T);
+}
+
+// Input:
+//  A1: list of high set of P-reps to be completed by the remaining P-reps
+//  L1: the list A1, completed with the list of lower P-reps.
+// Output:
+//  A: list A1 completed with all possible parts of the remaining parts of L1 with the
+//      condition of building locally closed sets.
+static proc CompleteA(list A1,list L1)
+{
+  int modif; int i; int j; int k; int l;
+  ideal pij; ideal pk; ideal pijkl; ideal pkl;
+  int n; list nl; int te; int clas; list vvv; list AAA;
+  list Lp; int m; ideal Pij;
+  list A0;
+  modif=1;
+  list A=A1;
+  while(modif==1)
+  {
+      modif=0;
+      A0=A;
+      for(i=1;i<=size(A);i++)
+      {
+          for(j=1;j<=size(A[i][2]); j++)
+          {
+              pij=A[i][2][j];
+             for(k=1;k<=size(L1);k++)
+             {
+                 if(k!=i)
+                 {
+                     pk=L1[k][1];
+                     if(idcontains(pij,pk)==1)
+                     {
+                         te=0;
+                         kill nl;
+                         list nl; l=1;
+                         while((te==0) and (l<=size(L1[k][2])))
+                         {
+                              pkl=L1[k][2][l];
+                              if((equalideals(pij,pkl)==1) or (idcontains(pij,pkl)==1)) {te=1;}
+                              l++;
+                         }   // end while ((te=0) and (l>...
+                         //"T_te="; te; pause();
+                         if(te==0)
+                         {
+                           modif=1;
+                           kill Pij; ideal Pij=1;
+                           for(l=1; l<=size(L1[k][2]);l++)
+                           {
+                             pkl=L1[k][2][l];
+                             pijkl=std(pij+pkl);
+                             Pij=intersect(Pij,pijkl);
+                           }
+                           kill Lp; list Lp;
+                           Lp=minGTZ(Pij);
+                           for(m=1;m<=size(Lp);m++)
+                            {
+                               nl[size(nl)+1]=Lp[m];
+                            }
+                            A[i][2]=delete(A[i][2], j);
+                            for(n=1;n<=size(nl);n++)
+                            {
+                              A[i][2][size(A[i][2])+1]=nl[n];
+                            }
+                          } // end if(te==0)
+                     } // end if(idcontains(pij,pk)==1)
+                 }  // end if (k!=i)
+             } // end for k
+         }  // end for j
+         kill vvv; list vvv;
+         if(modif==1)
+         // Select the maximal ideals of the set A[I][2][j]
+         {
+             kill nl; list nl;
+             nl=selectminideals(A[i][2]);
+             kill AAA; list AAA;
+             for(j=1;j<=size(nl);j++)
+             {
+               AAA[size(AAA)+1]=A[i][2][nl[j]];
+             }
+             A[i][2]=AAA;
+         } // end if(modif=1)
+      }  // end for i
+      modif=1-equallistsAall(A,A0);
+  } // end while(modif==1)
+  return(A);
+}
+
+// Input:
+//   A: list of the top P-reps of one level
+//   B: list of remaining lower P-reps that have not yeen be possible to add
+// Output:
+//   Bn: list B where the elements that are completely included in A are removed and the parts that are
+//         included in A also.
+static proc ReduceB(list A,list B)
+{
+     int i; int j; list nl; list Bn; int te; int k; int elim;
+     ideal pC; ideal pD; ideal pCD; ideal pBC; list nB; int j0;
+     for(i=1;i<=size(B);i++)
+     {
+         j=1; te=0;
+         while((te==0) and (j<=size(A)))
+         {
+             if(idcontains(B[i][1],A[j][1])){te=1; j0=j;}
+             else{j++;}
+         }
+         pD=B[i][2][1];
+         for(k=2;k<=size(B[i][2]);k++){pD=intersect(pD,B[i][2][k]);}
+         pC=A[j0][2][1];
+         for(k=2;k<=size(A[j0][2]);k++) {pC=intersect(pC,A[j0][2][k]);}
+         pCD=std(pD+pC);
+         pBC=std(B[i][1]+pC);
+         elim=0;
+         if(idcontains(pBC,pCD)){elim=1;}   // B=delfromlist(B,i);i=i-1;
+         else
+         {
+              nB=Prep(pBC,pCD);
+              if(equalideals(nB[1][1],ideal(1))==0)
+              {
+                  for(k=1;k<=size(nB);k++){Bn[size(Bn)+1]=nB[k];}
+              }
+         }
+    }   // end for i
+    return(Bn);
+}
+
+// AddConsP: 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.
+proc AddConsP(list L)
+"USAGE:   AddConsP(L)
+      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]]  ]
+RETURN:
+     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.
+NOTE:     Operates in a ring R=Q[u_1,..,u_m]
+KEYWORDS: Constructible set, Canoncial form
+EXAMPLE:  AddConsP; shows an example"
+{
+  list LL; list A; list B; list L1; list T; int level=0; list h;
+  LL=L; int i;
+  while(size(LL)!=0)
+  {
+    level++;
+    L1=SepareAB(LL);
+    A=L1[1]; B=L1[2]; LL=L1[3];
+    A=CompleteA(A,LL);
+    for(i=1;i<=size(A);i++)
+    {
+      LL[i]=A[i];
+    }
+    h[1]=level; h[2]=A;
+    T[size(T)+1]=h;
+    LL=ReduceB(A,B);
+  }
+  return(T);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  if (defined(Grobcov::@P)){kill Grobcov::@P; kill Grobcov::@R; kill Grobcov::@RP;}
+  ring R=0,(x,y,z),lp;
+  short=0;
+
+  ideal P1=x;
+  ideal Q1=x,y;
+  ideal P2=y;
+  ideal Q2=y,z;
+
+  list L=list(Prep(P1,Q1)[1],Prep(P2,Q2)[1]);
+  L;
+  AddConsP(L);
+}
+
+// AddCons:  given a set of  locally closed sets by pairs of ideal, it builds the
+//       canonical P-representation of the corresponding constructible set, of its union,
+//       including levels it they are.
+//       It makes a call to AddConsP after transforming the input.
+// Input list of lists of pairs of ideals representing locally colsed sets:
+//     L=  [ [E1,N1], .. , [E_s,N_s] ]
+// Output: The canonical frepresentation of the constructible set union of the V(E_i) \ V(N_i)
+//     T=[  [level_1,[ [p_1,[p_11,..,p_1s_1]],.., [p_k,[p_k1,..,p_ks_k]] ]],, .. , [level_r,[..       ]]  ]
+proc AddCons(list L)
+"USAGE:   AddCons(L)
+      Calls internally AddConsP and allows a different input.
+      Input L: list of pairs of of ideals [ [P_1,Q_1], .., [Pr,Qr] ]
+      representing a set of locally closed setsV(P_i) \ V(Q_i)
+      to be added.
+RETURN:
+      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.
+NOTE: Operates in a ring R=Q[u_1,..,u_m]
+KEYWORDS: Constructible set, Canoncial form
+EXAMPLE: AddCons; shows an example"
+{
+  int i; list H; list P; int j;
+  for(i=1;i<=size(L);i++)
+  {
+    P=Prep(L[i][1],L[i][2]);
+    for(j=1;j<=size(P);j++)
+    {
+      H[size(H)+1]=P[j];
+    }
+  }
+  return(AddConsP(H));
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  if (defined(Grobcov::@P)){kill Grobcov::@P; kill Grobcov::@R; kill Grobcov::@RP;}
+  ring R=0,(x,y,z),lp;
+  short=0;
+
+  ideal P1=x;
+  ideal Q1=x,y;
+  ideal P2=y;
+  ideal Q2=y,z;
+
+  list L=list(list(P1,Q1), list(P2,Q2) );
+  L;
+
+  AddCons(L);
+}
+
+// AddLocus: auxilliary routine for locus0 that computes the components of the constructible:
+// Input:  the list of locally closed sets to be added, each with its type as third argument
+//     L=[ [LC[11],,,LC[1k_1],  .., [LC[r1],,,LC[rk_r] ] where
+//            LC[1]=[p1,[p11,..,p1k],typ]
+// Output:  the list of components of the constructible union of the L, with the type of the corresponding top
+//               and the level of the constructible
+//     L4= [[v1,p1,[p11,..,p1l],typ_1,level]_1 ,.. [vs,ps,[ps1,..,psl],typ_s,level_s]
+static proc AddLocus(list L)
+{
+//  int te0=0;
+//  def RR=basering;
+//  if(defined(@P)){te0=1;  def Rx=@R;  kill @P; setring RR;}
+  list L1; int i; int j;  list L2; list L3;
+  list l1; list l2;
+  intvec v;
+  for(i=1; i<=size(L); i++)
+  {
+    for(j=1;j<=size(L[i]);j++)
+    {
+      l1[1]=L[i][j][1];
+      l1[2]=L[i][j][2];
+      l2[1]=l1[1];
+      if(size(L[i][j])>2){l2[3]=L[i][j][3];}
+      v[1]=i; v[2]=j;
+      l2[2]=v;
+      L1[size(L1)+1]=l1;
+      L2[size(L2)+1]=l2;
+    }
+  }
+  L3=AddConsP(L1);
+  list L4; int level;
+  ideal p1; ideal pp1; int t; int k; int k0; string typ; list l4;
+  for(i=1;i<=size(L3);i++)
+  {
+    level=L3[i][1];
+    for(j=1;j<=size(L3[i][2]);j++)
+    {
+      p1=L3[i][2][j][1];
+      t=1; k=1;
+      while((t==1) and (k<=size(L2)))
+      {
+        pp1=L2[k][1];
+        if(equalideals(p1,pp1)){t=0; k0=k;}
+        k++;
+      }
+      if(t==0)
+      {
+        v=L2[k0][2];
+      }
+      else{"ERROR p1 NOT FOUND";}
+      l4[1]=v; l4[2]=p1; l4[3]=L3[i][2][j][2];  l4[5]=level;
+      if(size(L2[k0])>2){l4[4]=L2[k0][3];}
+      L4[size(L4)+1]=l4;
+    }
+  }
+  return(L4);
+}
+
+//********************* End AddCons **********************
 ;

Singular/LIB/grwalk.lib

@@ -255 +255 @@
 }
 
-proc gwalk(ideal Go, list #)
-"SYNTAX: gwalk(ideal i);
-         gwalk(ideal i, intvec v, intvec w);
+proc gwalk(ideal Go, int reduction,int printout, list #)
+"SYNTAX: gwalk(ideal i, int reduction, int printout);
+         gwalk(ideal i, int reduction, int printout, intvec v, intvec w);
 TYPE:    ideal
 PURPOSE: compute the standard basis of the ideal, calculated via
@@ -284 +284 @@
    //print("//** help ring = " + string(basering));
    ideal G = fetch(xR, Go);
-   G = system("Mwalk", G, curr_weight, target_weight,basering);
+   G = system("Mwalk", G, curr_weight, target_weight,basering,reduction,printout);
 
    setring xR;
@@ -300 +300 @@
   ring r = 32003,(z,y,x), lp;
   ideal I = y3+xyz+y2z+xz3, 3+xy+x2y+y2z;
-  gwalk(I);
+  gwalk(I,0,1);
 }
 
@@ -346 +346 @@
 }
 
-proc fwalk(ideal Go, list #)
-"SYNTAX: fwalk(ideal i);
-         fwalk(ideal i, intvec v, intvec w);
+proc fwalk(ideal Go, int reduction, int printout, list #)
+"SYNTAX: fwalk(ideal i,int reductioin);
+         fwalk(ideal i, int reduction intvec v, intvec w);
 TYPE:    ideal
 PURPOSE: compute the standard basis of the ideal w.r.t. the
@@ -372 +372 @@
 
    ideal G = fetch(xR, Go);
-   G = system("Mfwalk", G, curr_weight, target_weight);
+   G = system("Mfwalk", G, curr_weight, target_weight, reduction, printout);
 
    setring xR;
@@ -387 +387 @@
     ring r = 32003,(z,y,x), lp;
     ideal I = y3+xyz+y2z+xz3, 3+xy+x2y+y2z;
-    fwalk(I);
+    int reduction = 1;
+    int printout = 1;
+    fwalk(I,reduction,printout);
 }
 
@@ -437 +439 @@
 }
 
-proc pwalk(ideal Go, int n1, int n2, list #)
-"SYNTAX: pwalk(int d, ideal i, int n1, int n2);
-         pwalk(int d, ideal i, int n1, int n2, intvec v, intvec w);
+proc pwalk(ideal Go, int n1, int n2, int reduction, int printout, list #)
+"SYNTAX: pwalk(int d, ideal i, int n1, int n2, int reduction, int printout);
+         pwalk(int d, ideal i, int n1, int n2, int reduction, int printout, intvec v, intvec w);
 TYPE:    ideal
 PURPOSE: compute the standard basis of the ideal, calculated via
@@ -477 +479 @@
   ideal G = fetch(xR, Go);
 
+<<<<<<< HEAD
+  G = system("Mpwalk",G,n1,n2,curr_weight,target_weight,nP,reduction,printout);
+ 
+=======
   G = system("Mpwalk", G, n1, n2, curr_weight, target_weight,nP);
 
+>>>>>>> f533f6f7667328bccb271b19b2f603aaebe41596
   setring xR;
   //kill Go;
@@ -492 +499 @@
     ring r = 32003,(z,y,x), lp;
     ideal I = y3+xyz+y2z+xz3, 3+xy+x2y+y2z;
-    //I = std(I);
-    //ring rr = 32003,(z,y,x),lp;
-    //ideal I = fetch(r,I);
-    pwalk(I,2,2);
+    int reduction = 1;
+    int printout = 2;
+    pwalk(I,2,2,reduction,printout);
 }
 

Singular/LIB/homolog.lib

@@ -534 +534 @@
   { // compute resolution via sres command
     if( attrib(M,"isSB")==0 ) {
-      if (size(M)==0) { attrib(M,"isSB")=1; }
-      else { M=std(M); }
+      M=std(M);
     }
     list resl = sres(M,max+1);

Singular/LIB/modwalk.lib

@@ -16 +16 @@
 
 PROCEDURES:
+modpWalk(ideal,int,int,int[,int,int,int,int])    standard basis conversion of I in prime characteristic
 modWalk(ideal,int,int[,int,int,int,int]);        standard basis conversion of I using modular methods (chinese remainder)
 ";
@@ -29 +30 @@
 ////////////////////////////////////////////////////////////////////////////////
 
-static proc modpWalk(def II, int p, int variant, list #)
+proc modpWalk(def II, int p, int variant, int reduction, list #)
 "USAGE:  modpWalk(I,p,#); I ideal, p integer, variant integer
 ASSUME:  If size(#) > 0, then
@@ -80 +81 @@
     }
   }
+<<<<<<< HEAD
+=======
 
 //-------------------------  make i homogeneous  -----------------------------
@@ -106 +109 @@
   }
 
+>>>>>>> f533f6f7667328bccb271b19b2f603aaebe41596
 //-------------------------  compute a standard basis mod p  -----------------------------
-
   if(variant == 1)
   {
     if(size(#)>0)
     {
-      i = rwalk(i,radius,pert_deg,#);
-     // rwalk(i,radius,pert_deg,#); std(i);
+      i = rwalk(i,radius,pert_deg,reduction,#);
     }
     else
     {
-      i = rwalk(i,radius,pert_deg);
+      i = rwalk(i,radius,pert_deg,reduction);
     }
   }
@@ -124 +126 @@
     if(size(#) == 2)
     {
-      i = gwalk(i,#);
+      i = gwalk(i,reduction,#);
     }
     else
     {
-      i = gwalk(i);
+      i = gwalk(i,reduction);
     }
   }
@@ -135 +137 @@
     if(size(#) == 2)
     {
-      i = frandwalk(i,radius,#);
+      i = frandwalk(i,radius,reduction,#);
     }
     else
     {
-      i = frandwalk(i,radius);
+      i = frandwalk(i,radius,reduction);
     }
   }
@@ -157 +159 @@
     if(size(#) == 2)
     {
-     i=prwalk(i,radius,pert_deg,pert_deg,#);
+     i=prwalk(i,radius,pert_deg,pert_deg,reduction,#);
     }
     else
     {
-      i=prwalk(i,radius,pert_deg,pert_deg);
+      i=prwalk(i,radius,pert_deg,pert_deg,reduction);
     }
   }
@@ -168 +170 @@
     if(size(#) == 2)
     {
-      i=pwalk(i,pert_deg,pert_deg,#);
+      i=pwalk(i,pert_deg,pert_deg,reduction,#);
     }
     else
     {
+<<<<<<< HEAD
+      i=pwalk(i,pert_deg,pert_deg,reduction);
+=======
       i=pwalk(i,pert_deg,pert_deg);
     }
@@ -189 +194 @@
         kill HomR;
       }
-    }
-  }
+>>>>>>> f533f6f7667328bccb271b19b2f603aaebe41596
+    }
+  }
+
   setring R0;
   return(list(fetch(@r,i),p));
@@ -204 +211 @@
   ring ra = 0,x(1..4),(a(a),lp);
   ideal I = std(cyclic(4));
+  int reduction = 1;
   ring rb = 0,x(1..4),(a(b),lp);
   ideal I = imap(ra,I);
-  modpWalk(I,p,1,a,b);
+  modpWalk(I,p,1,reduction,a,b);
   std(I);
 }
@@ -212 +220 @@
 ////////////////////////////////////////////////////////////////////////////////
 
-proc modWalk(def II, int variant, list #)
+proc modWalk(def II, int variant, int reduction, list #)
 "USAGE:  modWalk(II); II ideal or list(ideal,int)
 ASSUME:  If variant =
@@ -487 +495 @@
   if(n2 > 4)
   {
-  //  L[5] = prime(random(an,en));
+    L[5] = prime(random(an,en));
   }
   if(printlevel >= 10)
@@ -504 +512 @@
     for(i=1; i<=size(L); i++)
     {
-      Arguments[i] = list(II,L[i],variant,list(curr_weight,target_weight));
+      Arguments[i] = list(II,L[i],variant,reduction,list(curr_weight,target_weight));
     }
   }
@@ -511 +519 @@
     for(i=1; i<=size(L); i++)
     {
-      Arguments[i] = list(II,L[i],variant);
+      Arguments[i] = list(II,L[i],variant,reduction);
     }
   }
@@ -528 +536 @@
 //-------------------  Now all leading ideals are the same  --------------------
 //-------------------  Lift results to basering via farey  ---------------------
-
     tt = timer; rt = rtimer;
     N = T2[1];
@@ -545 +552 @@
 
 //----------------  Test if we already have a standard basis of I --------------
-
     tt = timer; rt = rtimer;
-    pTest = pTestSB(I,J,L,variant);
-    //pTest = primeTestSB(I,J,L,variant);
+    pTest = primeTest(J, prime(random(1000000000,2134567879)));
     if(printlevel >= 10)
     {
@@ -596 +601 @@
     }
 //--------------  We do not already have a standard basis of I, therefore do the main computation for more primes  --------------
-
     T1 = H;
     T2 = N;
@@ -613 +617 @@
       for(i=j; i<=size(L); i++)
       {
-        //Arguments[i-j+1] = list(II,L[i],variant,list(curr_weight,target_weight));
-        Arguments[size(Arguments)+1] = list(II,L[i],variant,list(curr_weight,target_weight));
+        Arguments[size(Arguments)+1] = list(II,L[i],variant,reduction,list(curr_weight,target_weight));
       }
     }
@@ -621 +624 @@
       for(i=j; i<=size(L); i++)
       {
-        //Arguments[i-j+1] = list(II,L[i],variant);
-        Arguments[size(Arguments)+1] = list(II,L[i],variant);
+        Arguments[size(Arguments)+1] = list(II,L[i],variant,reduction);
       }
     }
@@ -632 +634 @@
     for(i=1; i<=size(PP); i++)
     {
-      //P[size(P) + 1] = PP[i];
       T1[size(T1) + 1] = PP[i][1];
       T2[size(T2) + 1] = bigint(PP[i][2]);
@@ -649 +650 @@
   echo = 2;
   ring R=0,(x,y,z),lp;
-  ideal I=-x+y2z-z,xz+1,x2+y2-1;
-  // I is a standard basis in dp
-  ideal J = modWalk(I,1);
+  ideal I= y3+xyz+y2z+xz3, 3+xy+x2y+y2z;
+  int reduction = 0;
+  ideal J = modWalk(I,1,1);
   J;
 }
@@ -772 +773 @@
   return(J);
 }
-//////////////////////////////////////////////////////////////////////////////////
-static proc primeTestSB(def II, ideal J, list L, int variant, list #)
-"USAGE:  primeTestSB(I,J,L,variant,#); I,J ideals, L intvec of primes, variant int
-RETURN:  1 (resp. 0) if for a randomly chosen prime p that is not in L
-         J mod p is (resp. is not) a standard basis of I mod p
-EXAMPLE: example primeTestSB; shows an example
-"
-{
-if(typeof(II) == "ideal")
-  {
-  ideal I = II;
-  int radius = 2;
-  }
-if(typeof(II) == "list")
-  {
-  ideal I = II[1];
-  int radius = II[2];
-  }
-
-int i,j,k,p;
-def R = basering;
-list r = ringlist(R);
-
-while(!j)
-  {
-  j = 1;
-  p = prime(random(1000000000,2134567879));
-  for(i = 1; i <= size(L); i++)
-    {
-    if(p == L[i])
-      {
-      j = 0;
-      break;
-      }
-    }
-  if(j)
-    {
-    for(i = 1; i <= ncols(I); i++)
-      {
-      for(k = 2; k <= size(I[i]); k++)
-        {
-        if((denominator(leadcoef(I[i][k])) mod p) == 0)
-          {
-          j = 0;
-          break;
-          }
-        }
-      if(!j)
-        {
-        break;
-        }
-      }
-    }
-  if(j)
-    {
-    if(!primeTest(I,p))
-      {
-      j = 0;
-      }
-    }
-  }
-r[1] = p;
-def @R = ring(r);
-setring @R;
-ideal I = imap(R,I);
-ideal J = imap(R,J);
-attrib(J,"isSB",1);
-
-int t = timer;
-j = 1;
-if(isIncluded(I,J) == 0)
-  {
-  j = 0;
-  }
-if(printlevel >= 11)
-  {
-  "isIncluded(I,J) takes "+string(timer - t)+" seconds";
-  "j = "+string(j);
-  }
-t = timer;
-if(j)
-  {
-  if(size(#) > 0)
-    {
-    ideal K = modpWalk(I,p,variant,#)[1];
-    }
-  else
-    {
-    ideal K = modpWalk(I,p,variant)[1];
-    }
-  t = timer;
-  if(isIncluded(J,K) == 0)
-    {
-    j = 0;
-    }
-  if(printlevel >= 11)
-    {
-    "isIncluded(K,J) takes "+string(timer - t)+" seconds";
-    "j = "+string(j);
-    }
-  }
-setring R;
-
-return(j);
-}
-example
-{ "EXAMPLE:"; echo = 2;
-   intvec L = 2,3,5;
-   ring r = 0,(x,y,z),lp;
-   ideal I = x+1,x+y+1;
-   ideal J = x+1,y;
-   primeTestSB(I,I,L,1);
-   primeTestSB(I,J,L,1);
-}
-
-////////////////////////////////////////////////////////////////////////////////
-static proc pTestSB(ideal I, ideal J, list L, int variant, list #)
-"USAGE:  pTestSB(I,J,L,variant,#); I,J ideals, L intvec of primes, variant int
-RETURN:  1 (resp. 0) if for a randomly chosen prime p that is not in L
-         J mod p is (resp. is not) a standard basis of I mod p
-EXAMPLE: example pTestSB; shows an example
-"
-{
-   int i,j,k,p;
-   def R = basering;
-   list r = ringlist(R);
-
-   while(!j)
-   {
-      j = 1;
-      p = prime(random(1000000000,2134567879));
-      for(i = 1; i <= size(L); i++)
-      {
-         if(p == L[i]) { j = 0; break; }
-      }
-      if(j)
-      {
-         for(i = 1; i <= ncols(I); i++)
-         {
-            for(k = 2; k <= size(I[i]); k++)
-            {
-               if((denominator(leadcoef(I[i][k])) mod p) == 0) { j = 0; break; }
-            }
-            if(!j){ break; }
-         }
-      }
-      if(j)
-      {
-         if(!primeTest(I,p)) { j = 0; }
-      }
-   }
-   r[1] = p;
-   def @R = ring(r);
-   setring @R;
-   ideal I = imap(R,I);
-   ideal J = imap(R,J);
-   attrib(J,"isSB",1);
-
-   int t = timer;
-   j = 1;
-   if(isIncluded(I,J) == 0) { j = 0; }
-
-   if(printlevel >= 11)
-   {
-      "isIncluded(I,J) takes "+string(timer - t)+" seconds";
-      "j = "+string(j);
-   }
-
-   t = timer;
-   if(j)
-   {
-      if(size(#) > 0)
-      {
-         ideal K = modpStd(I,p,variant,#[1])[1];
-      }
-      else
-      {
-         ideal K = groebner(I);
-      }
-      t = timer;
-      if(isIncluded(J,K) == 0) { j = 0; }
-
-      if(printlevel >= 11)
-      {
-         "isIncluded(J,K) takes "+string(timer - t)+" seconds";
-         "j = "+string(j);
-      }
-   }
-   setring R;
-   return(j);
-}
-example
-{ "EXAMPLE:"; echo = 2;
-   intvec L = 2,3,5;
-   ring r = 0,(x,y,z),dp;
-   ideal I = x+1,x+y+1;
-   ideal J = x+1,y;
-   pTestSB(I,I,L,2);
-   pTestSB(I,J,L,2);
-}
+
 ////////////////////////////////////////////////////////////////////////////////
 static proc mixedTest()

Singular/LIB/noether.lib

@@ -67 +67 @@
    for ( ii = 1; ii <= size(I); ii++ )
    {
-     Y=findvars(I[ii],1)[1];
+     Y=variables(I[ii]);
      l=rvar(Y[1][1]);
      if (size(Y[1])==1)
@@ -89 +89 @@
    else
    {
-     P2=findvars(ideal(P1[1..t]),1)[3];
+     P2=findvars(ideal(P1[1..t]))[3];
      for ( jj = 1; jj <= size(P2[1]); jj++ )
      {

Singular/LIB/normal.lib

@@ -1148 +1148 @@
           "// the ideal was the zero-ideal";
       }
-     // execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
-     //                 +ordstr(basering)+");");
          def newR7 = ring(gnirlist);
          setring newR7;
@@ -1597 +1595 @@
          // Now RR[2] must be 1, hence the test for normality was positive
          MB=SM[2];
-         //execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
-         //             +ordstr(basering)+");");
          def newR7 = ring(gnirlist);
          setring newR7;
@@ -1901 +1897 @@
       //ideal new2=SL[1],SM[2];
 
-      //execute("ring newR1="+charstr(basering)+",("+varstr(basering)+"),("