Index: ngular/LIB/bfct.lib
===================================================================
--- Singular/LIB/bfct.lib	(revision d499d16596daf2413557e9454d6d76712a4cc5e8)
+++ 	(revision )
@@ -1,1557 +1,0 @@
-//////////////////////////////////////////////////////////////////////////////
-version="$Id: bfct.lib,v 1.13 2009-01-16 12:53:09 Singular Exp $";
-category="Noncommutative";
-info="
-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 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, Sturmfels, 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]);            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
-linReduceIdeal(I,[s,t]);    reduces generators of ideal by linear reductions only
-linSyzSolve(I[,s]);         computes a linear dependency of the elements of ideal I
-
-AUXILIARY PROCEDURES:
-
-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
-";
-
-
-LIB "qhmoduli.lib"; // for Max
-LIB "dmod.lib";     // for SannfsBFCT etc
-LIB "dmodapp.lib";  // for initialIdealW etc
-
-
-///////////////////////////////////////////////////////////////////////////////
-// testing for consistency of the library:
-proc testbfctlib ()
-{
-  // tests all procs for consistency
-  "AUXILIARY PROCEDURES:";
-  example allPositive;
-  example scalarProd;
-  example vec2poly;
-  "MAIN PROCEDURES:";
-  example bfct;
-  example bfctSyz;
-  example bfctAnn;
-  example bfctOneGB;
-  example bfctIdeal;
-  example pIntersect;
-  example pIntersectSyz;
-  example linReduce;
-  example linReduceIdeal;
-  example linSyzSolve;
-}
-
-//--------------- auxiliary procedures ---------------------------------------------------------
-
-static proc gradedWeyl (intvec u,intvec v)
-"USAGE:  gradedWeyl(u,v); u,v intvecs
-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).
-@*       u+v has to be a non-negative intvec. 
-"
-{
-  int i;
-  def save = basering;
-  int n = nvars(save)/2;
-  if (nrows(u)<>n || nrows(v)<>n)
-  {
-    ERROR("weight vectors have wrong dimension");
-  } 
-  intvec uv,gr;
-  uv = u+v;
-  for (i=1; i<=n; i++)
-  {
-    if (uv[i]>=0)
-    {
-      if (uv[i]==0)
-      {
-	gr[i] = 0;
-      }
-      else
-      {
-	gr[i] = 1;
-      }
-    }
-    else
-    {
-      ERROR("the sum of the weight vectors has to be a non-negative intvec");
-    }
-  }
-  list l = ringlist(save);
-  list l2 = l[2];
-  matrix l6 = l[6];
-  for (i=1; i<=n; i++)
-  {
-    if (gr[i] == 1)  
-    {
-       l2[n+i] = "xi("+string(i)+")";
-       l6[i,n+i] = 0;
-    }
-  }
-  l[2] = l2;
-  l[6] = l6;
-  def G = ring(l);
-  return(G);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  ring @D = 0,(x,y,z,Dx,Dy,Dz),dp;
-  def D = Weyl();
-  setring D;
-  intvec u = -1,-1,1; intvec v = 2,1,1;
-  def G = gradedWeyl(u,v);
-  setring G; G;
-}
-
-
-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 allPositive; shows an example
-"
-{
-  int i;
-  for (i=1; i<=size(v); i++)
-  {
-    if (v[i]<=0)
-    {
-      return(0);
-      break;
-    }
-  }
-  return(1);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  intvec v = 1,2,3;
-  allPositive(v);
-  intvec w = 1,-2,3;
-  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 findFirst; shows an example
-"
-{
-  int j;
-  for (j=1; j<=size(l); j++)
-  {
-    if (l[j]==i)
-    {
-      return(j);
-      break;
-    }
-  }
-  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
-{
-  "EXAMPLE:"; echo = 2;
-  ring r = 0,x,dp;
-  list l = 1,2,3;
-  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
-"
-{
-  int i; int sp;
-  if (size(v)!=size(w))
-  {
-    ERROR("non-matching dimensions");
-  }
-  else
-  {
-    for (i=1; i<=size(v);i++)
-    {
-      sp = sp + v[i]*w[i];
-    }
-  }
-  return(sp);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  intvec v = 1,2,3;
-  intvec w = 4,5,6;
-  scalarProd(v,w);
-}
-
-//-------------- main procedures -------------------------------------------------------
-
-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 I by its elements
-PURPOSE: reduce a list of polys 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 redtail        = 0; // 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];
-      }
-    }
-  }
-  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];
-	}
-      }
-    }
-  }
-  int i,j,k;
-  int sI = ncols(I);
-  // 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;
-  // === 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])
-	{
-	  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);
-	    if (remembercoeffs <> 0)
-	    {
-	      v = v - c * M[i];
-	    }
-	  }
-	}
-      }
-    }
-  }
-  if (remembercoeffs <> 0)
-  {
-    list l = f,v;
-    return(l);
-  }
-  else
-  {
-    return(f);
-  }
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  ring r = 0,(x,y),dp;
-  ideal I = 1,y,xy;
-  poly f = 5xy+7y+3;
-  poly g = 7x+5y+3;
-  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.
-DISPLAY: If printlevel=1, progress debug messages will be printed,
-@*       if printlevel>=2, all the debug messages will be printed.
-EXAMPLE: example linSyzSolve; shows examples
-"
-{
-  int whichengine     = 0; // default
-  int enginespecified = 0; // default
-  if (size(#)>0)
-  {
-    if (typeof(#[1])=="int" || typeof(#[1])=="number")
-    {
-      whichengine = int( #[1]);
-      enginespecified = 1;
-    }
-  }
-  int ppl = printlevel - voice +2;
-  int sI = ncols(I);
-  // check if we are done
-  if (I[sI]==0)
-  {
-    vector v = gen(sI);
-  }
-  else
-  {
-    // ------- 1. introduce undefined coeffs ------------------
-    def save = basering;
-    int p = char(save);
-    if (enginespecified == 0)
-    {
-      if (p <> 0)
-      {
-	whichengine = 1;
-      }
-    }
-    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 ---
-    int i;   poly W;  ideal QQ;
-    for (i=1; i<=sI; i++)
-    {
-      W = W + @a(i)*I[i];
-    }
-    while (W!=0)
-    {
-      QQ = QQ,leadcoef(W);
-      W = W - lead(W);
-    }
-    // QQ consists of polynomial expressions in @a(i) of type number
-    setring @A;
-    ideal QQ = imap(@B,QQ);
-    // ------- 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 linSyzSolve");
-      // output zeroes
-      setring save;
-      kill @A,@aA,@B;
-      return(v);
-    }
-    // ------- 4. in order to get the numeric values -------
-    matrix AA = matrix(maxideal(1));
-    module MQQ = std(module(QQ));
-    AA = NF(AA,MQQ); // todo: we still receive NF warnings
-    dbprint(ppl, "AA after NF:",AA);
-    //    "AA after NF:"; print(AA);
-    for(i=1; i<=sI; i++)
-    {
-      AA = subst(AA,var(i),1);
-    }
-    dbprint(ppl, "AA after subst:",AA);
-    //    "AA after subst: "; print(AA);
-    vector v = (module(transpose(AA)))[1];
-    setring save;
-    vector v = imap(@A,v);
-    kill @A,@aA,@B;
-  }
-  return(v);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  ring r = 0,x,dp;
-  ideal I = x,2x;
-  linSyzSolve(I);
-  ideal J = x,x2;
-  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.
-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
-"
-{
-  // assume I is given in Groebner basis
-  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+1;
-  // ---case 1: I = basering---
-  if (size(I) == 1)
-  {
-    if (simplify(I,2)[1] == 1)
-    {
-      return(gen(2)); // = s
-    }
-  }
-  def save = basering;
-  int n = nvars(save);
-  int i,j,k;
-  // ---case 2: intersection is zero---
-  intvec degs = leadexp(s);
-  intvec possdegbounds;
-  list degI;
-  i = 1;
-  while (i <= ncols(I))
-  {
-    if (i == ncols(I)+1)
-    {
-      break;
-    }
-    degI[i] = leadexp(I[i]);
-    for (j=1; j<=n; j++)
-    {
-      if (degs[j] == 0)
-      {
-	if (degI[i][j] <> 0)
-	{
-	  break;
-	}
-      }
-      if (j == n)
-      {
-	k++;
-	possdegbounds[k] = Max(degI[i]);
-      }
-    }
-    i++;
-  }
-  int degbound = Min(possdegbounds);
-  dbprint(ppl,"a lower bound for the degree of the insection is:");
-  dbprint(ppl,degbound);
-  if (degbound == 0) // lm(s) does not appear in lm(I)
-  {
-    return(vector(0));
-  }
-  // ---case 3: intersection is non-trivial---
-  ideal redNI = 1;
-  vector v;
-  list l,ll;
-  l[1] = vector(0);
-  poly toNF,tobracket,newNF,rednewNF,oldNF,secNF;
-  i = 1;
-  dbprint(ppl+1,"pIntersect starts...");
-  while (1)
-  {
-    dbprint(ppl,"testing degree: "+string(i));
-    if (i>1)
-    {
-      oldNF = newNF;
-      tobracket = s^(i-1) - oldNF;
-      if (tobracket==0) // todo bug in bracket?
-      {
-	toNF = 0;
-      }
-      else
-      {
-	toNF = bracket(tobracket,secNF);
-      }
-      newNF = NF(toNF+oldNF*secNF,I);  // = NF(s^i,I)
-    }
-    else
-    {
-      newNF = NF(s,I);
-      secNF = newNF;
-    }
-    ll = linReduce(newNF,redNI,1);
-    rednewNF = ll[1];
-    l[i+1] = ll[2];
-    dbprint(ppl,"newNF is:", newNF);
-    dbprint(ppl,"rednewNF is:", rednewNF);
-    if (rednewNF != 0) // no linear dependency
-    {
-      redNI[i+1] = rednewNF;
-      i++;
-    }
-    else // there is a linear dependency, hence we are done
-    {
-      dbprint(ppl+1,"the degree of the generator of the intersection is:", i);
-      break;
-    }
-  }
-  dbprint(ppl,"used linear reductions:", l);
-  // we obtain the coefficients of the generator of the intersection by the used reductions:
-  ring @R = 0,(a(1..i+1)),dp;
-  setring @R;
-  list l = imap(save,l);
-  ideal C;
-  for (j=1;j<=i+1;j++)
-  {
-    C[j] = 0;
-    for (k=1;k<=j;k++)
-    {
-      C[j] = C[j]+l[j][k]*a(k);
-    }
-  }
-  for (j=i;j>=1;j--)
-  {
-    C[i+1]= subst(C[i+1],a(j),a(j)+C[j]);
-  }
-  matrix m = coeffs(C[i+1],maxideal(1));
-  vector v = gen(i+1);
-  for (j=1;j<=i+1;j++)
-  {
-    v = v + m[j,1]*gen(j);
-  }
-  setring save; 
-  v = imap(@R,v);
-  kill @R;
-  dbprint(ppl+1,"pIntersect finished");
-  return(v);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  ring r = 0,(x,y),dp;
-  poly f = x^2+y^3+x*y^2;
-  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:  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.
-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
-"
-{
-  // assume I is given in Groebner basis
-  ideal I = II;
-  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
-  int modengine    = 1; // default
-  int solveincharp = 0; // default
-  def save = basering;
-  if (size(#)>0)
-  {
-    if (typeof(#[1])=="int" || typeof(#[1])=="number")
-    {
-      solveincharp = int(#[1]);
-    }
-    if (size(#)>1)
-    {
-      if (typeof(#[2])=="int" || typeof(#[2])=="number")
-      {
-        whichengine = int(#[2]);
-      }
-      if (size(#)>2)
-      {
-	if (typeof(#[3])=="int" || typeof(#[3])=="number")
-	{
-	  modengine = int(#[3]);
-	}
-      }
-    }
-  }
-  int i,j;
-  vector v;
-  poly tobracket,toNF,newNF,p;
-  ideal NI = 1;
-  newNF = NF(s,I);
-  NI[2] = newNF;
-  if (solveincharp<>0)
-  {
-    list l = ringlist(save);
-    l[1] = solveincharp;
-    matrix l5 = l[5];
-    matrix l6 = l[6];
-    def @Rp = ring(l);
-    setring @Rp;
-    list l = ringlist(@Rp);
-    l[5] = fetch(save,l5);
-    l[6] = fetch(save,l6);
-    def Rp = ring(l);
-    setring Rp;
-    kill @Rp;
-    dbprint(ppl+1,"solving in ring ", Rp);
-    vector v;
-    map phi = save,maxideal(1);
-    poly s = phi(s);
-    ideal NI = 1;
-    setring save;
-  }
-  i = 1;
-  dbprint(ppl+1,"pIntersectSyz starts...");
-  dbprint(ppl+1,"with ideal I=", I);
-  while (1)
-  {
-    dbprint(ppl,"i:"+string(i));
-    if (i>1)
-    {
-      tobracket = s^(i-1)-NI[i];
-      if (tobracket!=0)
-      {
-	toNF = bracket(tobracket,NI[2]) + NI[i]*NI[2];
-      }
-      else
-      {
-	toNF = NI[i]*NI[2];
-      }
-      newNF =  NF(toNF,I);
-      NI[i+1] = newNF;
-    }
-    // look for a solution
-    dbprint(ppl,"linSyzSolve starts with: "+string(matrix(NI)));
-    if (solveincharp<>0) // modular method
-    {
-      setring Rp;
-      NI[i+1] = phi(newNF);
-      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 = linSyzSolve(NI,whichengine);
-	if (v==0)
-	{
-	  break;
-	}
-      }
-      else // no modular solution
-      {
-	setring save;
-	v = 0;
-      }
-    }
-    else // non-modular method
-    {
-      v = linSyzSolve(NI,whichengine);
-    }
-    matrix MM[1][nrows(v)] = matrix(v);
-    dbprint(ppl,"linSyzSolve ready  with: "+string(MM));
-    kill MM;
-    //  "linSyzSolve ready with"; print(v);
-    if (v!=0)
-    {
-      // a solution:
-      //check for the reality of the solution
-      p = 0;
-      for (j=1; j<=i+1; j++)
-      {
-        p = p + v[j]*NI[j];
-      }
-      if (p!=0)
-      {
-        dbprint(ppl,"linSyzSolve: bad solution!");
-      }
-      else
-      {
-        dbprint(ppl,"linSyzSolve: got solution!");
-	// "got solution!";
-        break;
-      }
-    }
-    // no solution:
-    i++;
-  }
-  dbprint(ppl+1,"pIntersectSyz finished");
-  return(v);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  ring r = 0,(x,y),dp;
-  poly f = x^2+y^3+x*y^2;
-  def D = initialMalgrange(f);
-  setring D;
-  inF;
-  poly s = t*Dt;
-  pIntersectSyz(s,inF);
-  int p = prime(20000);
-  pIntersectSyz(s,inF,p,0,0);
-}
-
-proc vec2poly (list #)
-"USAGE:  vec2poly(v [,i]);  v a vector or an intvec, i an optional int
-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
-"
-{
-  def save = basering;
-  int i,ringvar;
-  ringvar = 1; // default
-  if (size(#) > 0)
-  {
-    if (typeof(#[1])=="vector" || typeof(#[1])=="intvec")
-    {
-      def v = #[1];
-    }
-    else
-    {
-      ERROR("wrong input: expected vector/intvec expression");
-    }
-    if (size(#) > 1)
-    {
-      if (typeof(#[2])=="int" || typeof(#[2])=="number")
-      {
-	ringvar = int(#[2]);
-      }
-    }
-  }
-  if (ringvar > nvars(save))
-  {
-    ERROR("var out of range");
-  }
-  poly p;
-  for (i=1; i<=nrows(v); i++)
-  {
-    p = p + v[i]*(var(ringvar))^(i-1);
-  }
-  return(p);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  ring r = 0,(x,y),dp;
-  vector v = gen(1) + 3*gen(3) + 22/9*gen(4);
-  intvec iv = 3,2,1;
-  vec2poly(v,2);
-  vec2poly(iv);
-}
-
-static proc listofroots (list #)
-{
-  def save = basering;
-  int n = nvars(save);
-  int i;
-  poly p;
-  if (typeof(#[1])=="vector")
-  {
-    vector b = #[1];
-    for (i=1; i<=nrows(b); i++)
-    {
-      p = p + b[i]*(var(1))^(i-1);
-    }
-  }
-  else
-  {
-    p = #[1];
-  }
-  int substitution = int(#[2]);
-  ring S = 0,s,dp;
-  ideal J;
-  for (i=1; i<=n; i++)
-  {
-    J[i] = s;
-  }
-  map @m = save,J;
-  poly p = @m(p);
-  if (substitution == 1)
-  {
-    p = subst(p,s,-s-1);
-  }
-  // the rest of this proc is nicked from bernsteinBM from dmod.lib
-  list P = factorize(p);//with constants and multiplicities
-  ideal bs; intvec m;   //the Bernstein polynomial is monic, so we are not interested in constants
-  for (i=2; i<= size(P[1]); i++)  //we delete P[1][1] and P[2][1]
-  {
-    bs[i-1] = P[1][i];
-    m[i-1]  = P[2][i];
-  }
-  bs =  normalize(bs);
-  bs = -subst(bs,s,0);
-  setring save;
-  ideal bs = imap(S,bs);
-  kill S;
-  list BS = bs,m;
-  return(BS);
-}
-
-static proc bfctengine (poly f, int inorann, int whichengine, int methodord, int methodpIntersect, int pIntersectchar, int modengine, intvec u0)
-{
-  int ppl = printlevel - voice +2;
-  int i;
-  def save = basering;
-  int n = nvars(save);
-  if (size(variables(f)) == 0) // f is constant
-  {
-    return(list(ideal(0),intvec(0)));
-  }  
-  if (inorann == 0) // bfct using initial ideal
-  {
-    def D = initialMalgrange(f,whichengine,methodord,1,u0);
-    setring D;
-    ideal J = inF;
-    kill inF;
-    poly s = t*Dt;
-  }
-  else // bfct using Ann(f^s)
-  {
-    def D = SannfsBFCT(f,whichengine);
-    setring D;
-    ideal J = LD;
-    kill LD;
-  }
-  vector b;
-  // try it modular
-  if (methodpIntersect <> 0) // pIntersectSyz
-  {
-    if (pIntersectchar == 0) // pIntersectSyz::modular
-    {
-      int lb = 30000;
-      int ub = 536870909;
-      i = 1;
-      list usedprimes;
-      while (b == 0)
-      {
-	dbprint(ppl,"number of run in the loop: "+string(i));
-	int q = prime(random(lb,ub));
-	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);
-	}
-	i++;
-      }
-    }
-    else // pIntersectSyz::non-modular 
-    {
-      b = pIntersectSyz(s,J,0,whichengine);
-    }
-  }
-  else // pIntersect: linReduce
-  {
-    b = pIntersect(s,J);
-  }
-  setring save;
-  vector b = imap(D,b);
-  if (inorann == 0)
-  {
-    list l = listofroots(b,1);
-  }
-  else
-  {
-    list l = listofroots(b,0);
-  }
-  return(l);
-}
-
-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 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, 
-@*       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
-"
-{
-  int ppl = printlevel - voice +2;
-  int i;
-  int n = nvars(basering);
-  // in # we have two switches:
-  // one for the engine used for Groebner basis computations,
-  // one for  M() ordering or its realization
-  // in # can also be the optional weight vector
-  int whichengine  = 0; // default
-  int methodord    = 0; // default
-  intvec u0        = 0; // default
-  if (size(#)>0)
-  {
-    if (typeof(#[1])=="int" || typeof(#[1])=="number")
-    {
-      whichengine = int(#[1]);
-    }
-    if (size(#)>1)
-    {
-      if (typeof(#[2])=="int" || typeof(#[2])=="number")
-      {
-        methodord = int(#[2]);
-      }
-      if (size(#)>2)
-      {
-	if (typeof(#[3])=="intvec" && size(#[3])==n && allPositive(#[3])==1)
-        {
-	  u0 = #[3];
-	}
-      }
-    }
-  }
-  list b = bfctengine(f,0,whichengine,methodord,0,0,0,u0);
-  return(b);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  ring r = 0,(x,y),dp;
-  poly f = x^2+y^3+x*y^2;
-  bfct(f);
-  intvec v = 3,2;
-  bfct(f,1,0,v);
-}
-
-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 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 bfctSyz; shows examples
-"
-{
-  int ppl = printlevel - voice +2;
-  int i;
-  // in # we have four switches:
-  // one for the engine used for Groebner basis computations in char 0,
-  // one for  M() ordering or its realization
-  // one for a modular method when computing the intersection
-  // and one for the engine used for Groebner basis computations in char >0
-  // in # can also be the optional weight vector
-  def save = basering;
-  int n = nvars(save);
-  int whichengine = 0; // default
-  int methodord   = 0; // default
-  int pIntersectchar  = 0; // default
-  int modengine   = 1; // default
-  intvec u0       = 0; // default
-  if (size(#)>0)
-  {
-    if (typeof(#[1])=="int" || typeof(#[1])=="number")
-    {
-      whichengine = int(#[1]);
-    }
-    if (size(#)>1)
-    {
-      if (typeof(#[2])=="int" || typeof(#[2])=="number")
-      {
-        methodord = int(#[2]);
-      }
-      if (size(#)>2)
-      {
-        if (typeof(#[3])=="int" || typeof(#[3])=="number")
-        {
-	  pIntersectchar = int(#[3]);
-	}
-	if (size(#)>3)
-        {
-	  if (typeof(#[4])=="int" || typeof(#[4])=="number")
-          {
-	    modengine = int(#[4]);
-	  }
-	  if (size(#)>4)
-	  {
-	    if (typeof(#[5])=="intvec" && size(#[5])==n && allPositive(#[5])==1)
-	    {
-	      u0 = #[5];
-	    }
-	  }
-	}
-      }
-    }
-  }
-  list b = bfctengine(f,0,whichengine,methodord,1,pIntersectchar,modengine,u0);
-  return(b);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  ring r = 0,(x,y),dp;
-  poly f = x^2+y^3+x*y^2;
-  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 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 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
-"
-{
-  int ppl = printlevel - voice +2;
-  int i;
-  def save = basering;
-  int n = nvars(save)/2;
-  int whichengine = 0; // default
-  int methodord   = 0; // default
-  if (size(#)>0)
-  {
-    if (typeof(#[1])=="int" || typeof(#[1])=="number")
-    {
-      whichengine = int(#[1]);
-    }
-    if (size(#)>1)
-    {
-      if (typeof(#[2])=="int" || typeof(#[2])=="number")
-      {
-        methodord = int(#[2]);
-      }
-    }
-  }
-  ideal J = initialIdealW(I,-w,w,whichengine,methodord);
-  poly s;
-  for (i=1; i<=n; i++)
-  {
-    s = s + w[i]*var(i)*var(n+i);
-  }
-  vector b = pIntersect(s,J);
-  list l = listofroots(b,0);
-  return(l);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  ring @D = 0,(x,y,Dx,Dy),dp;
-  def D = Weyl();
-  setring D;
-  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 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 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 bfctOneGB; shows examples
-"
-{
-  int ppl = printlevel - voice +2;
-  def save = basering;
-  int n = nvars(save);
-  int noofvars = 2*n+4;
-  int i;
-  int whichengine = 0; // default
-  int methodord   = 0; // default
-  if (size(#)>0)
-  {
-    if (typeof(#[1])=="int" || typeof(#[1])=="number")
-    {
-      whichengine = int(#[1]);
-    }
-    if (size(#)>1)
-    {
-      if (typeof(#[2])=="int" || typeof(#[2])=="number")
-      {
-        methodord = int(#[2]);
-      }
-    }
-  }
-  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,"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; }
-  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 I = @m(I);
-  poly p = I[1];
-  list l = listofroots(p,1);
-  return(l);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  ring r = 0,(x,y),dp;
-  poly f = x^2+y^3+x*y^2;
-  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. 
-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
-"
-{
-  def save = basering;
-  int ppl = printlevel - voice + 2;
-  int whichengine = 0; // default
-  if (size(#)>0)
-  {
-    if (typeof(#[1])=="int" || typeof(#[1])=="number")
-    {
-      whichengine = int(#[1]);
-    }
-  }
-  list b = bfctengine(f,1,whichengine,0,1,0,0,0);
-  return(b);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  ring r = 0,(x,y),dp;
-  poly f = x^2+y^3+x*y^2;
-  bfctAnn(f);
-}
-
-/*
-//static proc hardexamples ()
-{
-  //  some hard examples
-  ring r1 = 0,(x,y,z,w),dp;
-  // ab34
-  poly ab34 = (z3+w4)*(3z2x+4w3y);
-  bfct(ab34);
-  // ha3
-  poly ha3 = xyzw*(x+y)*(x+z)*(x+w)*(y+z+w);
-  bfct(ha3);
-  // ha4
-  poly ha4 = xyzw*(x+y)*(x+z)*(x+w)*(y+z)*(y+w);
-  bfct(ha4);
-  // chal4: reiffen(4,5)*reiffen(5,4)
-  ring r2 = 0,(x,y),dp;
-  poly chal4 = (x4+xy4+y5)*(x5+x4y+y4);
-  bfct(chal4);
-  // (xy+z)*reiffen(4,5)
-  ring r3 = 0,(x,y,z),dp;
-  poly xyzreiffen45 = (xy+z)*(y4+yz4+z5);
-  bfct(xyzreiffen45);
-}
-*/
