Changeset 565e86 in git

Timestamp:

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

Author:

Viktor Levandovskyy <levandov@…>

Branches:

(u'spielwiese', '91fdef05f09f54b8d58d92a472e9c4a43aa4656f')

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

Location:

Singular/LIB

Files:

• bfct.lib (modified) (47 diffs)
• dmod.lib (modified) (45 diffs)
• dmodapp.lib (modified) (43 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);
 }
-
+*/

Singular/LIB/dmod.lib

@@ -1 +1 @@
 //////////////////////////////////////////////////////////////////////////////
-version="$Id: dmod.lib,v 1.32 2008-12-09 16:50:21 levandov Exp $";
+version="$Id: dmod.lib,v 1.33 2008-12-23 21:39:31 levandov Exp $";
 category="Noncommutative";
 info="
@@ -27 +27 @@
 @* OT) the classical Ann F^s algorithm from Oaku and Takayama (J. Pure
         Applied Math., 1999),
-@* LOT) Levandovskyy's modification of the Oaku-Takayama algorithm (unpublished)
+@* LOT) Levandovskyy's modification of the Oaku-Takayama algorithm (ISSAC 2007)
 @* BM) the Ann F^s algorithm by Briancon and Maisonobe (Remarques sur
         l'ideal de Bernstein associe a des polynomes, preprint, 2002)
 
 GUIDE:
-@* - Ann F^s = I = I(F^s) = LD in D(R)[s] can be computed by SannfsBM SannfsOT, SannfsLOT
+@* - Ann F^s = I = I(F^s) = LD in D(R)[s] can be computed by SannfsBM, SannfsOT, SannfsLOT
 @* - Ann^(1) F^s in D(R)[s] can be computed by Sannfslog
 @* - global Bernstein polynomial bs resp. BS in K[s] can be computed by bernsteinBM
-@* - Ann F^s0 = I(F^s0) = LD0 in D(R) can be computed by annfs0, annfsBM, annfsOT, annfsLOT
+@* - Ann F^s0 = I(F^s0) = LD0 in D(R) can be computed by annfs0, annfsBM, annfsOT, annfsLOT, annfs2
 @* - all the relevant data (LD, LD0, bs, PS) are computed by operatorBM
 
@@ -49 +49 @@
 annfsBMI(F[,eng]);      compute Ann F^s and Bernstein ideal for a poly F=f1*..*fP (multivariate algorithm of Briancon-Maisonobe)
 checkRoot(F,a[,S,eng]); check if a given rational is a root of the global Bernstein polynomial of F and compute its multiplicity
-
-SECONDARY PROCEDURES FOR D-MODULES:
-
 annfsBM(F[,eng]);          compute Ann F^s0 in D and Bernstein poly for a poly F (algorithm of Briancon-Maisonobe)
 annfsLOT(F[,eng]);         compute Ann F^s0 in D and Bernstein poly for a poly F (Levandovskyy modification of the Oaku-Takayama algorithm)
 annfsOT(F[,eng]);          compute Ann F^s0 in D and Bernstein poly for a poly F (algorithm of Oaku-Takayama)
 SannfsBM(F[,eng]);         compute Ann F^s in D[s] for a poly F (algorithm of Briancon-Maisonobe)
-SannfsBFCT(F[,eng]);      compute Ann F^s in D[s] for a poly F (algorithm of Briancon-Maisonobe)
+SannfsBFCT(F[,eng]);      compute Ann F^s in D[s] for a poly F (algorithm of Briancon-Maisonobe, other output ordering)
 SannfsLOT(F[,eng]);        compute Ann F^s in D[s] for a poly F (Levandovskyy modification of the Oaku-Takayama algorithm)
 SannfsOT(F[,eng]);         compute Ann F^s in D[s] for a poly F (algorithm of Oaku-Takayama)
-annfs0(I,F [,eng]);           compute Ann F^s0 in D and Bernstein poly from the known Ann F^s in D[s]
+annfs0(I,F [,eng]);          compute Ann F^s0 in D and Bernstein poly from the known Ann F^s in D[s]
+annfs2(I,F [,eng]);           compute Ann F^s0 in D and Bernstein poly from the known Ann F^s in D[s] by using a trick of Noro
+annfsRB(I,F [,eng]);          compute Ann F^s0 in D and Bernstein poly from the known Ann F^s in D[s] by reduceB strategy of Macaulay
 checkRoot1(I,F,a[,eng]);   check whether a rational is a root of the global Bernstein polynomial of F from the known Ann F^s in D[s]
 checkRoot2(I,F,a[,eng]);   check whether a rational is a root of the global Bernstein polynomial of F and compute its multiplicity from the known Ann F^s in D[s]
-checkFactor(I,F,qs[,eng]); check whether a polynomial qs in K[s] is a factor of the global Bernstein polynomial of F from the known Ann F^s in D[s]
+checkFactor(I,F,q[,eng]); check whether a polynomial q in K[s] is a factor of the global Bernstein polynomial of F from the known Ann F^s in D[s]
 
 AUXILIARY PROCEDURES:
@@ -72 +71 @@
 minIntRoot(P,fact); minimal integer root among the roots of a maximal ideal P
 varnum(s);          the number of the variable with the name s
+
 
 SEE ALSO: gmssing_lib, bfct_lib, dmodapp_lib
@@ -92 +92 @@
   "MAIN PROCEDURES:";
   example annfs;
-  example annfs0;
   example Sannfs;
   example Sannfslog;
@@ -101 +100 @@
   example annfsBMI;
   example checkRoot;
+  example annfs0;
+  example annfs2;
+  example annfsRB;
+  example annfs2;
   "SECONDARY D-MOD PROCEDURES:";
   example annfsBM;
@@ -108 +111 @@
   example SannfsLOT;
   example SannfsOT;
+  example SannfsBFCT;
   example checkRoot1;
   example checkRoot2;
@@ -130 +134 @@
 @*       - the ideal LD (which is a Groebner basis) is the needed D-module structure,
 @*       - the list  BS is the list of roots and multiplicities of a Bernstein polynomial of f.
-@*       If @code{printlevel}=1, progress debug messages will be printed,
-@*       if @code{printlevel}>=2, all the debug messages will be printed.
+DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
+@*          if @code{printlevel}>=2, all the debug messages will be printed.
 EXAMPLE: example annfs; shows examples
 "
@@ -252 +256 @@
 @*       In the output ring:
 @*       - the ideal LD is the needed D-module structure.
-@*       If @code{printlevel}=1, progress debug messages will be printed,
-@*       if @code{printlevel}>=2, all the debug messages will be printed.
+DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
+@*          if @code{printlevel}>=2, all the debug messages will be printed.
 EXAMPLE: example Sannfs; shows examples
 "
@@ -365 +369 @@
 @*       If eng <>0, @code{std} is used for Groebner basis computations,
 @*       otherwise, and by default @code{slimgb} is used.
-@*       If printlevel=1, progress debug messages will be printed,
-@*       if printlevel>=2, all the debug messages will be printed.
+DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
+@*          if @code{printlevel}>=2, all the debug messages will be printed.
 EXAMPLE: example Sannfslog; shows examples
 "
@@ -506 +510 @@
 // alternative code for SannfsBM, renamed from annfsBM to ALTannfsBM
 // is superfluos - will not be included in the official documentation
-proc ALTannfsBM (poly F, list #)
+static proc ALTannfsBM (poly F, list #)
 "USAGE:  ALTannfsBM(f [,eng]);  f a poly, eng an optional int
 RETURN:  ring
@@ -514 +518 @@
 @*       If eng <>0, @code{std} is used for Groebner basis computations,
 @*       otherwise, and by default @code{slimgb} is used.
-@*       If printlevel=1, progress debug messages will be printed,
-@*       if printlevel>=2, all the debug messages will be printed.
+DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
+@*          if @code{printlevel}>=2, all the debug messages will be printed.
 EXAMPLE: example ALTannfsBM; shows examples
 "
@@ -703 +707 @@
 NOTE:    If eng <>0, @code{std} is used for Groebner basis computations,
 @*       otherwise, and by default @code{slimgb} is used.
-@*       If printlevel=1, progress debug messages will be printed,
-@*       if printlevel>=2, all the debug messages will be printed.
+DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
+@*          if @code{printlevel}>=2, all the debug messages will be printed.
 EXAMPLE: example bernsteinBM; shows examples
 "
@@ -933 +937 @@
 @*       If eng <>0, @code{std} is used for Groebner basis computations,
 @*       otherwise, and by default @code{slimgb} is used.
-@*       If printlevel=1, progress debug messages will be printed,
-@*       if printlevel>=2, all the debug messages will be printed.
+DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
+@*          if @code{printlevel}>=2, all the debug messages will be printed.
 EXAMPLE: example annfsBM; shows examples
 "
@@ -1217 +1221 @@
 
 
-// try to replace s with -s-1 => data is shorter
+// replacing s with -s-1 => data is shorter
 // analogue of annfs0
 proc annfs2(ideal I, poly F, list #)
@@ -1223 +1227 @@
 RETURN:  ring
 PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, based on the
-output of procedures SannfsBM, SannfsOT or SannfsLOT
+output of Sannfs-like procedure
 NOTE:    activate this ring with the @code{setring} command. In this ring,
 @*       - the ideal LD (which is a Groebner basis) is the annihilator of f^s,
@@ -1229 +1233 @@
 @*       If eng <>0, @code{std} is used for Groebner basis computations,
 @*       otherwise and by default @code{slimgb} is used.
-@*       Uses the shorter form of expressions in the variable s (the idea of Noro).
-@*       If printlevel=1, progress debug messages will be printed,
-@*       if printlevel>=2, all the debug messages will be printed.
+@*       annfs2 uses the shorter form of expressions in the variable s (the idea of Noro).
+DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
+@*       if @code{printlevel}>=2, all the debug messages will be printed.
 EXAMPLE: example annfs2; shows examples
 "
@@ -1390 +1394 @@
 @*       If eng <>0, @code{std} is used for Groebner basis computations,
 @*       otherwise and by default @code{slimgb} is used.
-@*       Uses the shorter form of expressions in the variable s (the idea of Noro).
-@*       If printlevel=1, progress debug messages will be printed,
-@*       if printlevel>=2, all the debug messages will be printed.
+@*       This procedure follows the 'reduceB' strategy, used in Macaulay2.
+DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
+@*       if @code{printlevel}>=2, all the debug messages will be printed.
 EXAMPLE: example annfsRB; shows examples
 "
@@ -1561 +1565 @@
   poly F = x^3+y^3+z^3;
   printlevel = 0;
-  def A = SannfsBM(F);
-  setring A;
-  LD;
+  def A = SannfsBM(F); setring A;
+  LD; // s-parametric ahhinilator
   poly F = imap(r,F);
-  def B  = annfsRB(LD,F);
-  setring B;
+  def B  = annfsRB(LD,F); setring B;
   LD;
   BS;
@@ -1584 +1586 @@
 @*       If eng <>0, @code{std} is used for Groebner basis computations,
 @*       otherwise and by default @code{slimgb} is used.
-@*       If printlevel=1, progress debug messages will be printed,
-@*       if printlevel>=2, all the debug messages will be printed.
+DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
+@*       if @code{printlevel}>=2, all the debug messages will be printed.
 EXAMPLE: example operatorBM; shows examples
 "
@@ -1889 +1891 @@
 // for this, do special modulo with emphasis on the 2comp
 // of the transp matrix, there must appear 1 in the GB
+// then use lift to get the expression...
 // need: (c,<) ordering for such comp's
 
-// proc operatorModulo(poly F, ideal I, poly b)
-// "USAGE:  operatorModulo(f,I,b);  f a poly, I an ideal, b a poly
-// RETURN:  poly
-// PURPOSE: compute the B-operator from the poly F, ideal I = Ann f^s and Bernstein-Sato
-// polynomial b using modulo
-// NOTE:  The computations take place in the ring, similar to the one returned by Sannfs procedure.
-// @*       If printlevel=1, progress debug messages will be printed,
-// @*       if printlevel>=2, all the debug messages will be printed.
-// EXAMPLE: example operatorModulo; shows examples
-// "
-// {
-//   // with hom_kernel;
-//   matrix AA[1][2] = F,-b;
-//   //  matrix M = matrix(subst(I,s,s+1)*freemodule(2)); // ann f^{s+1}
-//   matrix N = matrix(I); // ann f^s
-//   //  matrix K = hom_kernel(AA,M,N);
-//   matrix K = modulo(AA,N);
-//   K = transpose(K);
-//   print(K);
-//   ideal J;
-//   int i;
-//   poly t; number n;
-//   for(i=1; i<=nrows(K); i++)
-//   {
-//     if (K[i,2]!=0)
-//     {
-//       if ( leadmonom(K[i,2]) == 1)
-//       {
-//      t = K[i,1];
-//      n = leadcoef(K[i,2]);
-//      t = t/n;
-//      //        J = J, K[i][2];
-//      break;
-//       }
-//     }
-//   }
-//   ideal J = groebner(subst(I,s,s+1)); // for NF
-//   t = NF(t,J);
-//   "candidate:"; t;
-//   J = subst(J,s,s-1);
-//   // test:
-//   if ( NF(t*F-b,J) !=0)
-//   {
-//     "Problem: PS does not work on F";
-//   }
-//   return(t);
-// }
-// example
-// {
-//   "EXAMPLE:"; echo = 2;
-//   //  ring r = 0,(x,y,z),Dp;
-//   //   poly F = x^3+y^3+z^3;
-//   LIB "dmod.lib"; option(prot); option(mem);
-//   ring r = 0,(x,y),Dp;
-//   //  poly F = x^3+y^3+x*y^2;
-//   poly F = x^4 + y^4 + x*y^4;
-//   def A = Sannfs(F); // here we get LD = ann f^s
-//   setring A;
-//   poly F = imap(r,F);
-//   def B = annfs0(LD,F); // to obtain BS polynomial
-//   list BS = imap(B,BS);
-//   poly b = fl2poly(BS,"s");
-//   LD = groebner(LD);
-//   poly PS = operatorModulo(F,LD,b);
-//   PS; // the operator, s.t. PS*F^{s+1} = bs*F^s mod LD
-//   reduce(PS*F-bs,LD); // check the property of PS
-// }
-
-proc fl2poly(list L, string s)
-"USAGE:  fl2poly(L,s); L a list, s a string
+/*
+proc operatorModulo(poly F, ideal I, poly b)
+"USAGE:  operatorModulo(f,I,b);  f a poly, I an ideal, b a poly
 RETURN:  poly
-PURPOSE: reconstruct a monic polynomial in one variable from its factorization
-ASSUME:  s is a string with the name of some variable and L is supposed to consist of two entries:
-@*         L[1] of the type ideal with the roots of a polynomial
-@*         L[2] of the type intvec with the multiplicities of corr. roots
-EXAMPLE: example fl2poly; shows examples
+PURPOSE: compute the B-operator from the poly F, ideal I = Ann f^s and Bernstein-Sato
+polynomial b using modulo
+NOTE:  The computations take place in the ring, similar to the one returned by Sannfs procedure.
+@*       If printlevel=1, progress debug messages will be printed,
+@*       if printlevel>=2, all the debug messages will be printed.
+EXAMPLE: example operatorModulo; shows examples
 "
 {
-  if (varnum(s)==0)
-  {
-    ERROR("no such variable found"); return(0);
-  }
-  poly x = var(varnum(s));
-  poly P = 1;
-  int sl = size(L[1]);
-  ideal RR = L[1];
-  intvec IV = L[2];
-  for(int i=1; i<= sl; i++)
-  {
-    P = P*((x-RR[i])^IV[i]);
-  }
-  return(P);
+  // with hom_kernel;
+  matrix AA[1][2] = F,-b;
+  //  matrix M = matrix(subst(I,s,s+1)*freemodule(2)); // ann f^{s+1}
+  matrix N = matrix(I); // ann f^s
+  //  matrix K = hom_kernel(AA,M,N);
+  matrix K = modulo(AA,N);
+  K = transpose(K);
+  print(K);
+  ideal J;
+  int i;
+  poly t; number n;
+  for(i=1; i<=nrows(K); i++)
+  {
+    if (K[i,2]!=0)
+    {
+      if ( leadmonom(K[i,2]) == 1)
+      {
+        t = K[i,1];
+        n = leadcoef(K[i,2]);
+        t = t/n;
+        //        J = J, K[i][2];
+        break;
+      }
+    }
+  }
+  ideal J = groebner(subst(I,s,s+1)); // for NF
+  t = NF(t,J);
+  "candidate:"; t;
+  J = subst(J,s,s-1);
+  // test:
+  if ( NF(t*F-b,J) !=0)
+  {
+    "Problem: PS does not work on F";
+  }
+  return(t);
 }
 example
 {
   "EXAMPLE:"; echo = 2;
-  ring r = 0,(x,y,z,s),Dp;
-  ideal I = -1,-4/3,-5/3,-2;
-  intvec mI = 2,1,1,1;
-  list BS = I,mI;
-  poly p = fl2poly(BS,"s");
-  p;
-  factorize(p,2);
-}
+  //  ring r = 0,(x,y,z),Dp;
+  //   poly F = x^3+y^3+z^3;
+  LIB "dmod.lib"; option(prot); option(mem);
+  ring r = 0,(x,y),Dp;
+  //  poly F = x^3+y^3+x*y^2;
+  poly F = x^4 + y^4 + x*y^4;
+  def A = Sannfs(F); // here we get LD = ann f^s
+  setring A;
+  poly F = imap(r,F);
+  def B = annfs0(LD,F); // to obtain BS polynomial
+  list BS = imap(B,BS);
+  poly b = fl2poly(BS,"s");
+  LD = groebner(LD);
+  poly PS = operatorModulo(F,LD,b);
+  PS; // the operator, s.t. PS*F^{s+1} = bs*F^s mod LD
+  reduce(PS*F-bs,LD); // check the property of PS
+}
+*/
+
+
 
 proc annfsParamBM (poly F, list #)
@@ -2005 +1975 @@
 @*       If eng <>0, @code{std} is used for Groebner basis computations,
 @*       otherwise, and by default @code{slimgb} is used.
-@*       If printlevel=1, progress debug messages will be printed,
-@*       if printlevel>=2, all the debug messages will be printed.
+DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
+@*          if @code{printlevel}>=2, all the debug messages will be printed.
 EXAMPLE: example annfsParamBM; shows examples
 "
@@ -3445 +3415 @@
 @*       If eng <>0, @code{std} is used for Groebner basis 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 SannfsBFCT; shows examples
@@ -3637 +3607 @@
   poly F = x^3+y^3+z^3*w;
   printlevel = 0;
-  def A = SannfsBFCT(F);
-  setring A;
+  def A = SannfsBFCT(F); setring A;
   intvec v = 1,2,3,4,5,6,7,8;
   // are there polynomials, depending on s only?
   nselect(LD,v);
   // a fancier example:
-  def R = reiffen(4,5);
-  setring R;
+  def R = reiffen(4,5); setring R;
   v = 1,2,3,4;
+  RC; // the Reiffen curve in 4,5
   def B = SannfsBFCT(RC);
   setring B;
@@ -3885 +3854 @@
 }
 
-
+/*
 proc SannfsLOTold(poly F, list #)
 "USAGE:  SannfsLOT(f [,eng]);  f a poly, eng an optional int
@@ -4109 +4078 @@
 }
 
+*/
 
 proc annfsLOT(poly F, list #)
@@ -4158 +4128 @@
 RETURN:  ring
 PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, based on the
-output of procedures SannfsBM, SannfsOT or SannfsLOT
+output of Sannfs-like procedure
 NOTE:    activate this ring with the @code{setring} command. In this ring,
 @*       - the ideal LD (which is a Groebner basis) is the annihilator of f^s,
@@ -4299 +4269 @@
   poly F = x^3+y^3+z^3;
   printlevel = 0;
-  def A = SannfsBM(F);
+  def A = SannfsBM(F);   setring A;
   // alternatively, one can use SannfsOT or SannfsLOT
-  setring A;
   LD;
   poly F = imap(r,F);
-  def B  = annfs0(LD,F);
-  setring B;
+  def B  = annfs0(LD,F);  setring B;
   LD;
   BS;
@@ -4656 +4624 @@
 @*       We compute the real annihilator for any rational value of n (both generic and exceptional).
 @*       The implementation goes along the lines of Saito-Sturmfels-Takayama, Alg. 5.3.15
-@*       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 annfspecial; shows examples
@@ -4724 +4692 @@
 }
 
-static proc new_ex_annfspecial()
+/*
+//static proc new_ex_annfspecial()
 {
   //another example for annfspecial: x3+y3+z3
@@ -4740 +4709 @@
   annfspecial(LD,F,-1,1);   // exceptional root
 }
+*/
 
 proc minIntRoot(ideal P, int fact)
@@ -4817 +4787 @@
 "USAGE:  isHolonomic(M); M an ideal/module/matrix
 RETURN:  int, 1 if M is holonomic and 0 otherwise
+ASSUME: basering is a Weyl algebra
 PURPOSE: check the modules for the property of holonomy
 NOTE:    M is holonomic, if 2*dim(M) = dim(R), where R is a
@@ -4956 +4927 @@
 @*       'alg1' (default) - for the algorithm 1 of J. Martin-Morales (unpublished)
 @*       'alg2' - for the algorithm 2 of J. Martin-Morales (unpublished)
-@*       The output int is:
-@*       - if the algorithm 1 is chosen: 1 if -alpha is a root of the global Bernstein polynomial and 0 otherwise
-@*       - if the algorithm 2 is chosen: the multiplicity of -alpha as a root of the global Bernstein polynomial of f.
-@*                                       (If -alpha is not a root, the output is 0)
+@*       The output of type int is:
+@*       'alg1': 1 if -alpha is a root of the global Bernstein polynomial and 0 otherwise
+@*       'alg2':  the multiplicity of -alpha as a root of the global Bernstein polynomial of f.
+@*                (If -alpha is not a root, the output is 0)
 @*       If eng <>0, @code{std} is used for Groebner basis computations,
 @*       otherwise (and by default) @code{slimgb} is used.
-@*       If printlevel=1, progress debug messages will be printed,
-@*       if printlevel>=2, all the 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 checkRoot; shows examples
 "
@@ -5068 +5039 @@
 NOTE:    If eng <>0, @code{std} is used for Groebner basis computations,
 @*       otherwise (and by default) @code{slimgb} is used.
-@*       If printlevel=1, progress debug messages will be printed,
-@*       if printlevel>=2, all the 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 checkRoot1; shows examples
 "
@@ -5192 +5163 @@
 NOTE:    If eng <>0, @code{std} is used for Groebner basis computations,
 @*       otherwise (and by default) @code{slimgb} is used.
-@*       If printlevel=1, progress debug messages will be printed,
-@*       if printlevel>=2, all the 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 checkRoot2; shows examples
 "
@@ -5324 +5295 @@
 NOTE:    If eng <>0, @code{std} is used for Groebner basis computations,
 @*       otherwise (and by default) @code{slimgb} is used.
-@*       If printlevel=1, progress debug messages will be printed,
-@*       if printlevel>=2, all the 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 checkFactor; shows examples
 "
@@ -5461 +5432 @@
 }
 
-
-static proc exCheckGenericity()
+/// ****** EXAMPLES ************
+
+/*
+
+//static proc exCheckGenericity()
 {
   LIB "control.lib";
@@ -5484 +5458 @@
 }
 
-static proc exOT_17()
+//static proc exOT_17()
 {
   // Oaku-Takayama, p.208
@@ -5499 +5473 @@
 }
 
-static proc exOT_16()
+//static proc exOT_16()
 {
   // Oaku-Takayama, p.208
@@ -5513 +5487 @@
 }
 
-static proc ex_bcheck()
+//static proc ex_bcheck()
 {
   ring R = 0,(x,y),dp;
@@ -5526 +5500 @@
 }
 
-static proc ex_bcheck2()
+//static proc ex_bcheck2()
 {
   ring R = 0,(x,y),dp;
@@ -5536 +5510 @@
 }
 
-static proc ex_BMI()
+//static proc ex_BMI()
 {
   // a hard example
@@ -5549 +5523 @@
 }
 
-static proc ex2_BMI()
+//static proc ex2_BMI()
 {
   // this example was believed to be intractable in 2005 by Gago-Vargas, Castro and Ucha
@@ -5563 +5537 @@
 }
 
-static proc ex_operatorBM()
+//static proc ex_operatorBM()
 {
   ring r = 0,(x,y,z,w),Dp;
@@ -5578 +5552 @@
   reduce(PS*F-bs,LD); // check the property of PS
 }
+
+*/

Singular/LIB/dmodapp.lib

@@ -1 +1 @@
 //////////////////////////////////////////////////////////////////////////////
-version="$Id: dmodapp.lib,v 1.12 2008-12-09 17:52:35 levandov Exp $";
+version="$Id: dmodapp.lib,v 1.13 2008-12-23 21:39:31 levandov Exp $";
 category="Noncommutative";
 info="
 LIBRARY: dmodapp.lib     Applications of algebraic D-modules
-AUTHORS: Viktor Levandovskyy,     levandov@math.rwth-aachen.de
+AUTHORS: Viktor Levandovskyy,      levandov@math.rwth-aachen.de
 @*             Daniel Andres, daniel.andres@math.rwth-aachen.de
 
@@ -18 +18 @@
 rational function G/F from Quot(R). These can be computed via annPoly resp. annRat.
 
-@* - initial form and initial ideals in Weyl algebras with respect to a given weight vector can be computed with the help of inForm, initialmalgrange, initialideal.
+@* - initial form and initial ideals in Weyl algebras with respect to a given weight vector can be computed with the help of inForm, initialMalgrange, initialIdeal.
 
 @* - appelF1, appelF2 and appelF4 return ideals in parametric Weyl algebra, which
@@ -32 +32 @@
 DLoc0(I, F);  presentation of the localization of D/I w.r.t. f^s, based on the procedure SDLoc
 
-initialmalgrange(f[,s,t,u,v]); Groebner basis of the initial Malgrange ideal for a given poly
-initialideal(I,u,v[,s,t]); initial ideal of a given ideal w.r.t. given weights
+initialMalgrange(f[,s,t,u,v]); Groebner basis of the initial Malgrange ideal for a given poly
+initialIdeal(I,u,v[,s,t]); initial ideal of a given ideal w.r.t. given weights
 inForm(f,w);     initial form of a poly/ideal w.r.t. a given weight
 isFsat(I, F);       check whether the ideal I is F-saturated
@@ -44 +44 @@
 engine(I,i);     computes a Groebner basis with the algorithm given by i
 poly2list(f);   decompose a poly to a list of terms and exponents
+fl2poly(L,s);  reconstruct a monic univariate polynomial from its factorization
+insertGenerator(id,p[,k]); insert an element into an ideal/module
+deleteGenerator(id,k); delete the k-th element from an ideal/module
+
 
 SEE ALSO: dmod_lib, gmssing_lib, bfct_lib
@@ -58 +62 @@
 // charInfo(); ???
 
+
+///////////////////////////////////////////////////////////////////////////////
+// testing for consistency of the library:
 proc testdmodapp()
 {
-  example initialideal;
-  example initialmalgrange;
+  example initialIdeal;
+  example initialMalgrange;
   example DLoc;
   example DLoc0;
@@ -73 +80 @@
   example appelF4;
   example poly2list;
-}
-
-
-proc initialideal (ideal I, intvec u, intvec v, list #)
-"USAGE:  initialideal(I,u,v [,s,t]);  I an ideal, u,v intvecs, s,t optional ints
-RETURN:  an ideal, a Broebner basis of the initial ideal of the input ideal w.r.t. the weights u and v
-PURPOSE: compute the initial ideal
-NOTE:    Assume, I is an ideal in the n-th Weyl algebra where the sequence of the
-@*       indeterminates is x(1),...,x(n),D(1),...,D(n).
-@*       Further assume that u is the weight for the x(i) and v the weight for the D(i).
+  example fl2poly;
+  example insertGenerator;
+  example deleteGenerator;
+}
+
+
+proc initialIdeal (ideal I, intvec u, intvec v, list #)
+"USAGE:  initialIdeal(I,u,v [,s,t]);  I ideal, u,v intvecs, s,t optional ints
+RETURN:  ideal, GB of initial ideal of the input ideal wrt the weights u and v
+PURPOSE: computes the initial ideal
+NOTE:    Assume, I is an ideal in the n-th Weyl algebra D, where the sequence
+@*       of the indeterminates is x(1),...,x(n),D(1),...,D(n).
+@*       Further assume that u is the weight for the x(i) and v the weight for
+@*       the D(i).
 @*       Note that the returned ideal is not a D-ideal but an ideal in the associated
 @*       graded ring while the Groebner basis is a subset of D.
@@ -89 +100 @@
 @*       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,
+DISPLAY: If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
-EXAMPLE: example initialideal; shows examples
+EXAMPLE: example initialIdeal; shows examples
 "
 {
@@ -113 +124 @@
     }
   }
-  def D = initialidealengine("initialideal", whichengine, methodord, I, u, v);
-  ideal inF = fetch(D,inF);
+  def D = initialIdealEngine("initialIdeal", whichengine, methodord, I, u, v);
+  ideal inF = fetch(D,inF); attrib(inF,"isSB",1);
   return(inF);
 }
@@ -125 +136 @@
   intvec u = -1; intvec v = 2;
   ideal I = x^2*Dx^2,x*Dx^4;
-  ideal J = initialideal(I,u,v); J;
-}
-
-proc initialmalgrange (poly f,list #)
-"USAGE:  initialmalgrange(f, [,s,t,u,v]);  f a poly, s,t,u optional ints, v an optional intvec
-RETURN:  a ring, the Weyl algebra induced by the basering, extended in the indeterminates t and Dt
-PURPOSE: compute the initial Malgrange ideal of a given poly w.r.t. the weight vector
-@*       (-1,0...,0,1,0,...,0) such that the weight of t is -1 and the weight of Dt is 1.
+  ideal J = initialIdeal(I,u,v); J;
+}
+
+proc initialMalgrange (poly f,list #)
+"USAGE:  initialMalgrange(f, [,s,t,u,v]);  f poly, s,t,u optional ints, v optional intvec
+RETURN:  ring, the Weyl algebra induced by basering, extended with vars t and Dt,
+@*       containing the ideal \"inF\", being the initial ideal of the Malgrange
+@*       ideal of f.
+PURPOSE: computes the initial Malgrange ideal of a given poly wrt the weight
+@*       vector (-1,0...,0,1,0,...,0) such that the weight of t is -1 and the
+@*       weight of Dt is 1.
 NOTE:    Activate the output ring with the @code{setring} command.
 @*       Varnames of the basering should not include t and Dt.
-@*       In the ouput ring, inF is the initial Malgrange ideal.
-@*       If s<>0, @code{std} is used for Groebner basis computations, otherwise,
-@*        and by default, @code{slimgb} is used.
+@*       If s<>0, @code{std} is used for Groebner basis 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 u<>0, the order of the variables is reversed.
-@*       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 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 initialmalgrange; shows examples
+EXAMPLE: example initialMalgrange; shows examples
 "
 {
@@ -175 +188 @@
         if (size(#)>3)
         {
-          if (typeof(#[4])=="intvec" && size(#[4])==n && ispositive(#[4])==1)
+          if (typeof(#[4])=="intvec" && size(#[4])==n && allPositive(#[4])==1)
           {
             u0 = #[4];
@@ -187 +200 @@
     u0 = 1:n;
   }
-  def D = initialidealengine("initialmalgrange",whichengine, methodord, f, 0, 0, u0, reversevars);
-  setring D;
+  def D = initialIdealEngine("initialMalgrange",whichengine, methodord, f, 0, 0, u0, reversevars);
+  setring save;
   return(D);
 }
@@ -196 +209 @@
   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;
   setring r;
   intvec v = 3,2;
-  def D2 = initialmalgrange(f,1,0,1,v);
+  def D2 = initialMalgrange(f,1,0,1,v);
   setring D2;
   inF;
 }
 
-proc initialidealengine(string calledfrom, int whichengine, int blockord, list #)
+static proc initialIdealEngine(string calledfrom, int whichengine, int blockord, list #)
 {
   //#[1] = f or I
@@ -218 +231 @@
   n = nvars(save);
   intvec uv;
-  if (calledfrom == "initialideal")
+  if (calledfrom == "initialIdeal")
   {
     ideal I = #[1];
@@ -227 +240 @@
     noofvars = 2*n+1;
   }
-  else
+  else // initialMalgrange
   {
     poly f = #[1];
-    if (calledfrom == "initialmalgrange")
-    {
-      uv[n+2] = 1;
-      noofvars = 2*n+3;
-      intvec u0 = #[4];
-      int reversevars = #[5];
-      ring r = 0,(x(1..n)),wp(u0);
-    }
-    else // bfctonestep
-    {
-      uv[n+3] = 1;
-      noofvars = 2*n+4;
-      ring r = 0,(x(1..n)),dp;
-    }
+    uv[n+2] = 1;
+    noofvars = 2*n+3;
+    intvec u0 = #[4];
+    int reversevars = #[5];
+    ring r = 0,(x(1..n)),wp(u0);
     poly f = fetch(save,f);
     uv[1] = -1; uv[noofvars] = 0;
@@ -249 +253 @@
   //  for the ordering
   intvec @a;
-  if (calledfrom == "initialmalgrange")
+  if (calledfrom == "initialMalgrange")
   {
     int d = deg(f);
@@ -264 +268 @@
     @a = weighttx,weightD,1;
   }
-  else
+  else // initialIdeal
   {
     @a = 1:noofvars;
-    if (calledfrom == "bfctonestep")
-    {
-      @a[2] = 2;
-      intvec @a2 = @a; @a2[2] = 0; @a2[2*n+4] = 0;
-    }
   }
   if (blockord == 0) // default: blockordering
   {
-    if (calledfrom == "initialmalgrange")
+    if (calledfrom == "initialMalgrange")
     {
       ring Dh = 0,(t,x(n..1),Dt,D(n..1),h),(a(@a),a(uv),dp(noofvars-1),lp(1));
     }
-    else
-    {
-      if (calledfrom == "initialideal")
-      {
-        ring Dh = 0,(x(1..n),D(1..n),h),(a(@a),dp(noofvars-1),lp(1));
-      }
-      else // bfctonestep
-      {
-        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 // initialIdeal
+    {
+      ring Dh = 0,(x(1..n),D(1..n),h),(a(@a),dp(noofvars-1),lp(1));
     }
   }
@@ -302 +294 @@
     dbprint(ppl,"weights for ordering:",transpose(@a));
     dbprint(ppl,"the ordering matrix:",@Ord);
-    if (calledfrom == "initialmalgrange")
+    if (calledfrom == "initialMalgrange")
     {
       ring Dh = 0,(t,x(n..1),Dt,D(n..1),h),(a(@a),M(@Ord));
     }
-    else
-    {
-      if (calledfrom == "initialideal")
-      {
-        ring Dh = 0,(x(1..n),D(1..n),h),(a(@a),M(@Ord));
-      }
-      else // bfctonestep
-      {
-        ring Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),M(@Ord));
-      }
+    else // initialIdeal
+    {
+      ring Dh = 0,(x(1..n),D(1..n),h),(a(@a),M(@Ord));
     }
   }
@@ -321 +306 @@
   // non-commutative relations
   matrix @relD[noofvars][noofvars];
-  if (calledfrom == "initialmalgrange")
-  {
-    for (i=1; i<=n+1; i++)
-    {
-      @relD[i,n+1+i] = h^(d+1);
-    }
-  }
-  else
-  {
-    if (calledfrom == "initialideal")
-    {
-      for (i=1; i<=n; i++)
-      {
-        @relD[i,n+i] = h^2;
-      }
-    }
-    else // bfctonestep
-    {
-      @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;
-      }
+  if (calledfrom == "initialMalgrange")
+  {
+    for (i=1; i<=n+1; i++) { @relD[i,n+1+i] = h^(d+1); }
+  }
+  else // initialIdeal
+  {
+    for (i=1; i<=n; i++)
+    {
+      @relD[i,n+i] = h^2;
     }
   }
@@ -353 +322 @@
   dbprint(ppl,"computing in ring",DDh);
   ideal I;
-  if (calledfrom == "initialideal")
+  if (calledfrom == "initialIdeal")
   {
     I = fetch(save,I);
@@ -367 +336 @@
     {
       I = I, D(i)+diff(f,x(i))*Dt;
-    }
-    if (calledfrom == "bfctonestep")
-    {
-      I = I,t*Dt-s;
     }
   }
@@ -380 +345 @@
   dbprint(ppl, "I after dehomogenization is" ,I);
   I = inForm(I,uv); // we are only interested in the initial form w.r.t. uv
-  if (calledfrom == "initialmalgrange")
+  if (calledfrom == "initialMalgrange")
   {
     // keep the names of the variables of the basering
@@ -425 +390 @@
     }
   }
-  else
-  {
-    if (calledfrom == "initialideal")
-    {
-      ring @D = 0,(x(1..n),D(1..n)),dp;
-      def D = Weyl(@D);
-      setring D;
-      ideal I = imap(DDh,I);
-      kill Dh,DDh;
-    }
-  }
-  if (calledfrom <> "bfctonestep")
-  {
-    dbprint(ppl, "starting cosmetic Groebner basis computation with engine:", whichengine);
-    I = engine(I, whichengine);
-    dbprint(ppl,"finished cosmetic Groebner basis computation:");
-    dbprint(ppl,"the initial ideal is:", I);
-  }
-  ideal inF = I;
+  else // initialIdeal
+  {
+    ring @D = 0,(x(1..n),D(1..n)),dp;
+    def D = Weyl(@D);
+    setring D;
+    ideal I = imap(DDh,I);
+    kill Dh,DDh;
+  }
+  dbprint(ppl, "starting cosmetic Groebner basis computation with engine:", whichengine);
+  I = engine(I, whichengine);
+  dbprint(ppl,"finished cosmetic Groebner basis computation:");
+  dbprint(ppl,"the initial ideal is:", I);
+  ideal inF = I; attrib(inF,"isSB",1);
   export(inF);
   return(basering);
 }
 
-proc inForm (list #)
-"USAGE:  inForm(I,w);  I an ideal, w an intvec
-RETURN:  the initial form of I w.r.t. the weight vector w
-PURPOSE: compute the initial form of an  ideal w.r.t a given weight vector
-NOTE:  the size of the weight vector must be equal to the number of variables of the basering
+proc inForm (ideal I, intvec w)
+"USAGE:  inForm(I,w);  I ideal, w intvec
+RETURN:  the initial form of I wrt the weight vector w
+PURPOSE: computes the initial form of an ideal wrt a given weight vector
+NOTE:  the size of the weight vector must be equal to the number of variables
+@*     of the basering.
 EXAMPLE: example inForm; shows examples
 "
 {
-  if (size(#)<2)
-  {
-    ERROR("inForm needs two arguments of type poly/ideal,intvec");
-  }
-  if (typeof(#[1])=="poly" || typeof(#[1])=="ideal")
-  {
-    ideal I = #[1];
-  }
-  else
-  {
-    ERROR("first argument must be of type poly or ideal");
-  }
-  if (typeof(#[2])=="intvec")
-  {
-    def w = #[2];
-  }
-  else
-  {
-    ERROR("second argument must be of type intvec");
-  }
   if (size(w) != nvars(basering))
   {
     ERROR("weight vector has wrong dimension");
   }
+  if (I == 0)
+  {
+    return(I);
+  }  
   int j,i,s,m;
   list l;
@@ -487 +431 @@
   {
     l = poly2list(I[j]);
-    m = scalarprod(w,l[1][1]);
-    g = 0;
+    m = scalarProd(w,l[1][1]);
+    g = l[1][2];
     for (i=2; i<=size(l); i++)
     {
-      s = scalarprod(w,l[i][1]);
-      m = Max(list(s,m));
-    }
-    for (i=1; i<=size(l); i++)
-    {
-      if (scalarprod(w,l[i][1]) == m)
+      s = scalarProd(w,l[i][1]);
+      if (s == m)
       {
-        g = g+l[i][2];
+        g = g + l[i][2];
       }
+      else
+      {
+        if (s > m)
+        {
+          m = s;
+          g = l[i][2];
+        }
+      }
     }
     J[j] = g;
   }
-  if (typeof(#[1])=="poly")
-  {
-    return(J[1]); // if the input was a poly, return a poly
-  }
-  else
-  {
-    return(J);    // otherwise, return an ideal
-  }
+  return(J);
 }
 example
@@ -527 +468 @@
   inForm(I,w2);
   inForm(I,w3);
-  inForm(F,w1);
-}
-
-static proc charVariety(ideal I)
+}
+
+/*
+
+proc charVariety(ideal I)
 "USAGE:  charVariety(I);  I an ideal
 RETURN:  ring
@@ -615 +557 @@
 }
 
+/*
+
 // TODO
-static proc charInfo(ideal I)
+
+/*
+proc charInfo(ideal I)
 "USAGE:  charInfo(I);  I an ideal
 RETURN:  ring
@@ -638 +584 @@
   // run primdec
 }
+*/
 
 
@@ -666 +613 @@
   "EXAMPLE:"; echo = 2;
   ring r = 0,x,dp;
-  def A = appelF1(); //(1,2,3,4);
+  def A = appelF1(); 
   setring A;
   IAppel1;
@@ -696 +643 @@
   "EXAMPLE:"; echo = 2;
   ring r = 0,x,dp;
-  def A = appelF2(); //(1,2,3,4);
+  def A = appelF2(); 
   setring A;
   IAppel2;
@@ -726 +673 @@
   "EXAMPLE:"; echo = 2;
   ring r = 0,x,dp;
-  def A = appelF4(); //(1,2,3,4);
+  def A = appelF4(); 
   setring A;
   IAppel4;
@@ -740 +687 @@
   list l;
   int i = 1;
-  if (f==0) // just for the zero polynomial
+  if (f == 0) // just for the zero polynomial
   {
     l[1] = list(leadexp(f), lead(f));
   }
-  while (f!=0)
+  else { l[size(f)] = list(); } // memory pre-allocation
+  while (f != 0)
   {
     l[i] = list(leadexp(f), lead(f));
-    f = f-lead(f);
+    f = f - lead(f);
     i++;
   }
@@ -803 +751 @@
 "USAGE:  DLoc(I, F);  I an ideal, F a poly
 RETURN: nothing (exports objects instead)
-ASSUME: the basering is a Weyl algebra and I is F-saturated
+ASSUME: the basering is a Weyl algebra
 PURPOSE: compute the presentation of the localization of D/I w.r.t. f^s
 NOTE:    In the basering, the following objects are exported:
@@ -860 +808 @@
 {
   /* assume: to be run in the output ring of SDLoc */
-  /* todo: add F, eliminate vars*Dvars, factorize BS */
+  /* doing: add F, eliminate vars*Dvars, factorize BS */
   /* analogue to annfs0 */
   def @R2 = basering;
@@ -1292 +1240 @@
 NOTE:    activate the output ring with the @code{setring} command.
 @*       In the output ring:
-@*       - the ideal RLD (which is given in a Groebner basis) is the annihilator.
+@*       - the ideal LD (which is given in a Groebner basis) is the annihilator.
 @*       If @code{printlevel}=1, progress debug messages will be printed,
 @*       if @code{printlevel}>=2, all the debug messages will be printed.
@@ -1369 +1317 @@
   kill @R1;   kill @R2;
   dbprint(ppl,"// -4-[annRat] running modulo");
-  ideal RLD  = modulo(G,LL);
+  ideal LD  = modulo(G,LL);
   dbprint(ppl,"// -4-[annRat] running GB on the final result");
-  RLD  = engine(RLD,0);
-  export RLD;
+  LD  = engine(LD,0);
+  export LD;
   return(@R3);
 }
@@ -1382 +1330 @@
   def B = annRat(g,f);
   setring B;
-  RLD;
+  LD;
   // Now, compare with the output of Macaulay2:
   ideal tst = 3*x*Dx + 2*y*Dy + 1, y^3*Dy^2 - x^2*Dy^2 + 6*y^2*Dy + 6*y, 
 9*y^2*Dx^2*Dy - 4*y*Dy^3 + 27*y*Dx^2 + 2*Dy^2, 9*y^3*Dx^2 - 4*y^2*Dy^2 + 10*y*Dy -10; 
  option(redSB); option(redTail);
- RLD = groebner(RLD);
+ LD = groebner(LD);
  tst = groebner(tst);
- print(matrix(NF(RLD,tst)));  print(matrix(NF(tst,RLD)));
+ print(matrix(NF(LD,tst)));  print(matrix(NF(tst,LD)));
  // So, these two answers are the same
 }
 
-static proc ex_annRat()
-{
-  // more complicated example
+/*
+//static proc ex_annRat()
+{
+  // more complicated example for annRat
   ring r = 0,(x,y,z),dp;
   poly f = x3+y3+z3; // mir = -2
@@ -1402 +1351 @@
   setring A;
 }
+*/
 
 proc annPoly(poly f)
@@ -1409 +1359 @@
 NOTE:    activate the output ring with the @code{setring} command.
 @*       In the output ring:
-@*       - the ideal RLD (which is given in a Groebner basis) is the annihilator.
+@*       - the ideal LD (which is given in a Groebner basis) is the annihilator.
 @*       If @code{printlevel}=1, progress debug messages will be printed,
 @*       if @code{printlevel}>=2, all the debug messages will be printed.
@@ -1455 +1405 @@
 }
 
-static proc exCusp()
+/* DIFFERENT EXAMPLES 
+
+//static proc exCusp()
 {
   "EXAMPLE:"; echo = 2;
@@ -1474 +1426 @@
 }
 
-static proc exWalther1()
+//static proc exWalther1()
 {
   // p.18 Rem 3.10
@@ -1494 +1446 @@
 }
 
-static proc exWalther2()
+//static proc exWalther2()
 {
   // p.19 Rem 3.10 cont'd
@@ -1518 +1470 @@
 }
 
-static proc exWalther3()
+//static proc exWalther3()
 {
   // can check with annFs too :-)
@@ -1550 +1502 @@
 }
 
+*/
+
 proc engine(ideal I, int i) 
 "USAGE:  engine(I,i);  I an ideal, i an int
@@ -1581 +1535 @@
   engine(I,1); // uses std
 }
+
+proc insertGenerator (list #)
+"USAGE:  insertGenerator(id,p[,k]);  id an ideal/module, p a poly/vector, k an optional int
+RETURN:  same as id
+PURPOSE: insert an element into an ideal or a module
+NOTE:    If k is given, p is inserted at position k, otherwise (and by default),
+@*       p is inserted at the beginning.
+EXAMPLE: example insertGenerator; shows examples
+"
+{
+  if (size(#) < 2)
+  {
+    ERROR("insertGenerator has to be called with at least 2 arguments (ideal/module,poly/vector)");
+  }
+  string inp1 = typeof(#[1]);
+  if (inp1 == "ideal" || inp1 == "module")
+  {
+    if (inp1 == "ideal") { ideal id = #[1]; }
+    else { module id = #[1]; }
+  }
+  else { ERROR("first argument has to be of type ideal or module"); }
+  string inp2 = typeof(#[2]);
+  if (inp2 == "poly" || inp2 == "vector")
+  {
+    if (inp2 =="poly") { poly f = #[2]; }
+    else
+    {
+      if (inp1 == "ideal") 
+      {
+        ERROR("second argument has to be a poly if first argument is an ideal");
+      }
+      else { vector f = #[2]; }
+    }
+  }
+  else { ERROR("second argument has to be of type poly/vector"); }
+  int n = ncols(id);
+  int k = 1; // default
+  if (size(#)>=3)
+  {
+    if (typeof(#[3]) == "int")
+    {
+      k = #[3];
+      if (k<=0)
+      {
+        ERROR("third argument has to be positive");
+      }
+    }
+    else { ERROR("third argument has to be of type int"); }
+  }
+  execute(inp1 +" J;");
+  if (k == 1) { J = f,id; }
+  else
+  {
+    if (k>n)
+    {
+      J = id;
+      J[k] = f;
+    }
+    else // 1<k<=n
+    {
+      J[1..k-1] = id[1..k-1];
+      J[k] = f;
+      J[k+1..n+1] = id[k..n];
+    }
+  }
+  return(J);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  ring r = 0,(x,y,z),dp;
+  ideal I = x^2,z^4;
+  insertGenerator(I,y^3);
+  insertGenerator(I,y^3,2);
+  module M = I;
+  insertGenerator(M,[x^3,y^2,z],2);
+}
+
+proc deleteGenerator (list #)
+"USAGE:  deleteGenerator(id,k);  id an ideal/module, k an int
+RETURN:  same as id
+PURPOSE: deletes the k-th element from an ideal or a module
+EXAMPLE: example insertGenerator; shows examples
+"
+{
+  if (size(#) < 2)
+  {
+    ERROR("deleteGenerator has to be called with 2 arguments (ideal/module,int)");
+  }
+  string inp1 = typeof(#[1]);
+  if (inp1 == "ideal" || inp1 == "module")
+  {
+    if (inp1 == "ideal") { ideal id = #[1]; }
+    else { module id = #[1]; }
+  }
+  else { ERROR("first argument has to be of type ideal or module"); }
+  string inp2 = typeof(#[2]);
+  if (inp2 == "int" || inp2 == "number") { int k = int(#[2]); }
+  else { ERROR("second argument has to be of type int"); }
+  int n = ncols(id);
+  if (n == 1) { ERROR(inp1+" must have more than one generator"); }
+  if (k<=0 || k>n) { ERROR("second argument has to be in the range 1,...,"+string(n)); }
+  execute(inp1 +" J;");
+  if (k == 1) { J = id[2..n]; }
+  else
+  {
+    if (k == n) { J = id[1..n-1]; }
+    else
+    {
+      J[1..k-1] = id[1..k-1];
+      J[k..n-1] = id[k+1..n];
+    }
+  }
+  return(J);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  ring r = 0,(x,y,z),dp;
+  ideal I = x^2,y^3,z^4;
+  deleteGenerator(I,2);
+  module M = [x,y,z],[x2,y2,z2],[x3,y3,z3];;
+  deleteGenerator(M,2);
+}
+
+proc fl2poly(list L, string s)
+"USAGE:  fl2poly(L,s); L a list, s a string
+RETURN:  poly
+PURPOSE: reconstruct a monic polynomial in one variable from its factorization
+ASSUME:  s is a string with the name of some variable and L is supposed to consist of two entries:
+@*         L[1] of the type ideal with the roots of a polynomial
+@*         L[2] of the type intvec with the multiplicities of corr. roots
+EXAMPLE: example fl2poly; shows examples
+"
+{
+  if (varnum(s)==0)
+  {
+    ERROR("no such variable found"); return(0);
+  }
+  poly x = var(varnum(s));
+  poly P = 1;
+  int sl = size(L[1]);
+  ideal RR = L[1];
+  intvec IV = L[2];
+  for(int i=1; i<= sl; i++)
+  {
+    P = P*((x-RR[i])^IV[i]);
+  }
+  return(P);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  ring r = 0,(x,y,z,s),Dp;
+  ideal I = -1,-4/3,-5/3,-2;
+  intvec mI = 2,1,1,1;
+  list BS = I,mI;
+  poly p = fl2poly(BS,"s");
+  p;
+  factorize(p,2);
+}