Changeset 565e866 in git for Singular/LIB/bfct.lib

Timestamp:

Dec 23, 2008, 10:39:31 PM (15 years ago)

Author:

Viktor Levandovskyy <levandov@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

199ae72040d35b09813089e33464df761b5aad78

Parents:

6045c31371ebacbf1acfc2f082225f05016f06b2

Message:

*levandov: numerous changes, docu fixes, exchange of procs between libs etc


git-svn-id: file:///usr/local/Singular/svn/trunk@11268 2c84dea3-7e68-4137-9b89-c4e89433aadc

File:

• Singular/LIB/bfct.lib (modified) (47 diffs)

Singular/LIB/bfct.lib

@@ -1 +1 @@
 //////////////////////////////////////////////////////////////////////////////
-version="$Id: bfct.lib,v 1.7 2008-12-09 16:50:21 levandov Exp $";
+version="$Id: bfct.lib,v 1.8 2008-12-23 21:39:31 levandov Exp $";
 category="Noncommutative";
 info="
-LIBRARY: bfct.lib     M. Noro's Algorithm for Bernstein-Sato polynomial
+LIBRARY: bfct.lib     Algorithms for b-functions and Bernstein-Sato polynomial
 AUTHORS: Daniel Andres, daniel.andres@math.rwth-aachen.de
 @* Viktor Levandovskyy,      levandov@math.rwth-aachen.de
 
 THEORY: Given a polynomial ring R = K[x_1,...,x_n] and a polynomial F in R,
-@*      one is interested in the global Bernstein-Sato polynomial b(s) in K[s],
-@*      defined to be the monic polynomial, satisfying a functional identity
-@*             L * f^{s+1} = b(s) f^s,   for some operator L in D[s].
-@* Here, D stands for an n-th Weyl algebra K<x_1,...,x_n,d_1,...,d_n | d_j x_j = x_j d_j +1>
-@* One is interested in the following data:
-@*   global Bernstein-Sato polynomial in K[s] and
-@*   the list of all roots of b(s), which are known to be rational, with their multiplicities.
+@*      one is interested in the global b-Function (also known as Bernstein-Sato
+@*      polynomial) b(s) in K[s], defined to be the monic polynomial, satisfying a functional
+@*      identity L * F^{s+1} = b(s) F^s,   for some operator L in D[s]. Here, D stands for an
+@*      n-th Weyl algebra K<x_1,...,x_n,d_1,...,d_n | d_j x_j = x_j d_j +1>.
+@*      One is interested in the following data:
+@*      - Bernstein-Sato polynomial b(s) in K[s],
+@*      - the list of its roots, which are known to be rational
+@*      - the multiplicities of the roots.
+@*      References: Saito, Strurmfels, Takayama: Groebner Deformations of Hypergeometric
+@*      Differential Equations (2000), Noro: An Efficient Modular Algorithm for Computing
+@*      the Global b-function, (2002).
+
 
 MAIN PROCEDURES:
 
-bfct(f[,s,t,v]);            compute the global Bernstein-Sato polynomial of a given poly
-bfctsyz(f[,r,s,t,u,v]);     compute the global Bernstein-Sato polynomial of a given poly
-bfctann(f[,s]);             compute the global Bernstein-Sato polynomial of a given poly
-bfctonestep(f[,s,t]);       compute the global Bernstein-Sato polynomial of a given poly
-bfctideal(I,w[,s,t]);       compute the global b-function of a given ideal w.r.t. a given weight
-pintersect(f,I);            compute the intersection of the ideal generated by a given poly with a given ideal
-pintersectsyz(f,I[,p,s,t]); compute the intersection of the ideal generated by a given poly with a given ideal
-linreduce(f,I[,s]);         reduce a poly by linear reductions only
-ncsolve(I[,s]);             find and compute a linear dependency of the elements of an ideal
+bfct(f[,s,t,v]);            computes the Bernstein-Sato polynomial of poly f
+bfctSyz(f[,r,s,t,u,v]);     computes the Bernstein-Sato polynomial of poly f
+bfctAnn(f[,s]);             computes the Bernstein-Sato polynomial of poly f
+bfctOneGB(f[,s,t]);         computes the Bernstein-Sato polynomial of poly f
+bfctIdeal(I,w[,s,t]);       computes the global b-function of ideal I w.r.t. the weight w
+pIntersect(f,I);            intersection of the ideal I with the subalgebra K[f] for a poly f
+pIntersectSyz(f,I[,p,s,t]);  intersection of the ideal I with the subalgebra K[f] for a poly f
+linReduce(f,I[,s]);         reduces a poly by linear reductions only
+linSyzSolve(I[,s]);         computes a linear dependency of the elements of ideal I
 
 AUXILIARY PROCEDURES:
 
-ispositive(v);   check whether all entries of an intvec are positive 
-isin(l,i);       check whether an element is a member of a list 
-scalarprod(v,w); compute the standard scalar product of two intvecs
-vec2poly(v[,i]); convert a coefficient vector to a poly
+allPositive(v);  checks whether all entries of an intvec are positive 
+scalarProd(v,w); computes the standard scalar product of two intvecs
+vec2poly(v[,i]); constructs an univariate poly with given coefficients
 
 SEE ALSO: dmod_lib, dmodapp_lib, gmssing_lib
@@ -40 +44 @@
 
 LIB "qhmoduli.lib"; // for Max
-LIB "dmod.lib"; // for SannfsBFCT etc
-LIB "dmodapp.lib";  // for initialideal etc
-
-
+LIB "dmod.lib";     // for SannfsBFCT etc
+LIB "dmodapp.lib";  // for initialIdeal etc
+
+
+///////////////////////////////////////////////////////////////////////////////
+// testing for consistency of the library:
 proc testbfctlib ()
 {
   // tests all procs for consistency
   "AUXILIARY PROCEDURES:";
-  example ispositive;
-  example isin;
-  example scalarprod;
+  example allPositive;
+  example scalarProd;
   example vec2poly;
   "MAIN PROCEDURES:";
   example bfct;
-  example bfctsyz;
-  example bfctann;
-  example bfctonestep;
-  example bfctideal;
-  example pintersect;
-  example pintersectsyz;
-  example linreduce;
-  example ncsolve;
-}
-
+  example bfctSyz;
+  example bfctAnn;
+  example bfctOneGB;
+  example bfctIdeal;
+  example pIntersect;
+  example pIntersectSyz;
+  example linReduce;
+  example linReduceIdeal;
+  example linSyzSolve;
+}
 
 //--------------- auxiliary procedures ---------------------------------------------------------
@@ -71 +76 @@
 RETURN:  a ring, the associated graded ring of the basering w.r.t. u and v
 PURPOSE: compute the associated graded ring of the basering w.r.t. u and v
+ASSUME: basering is a Weyl algebra
 EXAMPLE: example gradedWeyl; shows examples
 NOTE:    u[i] is the weight of x(i), v[i] the weight of D(i).
@@ -122 +128 @@
 {
   "EXAMPLE:"; echo = 2;
-  LIB "bfct.lib";
   ring @D = 0,(x,y,z,Dx,Dy,Dz),dp;
   def D = Weyl();
@@ -132 +137 @@
 
 
-proc ispositive (intvec v)
-"USAGE:  ispositive(v);  v an intvec
-RETURN:  1 if all components of v are positive, or 0 otherwise
+proc allPositive (intvec v)
+"USAGE:  allPositive(v);  v an intvec
+RETURN:  int, 1 if all components of v are positive, or 0 otherwise
 PURPOSE: check whether all components of an intvec are positive
-EXAMPLE: example ispositive; shows an example
+EXAMPLE: example allPositive; shows an example
 "
 {
@@ -154 +159 @@
   "EXAMPLE:"; echo = 2;
   intvec v = 1,2,3;
-  ispositive(v);
+  allPositive(v);
   intvec w = 1,-2,3;
-  ispositive(w);
-}
-
-proc isin (list l, i)
-"USAGE:  isin(l,i);  l a list, i an argument of any type
-RETURN:  an int, the position of the first appearance of i in l, or 0 if i is not a member of l
+  allPositive(w);
+}
+
+static proc findFirst (list l, i)
+"USAGE:  findFirst(l,i);  l a list, i an argument of any type
+RETURN:  int, the position of the first appearance of i in l, or 0 if i is not a member of l
 PURPOSE: check whether the second argument is a member of a list
-EXAMPLE: example isin; shows an example
+EXAMPLE: example findFirst; shows an example
 "
 {
@@ -177 +182 @@
   return(0);
 }
+//   isin(list(1, 2, list(1)),2);       // seems to be a bit buggy,
+//   isin(list(1, 2, list(1)),3);       // but works for the purposes
+//   isin(list(1, 2, list(1)),list(1)); // of this library
+//   isin(l,list(2));
 example
 {
@@ -182 +191 @@
   ring r = 0,x,dp;
   list l = 1,2,3;
-  isin(l,4);
-  isin(l,2);
-}
-
-proc scalarprod (intvec v, intvec w)
-"USAGE:  scalarprod(v,w);  v,w intvecs
-RETURN:  an int, the standard scalar product of v and w
-PURPOSE: compute the scalar product of two intvecs
-NOTE:    the arguments must have the same size
-EXAMPLE: example scalarprod; shows examples
+  findFirst(l,4);
+  findFirst(l,2);
+}
+
+proc scalarProd (intvec v, intvec w)
+"USAGE:  scalarProd(v,w);  v,w intvecs
+RETURN:  int, the standard scalar product of v and w
+PURPOSE: computes the scalar product of two intvecs
+ASSUME:  the arguments are of the same size
+EXAMPLE: example scalarProd; shows examples
 "
 {
@@ -213 +222 @@
   intvec v = 1,2,3;
   intvec w = 4,5,6;
-  scalarprod(v,w);
-}
-
+  scalarProd(v,w);
+}
 
 //-------------- main procedures -------------------------------------------------------
 
-
-proc linreduce(poly f, ideal I, list #)
-"USAGE:  linreduce(f, I [,s,t]);  f a poly, I an ideal, s,t optional ints
-RETURN:  a poly obtained by linear reductions of the given poly with the given ideal
-PURPOSE: reduce a poly only by linear reductions
-NOTE:    If s<>0, a list consisting of the reduced poly and the coefficient vector of the used 
-@*       reductions is returned.
-@*       If t<>0, only leading monomials are reduced, otherwise, and by default, all monomials
-@*       are reduced, if possible.
-EXAMPLE: example linreduce; shows examples
+proc linReduceIdeal(ideal I, list #)
+"USAGE:  linReduceIdeal(I [,s,t]); I an ideal, s,t optional ints
+RETURN:  ideal/list, linear reductum (over field) of f by elements from I
+PURPOSE: reduce a poly only by linear reductions (no monomial multiplications)
+NOTE:    If s<>0, a list consisting of the reduced ideal and the coefficient 
+@*       vectors of the used reductions given as module is returned.
+@*       Otherwise (and by default), only the reduced ideal is returned.
+@*       If t=0 (and by default) all monomials are reduced (if possible),
+@*       otherwise, only leading monomials are reduced.
+DISPLAY: If @code{printlevel}>=1, all debug messages will be printed.
+EXAMPLE: example linReduceIdeal; shows examples
 "
 {
+  // #[1] = remembercoeffs
+  // #[2] = redtail
   int ppl = printlevel - voice + 2;
   int remembercoeffs = 0; // default
-  int redlm          = 0; // default
+  int redtail        = 0; // default
   if (size(#)>0)
   {
@@ -244 +255 @@
       if (typeof(#[2])=="int" || typeof(#[2])=="number")
       {
-        redlm = #[2];
+        redtail = #[2];
+      }
+    }
+  }
+  int sI = ncols(I);
+  int sZeros = sI - size(I);
+  dbprint(ppl,"ideal contains "+string(sZeros)+" zero(s)");
+  int i,j;
+  ideal J,lmJ,ordJ;
+  list lJ = sort(I);
+  module M; // for the coefficients
+  if (sZeros > 0) // I contains zeros
+  {
+    if (remembercoeffs <> 0)
+    {
+      j = 0;
+      for (i=1; i<=sI; i++)
+      {
+        if (I[i] == 0)
+        {
+          j++;
+          J[j] = 0;
+          ordJ[j] = -1;
+          M[j] = gen(i);
+        }
+        else
+        {
+          M[i+sZeros-j] = gen(lJ[2][i-j]+j);
+        }
+      }
+    }
+    else
+    {
+      for (i=1; i<=sZeros; i++)
+      {
+        J[i] = 0;
+        ordJ[i] = -1;
+      }
+    }
+    I = J,lJ[1];
+  }
+  else
+  {
+    I = lJ[1];
+    if (remembercoeffs <> 0)
+    {
+      for (i=1; i<=size(lJ[1]); i++) { M[i+sZeros] = gen(lJ[2][i]); }
+    }
+  }
+  dbprint(ppl,"initially sorted ideal:", I);
+  if (remembercoeffs <> 0) { dbprint(ppl," used permutations:", M); }
+  poly lm,c,redpoly,maxlmJ;
+  J[sZeros+1] = I[sZeros+1];
+  lm = I[sZeros+1];
+  maxlmJ = leadmonom(lm);
+  lmJ[sZeros+1] = maxlmJ;
+  int ordlm,reduction,maxordJ;
+  maxordJ = ord(lm);
+  ordJ[sZeros+1] = maxordJ;
+  vector v;
+  for (i=sZeros+2; i<=sI; i++)
+  {
+    redpoly = I[i];
+    lm = leadmonom(redpoly);
+    ordlm = ord(lm);
+    reduction = 1;
+    if (remembercoeffs <> 0) { v = M[i]; }
+    while (reduction == 1) // while there was a reduction
+    {
+      reduction = 0;
+      for (j=sZeros+1; j<i; j++)
+      {
+        if (lm == 0) { break; }
+        if (ordlm > maxordJ) { break; }
+        if (ordlm == ordJ[j])
+        {           
+          if (lm > maxlmJ) { break; }
+          if (lm == lmJ[j])
+          {
+            dbprint(ppl,"reducing " + string(redpoly));
+            dbprint(ppl," with " + string(J[j]));
+            c = leadcoef(redpoly)/leadcoef(J[j]);
+            redpoly = redpoly - c*J[j];
+            dbprint(ppl," to " + string(redpoly));
+            lm = leadmonom(redpoly);
+            ordlm = ord(lm);
+            if (remembercoeffs <> 0) { M[i] = M[i] - c * M[j]; }
+            reduction = 1;
+          }
+        }
+      }
+    }
+    for (j=sZeros+1; j<i; j++)
+    {
+      if (redpoly < J[j]) { break; }
+    }
+    J = insertGenerator(J,redpoly,j);
+    lmJ = insertGenerator(lmJ,lm,j);
+    ordJ = insertGenerator(ordJ,poly(ordlm),j);
+    if (remembercoeffs <> 0) 
+    {
+      v = M[i];
+      M = deleteGenerator(M,i);
+      M = insertGenerator(M,v,j);
+    }
+    maxordJ = ord(J[i]);
+    maxlmJ = lmJ[i];
+  }
+  if (redtail <> 0)
+  {
+    dbprint(ppl,"finished reducing leading monomials:",J);
+    if (remembercoeffs <> 0) { dbprint(ppl,"used reductions:",M); }
+    int k;
+    for (i=sZeros+1; i<=sI; i++)
+    {
+      lm = lmJ[i];
+      for (j=i+1; j<=sI; j++)
+      {
+        for (k=2; k<=size(J[j]); k++) // run over all terms in J[j]
+        {
+          if (ordJ[i] == ord(J[j][k]))
+          {    
+            if (lm == normalize(J[j][k]))
+            {
+              c = leadcoef(J[j][k])/leadcoef(J[i]);
+              dbprint(ppl,"reducing " + string(J[j]));
+              dbprint(ppl," with " + string(J[i]));
+              J[j] = J[j] - c*J[i];
+              dbprint(ppl," to " + string(J[j]));
+              if (remembercoeffs <> 0) { M[j] = M[j] - c * M[i]; }
+            }
+          }
+        }
+      }
+    }
+  }
+  if (remembercoeffs <> 0) { return(list(J,M)); }
+  else { return(J); }
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  ring r = 0,(x,y),dp;
+  ideal I = 3,x+9,y4,y4+7x+2;
+  linReduceIdeal(I);
+  linReduceIdeal(I,0,1);
+  list l = linReduceIdeal(I,1,1); l;
+  module M = I;
+  l[1] - ideal(M*l[2]);
+}
+
+proc linReduce(poly f, ideal I, list #)
+"USAGE:  linReduce(f, I [,s,t]); f a poly, I an ideal, s,t optional ints
+RETURN:  poly/list, linear reductum (over field) of f by elements from I
+PURPOSE: reduce a poly only by linear reductions (no monomial multiplications)
+NOTE:    If s<>0, a list consisting of the reduced poly and the coefficient 
+@*       vector of the used reductions is returned, otherwise (and by default)
+@*       only reduced poly is returned.
+@*       If t=0 (and by default) all monomials are reduced (if possible),
+@*       otherwise, only leading monomials are reduced.
+DISPLAY: If @code{printlevel}>=1, all debug messages will be printed.
+EXAMPLE: example linReduce; shows examples
+"
+{
+  int ppl = printlevel - voice + 2;
+  int remembercoeffs = 0; // default
+  int redtail        = 0; // default
+  int prepareideal   = 1; // default
+  if (size(#)>0)
+  {
+    if (typeof(#[1])=="int" || typeof(#[1])=="number")
+    {
+      remembercoeffs = #[1];
+    }
+    if (size(#)>1)
+    {
+      if (typeof(#[2])=="int" || typeof(#[2])=="number")
+      {
+        redtail = #[2];
+      }
+      if (size(#)>2)
+      {
+        if (typeof(#[3])=="int" || typeof(#[3])=="number")
+        {
+          prepareideal = #[3];
+        }
       }
     }
@@ -250 +446 @@
   int i,j,k;
   int sI = ncols(I);
-  ideal lmI,lcI;
+  // pre-reduce I:
+  module M;
+  if (prepareideal <> 0)
+  {
+    if (remembercoeffs <> 0)
+    {
+      list sortedI = linReduceIdeal(I,1,redtail);
+      I = sortedI[1];
+      M = sortedI[2];
+      dbprint(ppl,"prepared ideal:",I);
+      dbprint(ppl," with operations:",M);
+    }
+    else
+    {
+      I = linReduceIdeal(I,0,redtail);
+    }
+  }
+  else
+  {
+    if (remembercoeffs <> 0)
+    {
+      for (i=1; i<=sI; i++)
+      {
+        M[i] = gen(i);
+      }
+    }
+  }
+  ideal lmI,lcI,ordI;
   for (i=1; i<=sI; i++)
   {
     lmI[i] = leadmonom(I[i]);
     lcI[i] = leadcoef(I[i]);
+    ordI[i] = ord(lmI[i]);
   }
   vector v;
   poly c;
-  int reduction = 1;
-  if (redlm == 0)
-  {
-    ideal monomf;
-    for (k=1; k<=size(f); k++)
-    {
-      monomf[k] = normalize(f[k]);
-    }
-    while (reduction == 1) // while there was a reduction
-    {
-      reduction = 0;
-      for (i=sI; i>=1; i--)
-      {
-        dbprint(ppl,"testing ideal entry:",i);
-        for (j=1; j<=size(f); j++)
+  // === reduce leading monomials ===
+  poly lm = leadmonom(f);
+  int ordf = ord(lm);
+  for (i=sI; i>=1; i--)  // I is sorted: smallest lm's on top
+  {
+    if (lm == 0)
+    {
+      break;
+    }
+    else
+    {
+      if (ordf == ordI[i])
+      {
+        if (lm == lmI[i])
         {
-          if (monomf[j] == lmI[i])
+          c = leadcoef(f)/lcI[i];
+          f = f - c*I[i];
+          lm = leadmonom(f);
+          ordf = ord(lm);
+          if (remembercoeffs <> 0)
+          {
+            v = v - c * M[i];
+          }
+        }
+      }
+    }
+  }
+  // === reduce tails as well ===
+  if (redtail <> 0)
+  {
+    dbprint(ppl,"poly after reduction of leading monomials:",f);
+    for (i=1; i<=sI; i++)
+    {
+      dbprint(ppl,"testing ideal entry:",i);
+      for (j=1; j<=size(f); j++)
+      {
+        if (ord(f[j]) == ordI[i])
+        {
+          if (normalize(f[j]) == lmI[i])
           {
             c = leadcoef(f[j])/lcI[i];
             f = f - c*I[i];
             dbprint(ppl,"reducing poly to ",f);
-            monomf = 0;
-            for (k=1; k<=size(f); k++)
-            {
-              monomf[k] = normalize(f[k]);
-            }
-            reduction = 1;
             if (remembercoeffs <> 0)
             {
-              v = v - c * gen(i);
+              v = v - c * M[i];
             }
-            break;
-          }
-        }
-        if (reduction == 1)
-        {
-          break;
-        }
-      }
-    }
-  }
-  else // reduce only leading monomials
-  {
-    poly lm = leadmonom(f);
-    while (reduction == 1) // while there was a reduction
-    {
-      reduction = 0;
-      for (i=sI;i>=1;i--)
-      {
-        if (lm <> 0 && lm == lmI[i])
-        {
-          c = leadcoef(f)/lcI[i];
-          f = f - c*I[i];
-          lm = leadmonom(f);
-          reduction = 1;
-          if (remembercoeffs <> 0)
-          {
-            v = v - c * gen(i);
           }
         }
@@ -338 +551 @@
   poly f = 5xy+7y+3;
   poly g = 7x+5y+3;
-  linreduce(g,I);
-  linreduce(g,I,0,1);
-  linreduce(f,I,1);
-}
-
-proc ncsolve (ideal I, list #)
-"USAGE:  ncsolve(I[,s]);  I an ideal, s an optional int
-RETURN:  coefficient vector of a linear combination of 0 in the elements of I
-PURPOSE: compute a linear dependency between the elements of an ideal if such one exists
+  linReduce(g,I);
+  linReduce(g,I,0,1);
+  linReduce(f,I,1);
+  f = x3 + y2 + x2 + y + x;
+  I = x3 - y3, y3 - x2, x3 - y2, x2 - y, y2-x;
+  list l = linReduce(f, I, 1); 
+  l;
+  module M = I;
+  f - (l[1] - (M*l[2])[1,1]);
+}
+
+proc linSyzSolve (ideal I, list #)
+"USAGE:  linSyzSolve(I[,s]);  I an ideal, s an optional int
+RETURN:  vector, coefficient vector of a linear combination of 0 in the elements of I
+PURPOSE: compute a linear dependency between the elements of an ideal
+@*       if such one exists
 NOTE:    If s<>0, @code{std} is used for Groebner basis computations,
 @*       otherwise, @code{slimgb} is used.
 @*       By default, @code{slimgb} is used in char 0 and @code{std} in char >0.
-@*       If printlevel=1, progress debug messages will be printed,
+DISPLAY: If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
-EXAMPLE: example ncsolve; shows examples
+EXAMPLE: example linSyzSolve; shows examples
 "
 {
@@ -374 +594 @@
   else
   {
-    // 1. introduce undefined coeffs
+    // ------- 1. introduce undefined coeffs ------------------
     def save = basering;
     int p = char(save);
@@ -386 +606 @@
     ring @A = p,(@a(1..sI)),lp;
     ring @aA = (p,@a(1..sI)), (@z),dp;
+    // TODO: BUG! WHAT IF THE ORIGINAL RING ALREADY HAS SUCH VARIABLES/PARAMETERS!!!?
+    // TODO: YOU CAN OVERCOME THIS BY MEANS SIMILAR TO "chooseSafeVarName" FROM NEW "matrix.lib"
     def @B = save + @aA;
     setring @B;
     ideal I = imap(save,I);
-    // 2. form the linear system for the undef coeffs
+    // ------- 2. form the linear system for the undef coeffs ---
     int i;   poly W;  ideal QQ;
     for (i=1; i<=sI; i++)
@@ -403 +625 @@
     setring @A;
     ideal QQ = imap(@B,QQ);
-    // 3. this QQ is a polynomial ideal, so "solve" the system
-    dbprint(ppl, "ncsolve: starting Groebner basis computation with engine:", whichengine);
+    // ------- 3. this QQ is a polynomial ideal, so "solve" the system -----
+    dbprint(ppl, "linSyzSolve: starting Groebner basis computation with engine:", whichengine);
     QQ = engine(QQ,whichengine);
     dbprint(ppl, "QQ after engine:", QQ);
     if (dim(QQ) == -1)
     {
-      dbprint(ppl+1, "no solutions by ncsolve");
+      dbprint(ppl+1, "no solutions by linSyzSolve");
       // output zeroes
       setring save;
@@ -415 +637 @@
       return(v);
     }
-    // 4. in order to get the numeric values
+    // ------- 4. in order to get the numeric values -------
     matrix AA = matrix(maxideal(1));
     module MQQ = std(module(QQ));
@@ -439 +661 @@
   ring r = 0,x,dp;
   ideal I = x,2x;
-  ncsolve(I);
+  linSyzSolve(I);
   ideal J = x,x2;
-  ncsolve(J);
-}
-
-proc pintersect (poly s, ideal I)
-"USAGE:  pintersect(f, I);  f a poly, I an ideal
-RETURN:  coefficient vector of the monic generator of the intersection of the ideal generated by f with I
-PURPOSE: compute the intersection of an ideal with a principal ideal defined by f
+  linSyzSolve(J);
+}
+
+proc pIntersect (poly s, ideal I)
+"USAGE:  pIntersect(f, I);  f a poly, I an ideal
+RETURN:  vector, coefficient vector of the monic polynomial 
+PURPOSE: compute the intersection of ideal I with the subalgebra K[f]
+ASSUME:  I is given as Groebner basis.
 NOTE:    If the intersection is zero, this proc might not terminate.
-@*       I should be given as standard basis.
-@*       If printlevel=1, progress debug messages will be printed,
+DISPLAY: If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
-EXAMPLE: example pintersect; shows examples
+EXAMPLE: example pIntersect; shows examples
 "
 {
   // assume I is given in Groebner basis
-  attrib(I,"isSB",1); // set attribute for suppressing NF messages
+  if (attrib(I,"isSB") <> 1)
+  {
+    print("WARNING: The input has no SB attribute!");
+    print(" Treating it as if it were a Groebner basis and proceeding...");
+    attrib(I,"isSB",1); // set attribute for suppressing NF messages
+  }
   int ppl = printlevel-voice+2;
   // ---case 1: I = basering---
   if (size(I) == 1)
   {
-    if (simplify(I,1) == ideal(1))
+    if (simplify(I,2)[1] == 1)
     {
       return(gen(2)); // = s
@@ -512 +739 @@
   poly toNF,tobracket,newNF,rednewNF,oldNF,secNF;
   i = 1;
-  dbprint(ppl+1,"pintersect starts...");
+  dbprint(ppl+1,"pIntersect starts...");
   while (1)
   {
@@ -535 +762 @@
       secNF = newNF;
     }
-    ll = linreduce(newNF,redNI,1);
+    ll = linReduce(newNF,redNI,1);
     rednewNF = ll[1];
     l[i+1] = ll[2];
@@ -578 +805 @@
   v = imap(@R,v);
   kill @R;
-  dbprint(ppl+1,"pintersect finished");
+  dbprint(ppl+1,"pIntersect finished");
   return(v);
 }
@@ -586 +813 @@
   ring r = 0,(x,y),dp;
   poly f = x^2+y^3+x*y^2;
-  def D = initialmalgrange(f);
+  def D = initialMalgrange(f);
   setring D;
   inF;
-  pintersect(t*Dt,inF);
-}
-
-proc pintersectsyz (poly s, ideal II, list #)
-"USAGE:  pintersectsyz(f, I [,p,s,t]);  f a poly, I an ideal, p, t optial ints, p a prime number
-RETURN:  coefficient vector of the monic generator of the intersection of the ideal generated by f with I
-PURPOSE: compute the intersection of an ideal with a principal ideal defined by f
-NOTE:    If the intersection is zero, this proc might not terminate.
-@*       I should be given as standard basis.
-@*       If p>0 is given, this proc computes the generator of the intersection in char p first and
-@*       then only searches for a generator of the obtained degree in the basering.
-@*       Otherwise, it searched for all degrees.
-@*       This is done by computing syzygies.
+  pIntersect(t*Dt,inF);
+}
+
+proc pIntersectSyz (poly s, ideal II, list #)
+"USAGE:  pIntersectSyz(f, I [,p,s,t]);  f a poly, I an ideal, p, t optial ints, p a prime number
+RETURN:  vector, coefficient vector of the monic polynomial
+PURPOSE: compute the intersection of an ideal I with the subalgebra K[f]
+ASSUME:  I is given as Groebner basis.
+NOTE:    If the intersection is zero, this procedure might not terminate.
+@*       If p>0 is given, this proc computes the generator of the intersection in char p first
+@*       and then only searches for a generator of the obtained degree in the basering.
+@*       Otherwise, it searched for all degrees by computing syzygies.
 @*       If s<>0, @code{std} is used for Groebner basis computations in char 0,
 @*       otherwise, and by default, @code{slimgb} is used.
 @*       If t<>0 and by default, @code{std} is used for Groebner basis computations in char >0,
 @*       otherwise, @code{slimgb} is used.
-@*       If printlevel=1, progress debug messages will be printed,
+DISPLAY: If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
-EXAMPLE: example pintersectsyz; shows examples
+EXAMPLE: example pIntersectSyz; shows examples
 "
 {
   // assume I is given in Groebner basis
   ideal I = II;
-  attrib(I,"isSB",1); // set attribute for suppressing NF messages
+  if (attrib(I,"isSB") <> 1)
+  {
+    print("WARNING: The input has no SB attribute!");
+    print(" Treating it as if it were a Groebner basis and proceeding...");
+    attrib(I,"isSB",1); // set attribute for suppressing NF messages
+  }
   int ppl = printlevel-voice+2;
   int whichengine  = 0; // default
@@ -668 +899 @@
   }
   i = 1;
-  dbprint(ppl+1,"pintersectsyz starts...");
+  dbprint(ppl+1,"pIntersectSyz starts...");
   dbprint(ppl+1,"with ideal I=", I);
   while (1)
@@ -688 +919 @@
     }
     // look for a solution
-    dbprint(ppl,"ncsolve starts with: "+string(matrix(NI)));
+    dbprint(ppl,"linSyzSolve starts with: "+string(matrix(NI)));
     if (solveincharp<>0) // modular method
     {
       setring Rp;
       NI[i+1] = phi(newNF);
-      v = ncsolve(NI,modengine);
+      v = linSyzSolve(NI,modengine);
       if (v!=0) // there is a modular solution
       {
         dbprint(ppl,"got solution in char ",solveincharp," of degree " ,i);
         setring save;
-        v = ncsolve(NI,whichengine);
+        v = linSyzSolve(NI,whichengine);
         if (v==0)
         {
@@ -712 +943 @@
     else // non-modular method
     {
-      v = ncsolve(NI,whichengine);
+      v = linSyzSolve(NI,whichengine);
     }
     matrix MM[1][nrows(v)] = matrix(v);
-    dbprint(ppl,"ncsolve ready  with: "+string(MM));
+    dbprint(ppl,"linSyzSolve ready  with: "+string(MM));
     kill MM;
-    //  "ncsolve ready with"; print(v);
+    //  "linSyzSolve ready with"; print(v);
     if (v!=0)
     {
@@ -729 +960 @@
       if (p!=0)
       {
-        dbprint(ppl,"ncsolve: bad solution!");
+        dbprint(ppl,"linSyzSolve: bad solution!");
       }
       else
       {
-        dbprint(ppl,"ncsolve: got solution!");
+        dbprint(ppl,"linSyzSolve: got solution!");
         // "got solution!";
         break;
@@ -741 +972 @@
     i++;
   }
-  dbprint(ppl+1,"pintersectsyz finished");
+  dbprint(ppl+1,"pIntersectSyz finished");
   return(v);
 }
@@ -749 +980 @@
   ring r = 0,(x,y),dp;
   poly f = x^2+y^3+x*y^2;
-  def D = initialmalgrange(f);
+  def D = initialMalgrange(f);
   setring D;
   inF;
   poly s = t*Dt;
-  pintersectsyz(s,inF);
+  pIntersectSyz(s,inF);
   int p = prime(20000);
-  pintersectsyz(s,inF,p,0,0);
+  pIntersectSyz(s,inF,p,0,0);
 }
 
 proc vec2poly (list #)
 "USAGE:  vec2poly(v [,i]);  v a vector or an intvec, i an optional int
-RETURN:  a poly with coefficient vector v
-PURPOSE: convert a coefficient vector to a poly
-NOTE:    If i>0 is given, the returned poly is an element of K[var(i)], 
-@*       otherwise, and by default, @code{i=1} is used.
-@*       The first entry of v is the coefficient of 1.
+RETURN:  poly, an univariate poly in i-th variable with coefficients given by v
+PURPOSE: constructs an univariate poly in K[var(i)] with given coefficients,
+@*       such that the coefficient at var(i)^{j-1} is v[j].
+NOTE:    The optional argument i must be positive, by default i is 1.
 EXAMPLE: example vec2poly; shows examples
 "
@@ -858 +1088 @@
 }
 
-static proc bfctengine (poly f, int inorann, int whichengine, int methodord, int methodpintersect, int pintersectchar, int modengine, intvec u0)
+static proc bfctengine (poly f, int inorann, int whichengine, int methodord, int methodpIntersect, int pIntersectchar, int modengine, intvec u0)
 {
   int ppl = printlevel - voice +2;
@@ -866 +1096 @@
   if (inorann == 0) // bfct using initial ideal
   {
-    def D = initialmalgrange(f,whichengine,methodord,1,u0);
+    def D = initialMalgrange(f,whichengine,methodord,1,u0);
     setring D;
     ideal J = inF;
@@ -881 +1111 @@
   vector b;
   // try it modular
-  if (methodpintersect <> 0) // pintersectsyz
-  {
-    if (pintersectchar == 0) // pintersectsyz::modular
+  if (methodpIntersect <> 0) // pIntersectSyz
+  {
+    if (pIntersectchar == 0) // pIntersectSyz::modular
     {
       int lb = 30000;
@@ -893 +1123 @@
         dbprint(ppl,"number of run in the loop: "+string(i));
         int q = prime(random(lb,ub));
-        if (isin(usedprimes,q)==0) // if q was not already used
+        if (findFirst(usedprimes,q)==0) // if q was not already used
         {
           usedprimes = usedprimes,q;
           dbprint(ppl,"used prime is: "+string(q));
-          b = pintersectsyz(s,J,q,whichengine,modengine);
+          b = pIntersectSyz(s,J,q,whichengine,modengine);
         }
         i++;
       }
     }
-    else // pintersectsyz::non-modular 
-    {
-      b = pintersectsyz(s,J,0,whichengine);
-    }
-  }
-  else // pintersect: linreduce
-  {
-    b = pintersect(s,J);
+    else // pIntersectSyz::non-modular 
+    {
+      b = pIntersectSyz(s,J,0,whichengine);
+    }
+  }
+  else // pIntersect: linReduce
+  {
+    b = pIntersect(s,J);
   }
   setring save;
@@ -926 +1156 @@
 proc bfct (poly f, list #)
 "USAGE:  bfct(f [,s,t,v]);  f a poly, s,t optional ints, v an optional intvec
-RETURN:  list of roots of the Bernstein-Sato polynomial bs(f) and their multiplicies
-PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface, defined by f, according to the algorithm by Masayuki Noro
-NOTE:    In this proc, the initial Malgrange ideal is computed.
-@*       Further, a system of linear equations is solved by linear reductions.
-@*       If s<>0, @code{std} is used for Groebner basis computations,
+RETURN:  list of ideal and intvec
+PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s)
+@*       for the hypersurface defined by f. 
+ASSUME:  The basering is a commutative polynomial ring in char 0.
+BACKGROUND: In this proc, the initial Malgrange ideal is computed according to the algorithm 
+@*       by Masayuki Noro and then a system of linear equations is solved by linear reductions.
+NOTE:    In the output list, the ideal contains all the roots 
+@*       and the intvec their multiplicities.
+@*       If s<>0, @code{std} is used for GB computations,
 @*       otherwise, and by default, @code{slimgb} is used. 
 @*       If t<>0, a matrix ordering is used for Groebner basis computations,
 @*       otherwise, and by default, a block ordering is used.
-@*       If v is a positive weight vector, v is used for homogenization computations,
+@*       If v is a positive weight vector, v is used for homogenization computations, 
 @*       otherwise and by default, no weights are used.
-@*       If printlevel=1, progress debug messages will be printed,
+DISPLAY: If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
 EXAMPLE: example bfct; shows examples
@@ -965 +1199 @@
       if (size(#)>2)
       {
-        if (typeof(#[3])=="intvec" && size(#[3])==n && ispositive(#[3])==1)
+        if (typeof(#[3])=="intvec" && size(#[3])==n && allPositive(#[3])==1)
         {
           u0 = #[3];
@@ -985 +1219 @@
 }
 
-proc bfctsyz (poly f, list #)
-"USAGE:  bfctsyz(f [,r,s,t,u,v]);  f a poly, r,s,t,u optional ints, v an optional intvec
-RETURN:  list of roots of the Bernstein-Sato polynomial bs(f) and its multiplicies
-PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface, defined by f, according to the algorithm by Masayuki Noro
-NOTE:    In this proc, the initial Malgrange ideal is computed.
-@*       Further, a system of linear equations is solved by computing syzygies.
-@*       If r<>0, @code{std} is used for Groebner basis computations in characteristic 0,
+proc bfctSyz (poly f, list #)
+"USAGE:  bfctSyz(f [,r,s,t,u,v]);  f a poly, r,s,t,u optional ints, v an optional intvec
+RETURN:  list of ideal and intvec
+PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) 
+@*       for the hypersurface defined by f
+ASSUME:  The basering is a commutative polynomial ring in char 0.
+BACKGROUND: In this proc, the initial Malgrange ideal is computed according to the algorithm 
+@*       by Masayuki Noro and then a system of linear equations is solved by computing syzygies.
+NOTE:    In the output list, the ideal contains all the roots 
+@*       and the intvec their multiplicities.
+@*       If r<>0, @code{std} is used for GB computations in characteristic 0, 
 @*       otherwise, and by default, @code{slimgb} is used. 
-@*       If s<>0, a matrix ordering is used for Groebner basis computations,
-@*       otherwise, and by default, a block ordering is used.
-@*       If t<>0, the computation of the intersection is solely performed over charasteristic 0,
-@*       otherwise and by default, a modular method is used.
-@*       If u<>0 and by default, @code{std} is used for Groebner basis computations in characteristic >0,
-@*       otherwise, @code{slimgb} is used. 
-@*       If v is a positive weight vector, v is used for homogenization computations,
-@*       otherwise and by default, no weights are used.
-@*       If printlevel=1, progress debug messages will be printed,
+@*       If s<>0, a matrix ordering is used for GB computations, otherwise, 
+@*       and by default, a block ordering is used.
+@*       If t<>0, the computation of the intersection is solely performed over 
+@*       charasteristic 0, otherwise and by default, a modular method is used.
+@*       If u<>0 and by default, @code{std} is used for GB computations in 
+@*       characteristic >0, otherwise, @code{slimgb} is used. 
+@*       If v is a positive weight vector, v is used for homogenization 
+@*       computations, otherwise and by default, no weights are used.
+DISPLAY: If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
 EXAMPLE: example bfct; shows examples
@@ -1018 +1256 @@
   int whichengine = 0; // default
   int methodord   = 0; // default
-  int pintersectchar  = 0; // default
+  int pIntersectchar  = 0; // default
   int modengine   = 1; // default
   intvec u0       = 0; // default
@@ -1037 +1275 @@
         if (typeof(#[3])=="int" || typeof(#[3])=="number")
         {
-          pintersectchar = int(#[3]);
+          pIntersectchar = int(#[3]);
         }
         if (size(#)>3)
@@ -1047 +1285 @@
           if (size(#)>4)
           {
-            if (typeof(#[5])=="intvec" && size(#[5])==n && ispositive(#[5])==1)
+            if (typeof(#[5])=="intvec" && size(#[5])==n && allPositive(#[5])==1)
             {
               u0 = #[5];
@@ -1056 +1294 @@
     }
   }
-  list b = bfctengine(f,0,whichengine,methodord,1,pintersectchar,modengine,u0);
+  list b = bfctengine(f,0,whichengine,methodord,1,pIntersectchar,modengine,u0);
   return(b);
 }
@@ -1064 +1302 @@
   ring r = 0,(x,y),dp;
   poly f = x^2+y^3+x*y^2;
-  bfctsyz(f);
+  bfctSyz(f);
   intvec v = 3,2;
-  bfctsyz(f,0,1,1,0,v);
-}
-
-proc bfctideal (ideal I, intvec w, list #)
-"USAGE:  bfctideal(I,w[,s,t]);  I an ideal, w an intvec, s,t optional ints
-RETURN:  list of roots and their multiplicies of the global b-function of I w.r.t. the weight vector (-w,w)
-PURPOSE: compute the global b-function of an ideal according to the algorithm by M. Noro
-NOTE:    Assume, I is an ideal in the n-th Weyl algebra where the sequence of the
-@*       variables is x(1),...,x(n),D(1),...,D(n).
-@*       If s<>0, @code{std} is used for Groebner basis computations in characteristic 0,
+  bfctSyz(f,0,1,1,0,v);
+}
+
+proc bfctIdeal (ideal I, intvec w, list #)
+"USAGE:  bfctIdeal(I,w[,s,t]);  I an ideal, w an intvec, s,t optional ints
+RETURN:  list of ideal and intvec
+PURPOSE: computes the roots of the global b-function of I wrt the weight vector (-w,w).
+ASSUME:  The basering is the n-th Weyl algebra in char 0, where the sequence of
+@*       the variables is x(1),...,x(n),D(1),...,D(n).
+BACKGROUND:  In this proc, the initial ideal of I is computed according to the algorithm by
+@*       Masayuki Noro and then a system of linear equations is solved by linear reductions.
+NOTE:    In the output list, the ideal contains all the roots 
+@*       and the intvec their multiplicities.
+@*       If s<>0, @code{std} is used for GB computations in characteristic 0,
 @*       otherwise, and by default, @code{slimgb} is used. 
-@*       If t<>0, a matrix ordering is used for Groebner basis computations,
-@*       otherwise, and by default, a block ordering is used.
-@*       If printlevel=1, progress debug messages will be printed,
+@*       If t<>0, a matrix ordering is used for GB computations, otherwise,
+@*       and by default, a block ordering is used.
+DISPLAY: If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
 EXAMPLE: example bfctideal; shows examples
@@ -1104 +1346 @@
     }
   }
-  ideal J = initialideal(I,-w,w,whichengine,methodord);
+  ideal J = initialIdeal(I,-w,w,whichengine,methodord);
   poly s;
   for (i=1; i<=n; i++)
@@ -1110 +1352 @@
     s = s + w[i]*var(i)*var(n+i);
   }
-  vector b = pintersect(s,J);
+  vector b = pIntersect(s,J);
   list l = listofroots(b,0);
   return(l);
@@ -1120 +1362 @@
   def D = Weyl();
   setring D;
-  ideal I = 3*x^2*Dy+2*y*Dx,2*x*Dx+3*y*Dy+6;
+  ideal I = std(3*x^2*Dy+2*y*Dx,2*x*Dx+3*y*Dy+6); I;
   intvec w1 = 1,1;
   intvec w2 = 1,2;
   intvec w3 = 2,3;
-  bfctideal(I,w1);
-  bfctideal(I,w2,1);
-  bfctideal(I,w3,0,1);
-}
-
-
-proc bfctonestep (poly f,list #)
-"USAGE:  bfctonestep(f [,s,t]);  f a poly, s,t optional ints
-RETURN:  list of roots of the Bernstein-Sato polynomial bs(f) and its multiplicies
-PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface, defined by f, using only one Groebner basis computation
-NOTE:    If s<>0, @code{std} is used for the Groebner basis computation, otherwise,
+  bfctIdeal(I,w1);
+  bfctIdeal(I,w2,1);
+  bfctIdeal(I,w3,0,1);
+}
+
+proc bfctOneGB (poly f,list #)
+"USAGE:  bfctOneGB(f [,s,t]);  f a poly, s,t optional ints
+RETURN:  list of ideal and intvec
+PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) for the 
+@*       hypersurface defined by f, using only one GB computation
+ASSUME:  The basering is a commutative polynomial ring in char 0.
+BACKGROUND: In this proc, the initial Malgrange ideal is computed based on the 
+@*       algorithm by Masayuki Noro and combined with an elimination ordering.
+NOTE:    In the output list, the ideal contains all the roots 
+@*       and the intvec their multiplicities.
+@*       If s<>0, @code{std} is used for the GB computation, otherwise,
 @*       and by default, @code{slimgb} is used. 
-@*       If t<>0, a matrix ordering is used for Groebner basis computations,
+@*       If t<>0, a matrix ordering is used for GB computations,
 @*       otherwise, and by default, a block ordering is used.
-@*       If printlevel=1, progress debug messages will be printed,
+DISPLAY: If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
-EXAMPLE: example bfctonestep; shows examples
+EXAMPLE: example bfctOneGB; shows examples
 "
 {
@@ -1146 +1393 @@
   def save = basering;
   int n = nvars(save);
+  int noofvars = 2*n+4;
   int i;
   int whichengine = 0; // default
@@ -1163 +1411 @@
     }
   }
-  def DDh = initialidealengine("bfctonestep", whichengine, methodord, f);
+  intvec uv;
+  uv[n+3] = 1;
+  ring r = 0,(x(1..n)),dp;
+  poly f = fetch(save,f);
+  uv[1] = -1; uv[noofvars] = 0;
+  //  for the ordering
+  intvec @a; @a = 1:noofvars; @a[2] = 2;
+  intvec @a2 = @a; @a2[2] = 0; @a2[2*n+4] = 0;
+  if (methodord == 0) // default: block ordering
+  {
+    ring Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),a(uv),dp(noofvars-1),lp(1));
+  }
+  else // M() ordering
+  {
+    intmat @Ord[noofvars][noofvars];
+    @Ord[1,1..noofvars] = uv;
+    @Ord[2,1..noofvars] = 1:(noofvars-1);
+    for (i=1; i<=noofvars-2; i++)
+    {
+      @Ord[2+i,noofvars - i] = -1;
+    }
+    dbprint(ppl,"weights for ordering:",transpose(@a));
+    dbprint(ppl,"the ordering matrix:",@Ord);
+    ring Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),M(@Ord));
+  }
+  dbprint(ppl,"the ring Dh:",Dh);
+  // non-commutative relations
+  matrix @relD[noofvars][noofvars];
+  @relD[1,2] = t*h^2;// s*t = t*s+t*h^2
+  @relD[2,n+3] = Dt*h^2;// Dt*s = s*Dt+h^2
+  @relD[1,n+3] = h^2;
+  for (i=1; i<=n; i++)
+  {
+    @relD[i+2,n+3+i] = h^2;
+  }
+  dbprint(ppl,"nc relations:",@relD);
+  def DDh = nc_algebra(1,@relD);
   setring DDh;
-  dbprint(ppl, "the initial ideal:", string(matrix(inF)));
+  dbprint(ppl,"computing in ring",DDh);
+  ideal I;
+  poly f = imap(r,f);
+  kill r;
+  f = homog(f,h);
+  I = t - f, t*Dt - s;  // defining the Malgrange ideal
+  for (i=1; i<=n; i++)
+  {
+    I = I, D(i)+diff(f,x(i))*Dt;
+  }
+  dbprint(ppl, "starting Groebner basis computation with engine:", whichengine);
+  I = engine(I, whichengine);
+  dbprint(ppl, "finished Groebner basis computation");
+  dbprint(ppl, "I before dehomogenization is" ,I);
+  I = subst(I,h,1); // dehomogenization
+  dbprint(ppl, "I after dehomogenization is" ,I);
+  I = inForm(I,uv); // we are only interested in the initial form w.r.t. uv
+  dbprint(ppl, "the initial ideal:", string(matrix(I)));
   intvec tonselect = 1;
-  for (i=3; i<=2*n+4; i++)
-  {
-    tonselect = tonselect,i;
-  }
-  inF = nselect(inF,tonselect);
-  dbprint(ppl, "generators containing only s:", string(matrix(inF)));
-  inF = engine(inF, whichengine); // is now a principal ideal;
+  for (i=3; i<=2*n+4; i++) { tonselect = tonselect,i; }
+  I = nselect(I,tonselect);
+  dbprint(ppl, "generators containing only s:", string(matrix(I)));
+  I = engine(I, whichengine); // is now a principal ideal;
   setring save;
   ideal J; J[2] = var(1);
   map @m = DDh,J;
-  ideal inF = @m(inF);
-  poly p = inF[1];
+  ideal I = @m(I);
+  poly p = I[1];
   list l = listofroots(p,1);
   return(l);
@@ -1187 +1485 @@
   ring r = 0,(x,y),dp;
   poly f = x^2+y^3+x*y^2;
-  bfctonestep(f);
-  bfctonestep(f,1,1);
-}
-
-proc bfctann (poly f, list #)
-"USAGE:  bfctann(f [,r]);  f a poly, r an optional int
-RETURN:  list of roots of the Bernstein-Sato polynomial bs(f) and their multiplicies
-PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface, defined by f
-NOTE:    In this proc, ann(f^s) is computed.
-@*       If r<>0, @code{std} is used for Groebner basis computations,
+  bfctOneGB(f);
+  bfctOneGB(f,1,1);
+}
+
+proc bfctAnn (poly f, list #)
+"USAGE:  bfctAnn(f [,r]);  f a poly, r an optional int
+RETURN:  list of ideal and intvec
+PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) 
+@*       for the hypersurface defined by f
+ASSUME:  The basering is a commutative polynomial ring in char 0.
+BACKGROUND: In this proc, ann(f^s) is computed and then a system of linear
+@*       equations is solved by linear reductions.
+NOTE:    In the output list, the ideal contains all the roots 
+@*       and the intvec their multiplicities.
+@*       If r<>0, @code{std} is used for GB computations,
 @*       otherwise, and by default, @code{slimgb} is used. 
-@*       If printlevel=1, progress debug messages will be printed,
+DISPLAY: If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
 EXAMPLE: example bfctann; shows examples
@@ -1221 +1524 @@
   ring r = 0,(x,y),dp;
   poly f = x^2+y^3+x*y^2;
-  bfctann(f);
-}
-
-static proc hardexamples ()
+  bfctAnn(f);
+}
+
+/*
+//static proc hardexamples ()
 {
   //  some hard examples
@@ -1246 +1550 @@
   bfct(xyzreiffen45);
 }
-
+*/