Index: Singular/LIB/dmodvar.lib
===================================================================
--- Singular/LIB/dmodvar.lib	(revision 85688982c19c69995eb98b9b499e6f6517c6ede1)
+++ Singular/LIB/dmodvar.lib	(revision 85688982c19c69995eb98b9b499e6f6517c6ede1)
@@ -0,0 +1,843 @@
+////////////////////////////////////////////////////////////////////////
+version="$Id: dmod.lib,v 1.22 2007/11/27 11:08:51 Singular Exp $";
+category="Noncommutative";
+info="
+LIBRARY: dmodvar.lib     Algebraic D-modules for varieties
+AUTHORS: Daniel Andres,          daniel.andres@math.rwth-aachen.de
+@*       Viktor Levandovskyy,    levandov@math.rwth-aachen.de
+@*       Jorge Martin-Morales,   jorge@unizar.es
+
+THEORY: Let K be a field of characteristic 0. Given a polynomial ring
+@*      R = K[x_1,...,x_n] and a set of polynomial f_1,..., f_r in R, define
+@*      F = f_1 * ... * f_r and F^s:=f_1^s_1*...*f_r^s_r for symbolic discrete
+@*      (that is shiftable) variables s_1,..., s_r.
+@*      The module R[1/F]*F^s has a structure of a D<S>-module, where 
+@*      D<S> := D(R) tensored with S over K, where
+@*      - D(R) is an n-th Weyl algebra K<x_1,...,x_n,d_1,...,d_n | d_j x_j = x_j d_j +1>
+@*      - S is the universal enveloping algebra of gl_r, generated by s_{ij}, where s_{ii}=s_i.
+@* One is interested in the following data:
+@* - the left ideal Ann F^s in D<S>, usually denoted by LD in the output
+@* - global Bernstein polynomial in one variable s = s_1 + ...+ s_r, denoted by bs,
+@* - its minimal integer root s0, the list of all roots of bs, which are known
+@*   to be rational, with their multiplicities, which is denoted by BS
+@* - an r-tuple of operators in D<S>, denoted by PS, such that the functional equality
+@*     sum(k=1 to k=r) P_k*f_k*F^s = bs*F^s holds in R[1/F]*F^s.
+
+REFERENCES:
+@*      (BMS06) Budur, Mustata, Saito: Bernstein-Sato polynomials of arbitrary varieties (2006).
+@*      (ALM09) Andres, Levandovskyy, Martin-Morales : Principal Intersection and Bernstein-Sato Polynomial of an Affine Variety (2009).
+
+MAIN PROCEDURES:
+bfctVarIn(F[,L]);     compute the roots of the Bernstein-Sato polynomial b(s) of the variety V(F) using initial ideal approach
+bfctVarAnn(F[,L]);  compute the roots of the Bernstein-Sato polynomial b(s) of the variety V(F) using Sannfs approach
+SannfsVar(F[,O,e]); compute the annihilator of F^s in the ring D<S>
+
+AUXILIARY PROCEDURES:
+makeIF(F[,ORD]);    create the Malgrange ideal, associated with F = F[1],..,F[P]
+";
+
+// Static procs:
+//coDim(I);           compute the codimension of the leading ideal of I
+// dmodvarAssumeViolation()
+// ORDstr2list (ORD, NN)
+// smallGenCoDim(I,k)
+
+
+LIB "bfun.lib";    // for pIntersect
+LIB "dmodapp.lib"; // for isCommutative etc.
+
+
+///////////////////////////////////////////////////////////////////////////////
+
+// testing for consistency of the library:
+proc testdmodvarlib ()
+{
+  "AUXILIARY PROCEDURES:";
+  example makeIF;
+  "MAIN PROCEDURES:";
+  example bfctVarIn;
+  example bfctVarAnn;
+  example SannfsVar;
+}
+
+//   example coDim;
+
+///////////////////////////////////////////////////////////////////////////////
+
+static proc dmodvarAssumeViolation()
+{
+  // returns Boolean : yes/no [for assume violation]
+  // char K = 0
+  // no qring
+  if (  (size(ideal(basering)) >0) || (char(basering) >0) )
+  {
+    //    "ERROR: no qring is allowed";
+    return(1);
+  }
+  return(0);
+}
+
+// da: in smallGenCoDim(), rewritten using mstd business
+static proc coDim (ideal I)
+"USAGE:  coDim (I);  I an ideal
+RETURN:  int
+PURPOSE: computes the codimension of the ideal generated by the leading monomials
+@*       of the given generators of the ideal. This is also the codimension of
+@*       the ideal if it is represented by a standard basis.
+NOTE:    The codimension of an ideal I means the number of variables minus the
+@*       Krull dimension of the basering modulo I.
+EXAMPLE: example SannfsVar; shows examples
+"
+{
+  int n = nvars(basering);
+  int d = dim(I); // to insert: check whether I is in GB
+  return(n-d);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  ring R = 0,(x,y,z),Dp;
+  ideal I = x^2+y^3, z;
+  coDim(std(I));
+}
+
+static proc ORDstr2list (string ORD, int NN)
+{
+  /* convert an ordering defined in NN variables in the  */
+  /* string form into the same ordering in the list form */
+  string st;
+  st = "ring @Z = 0,z(1.." + string(NN) + "),";
+  st = st + ORD + ";";
+  execute(st); kill st;
+  list L = ringlist(@Z)[3];
+  kill @Z;
+  return(L);
+}
+
+proc SannfsVar (ideal F, list #)
+"USAGE:  SannfsVar(F [,ORD,eng]);  F an ideal, ORD an optional string, eng an optional int
+RETURN:  ring
+PURPOSE: compute the D<S>-module structure of D<S>*f^s where f = F[1]*..*F[P]
+and D<S> is the Weyl algebra D tensored with K<S>=U(gl_P), according to the
+generalized algorithm by Briancon and Maisonobe for affine varieties.
+NOTE:    activate this ring with the @code{setring} command.
+@*       In the ring D<S>, the ideal LD is the needed D<S>-module structure.
+@*       The value of ORD must be an elimination ordering in D<Dt,S> for Dt
+@*       written in the string form, otherwise the result may have no meaning.
+@*       By default ORD = '(a(1..(P)..1),a(1..(P+P^2)..1),dp)'.
+@*       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.
+EXAMPLE: example SannfsVar; shows examples
+"
+{
+  if (dmodvarAssumeViolation())
+  {
+    ERROR("Basering is inappropriate: characteristic>0 or qring present");
+  }
+  if (!isCommutative())
+  {
+    ERROR("Basering must be commutative");
+  }
+  def save = basering;
+  int N = nvars(basering);
+  int P = ncols(F);  //ncols better than size, since F[i] could be zero
+  // P is needed for default ORD
+  int i,j,k,l;
+  // st = "(a(1..(P)..1),a(1..(P+P^2)..1),dp)";
+  string st = "(a(" + string(1:P);
+  st = st + "),a(" + string(1:(P+P^2));
+  st = st + "),dp)";
+  // default values
+  string ORD = st;
+  int eng = 0;
+  if ( size(#)>0 )
+  {
+    if ( typeof(#[1]) == "string" )
+    {
+      ORD = string(#[1]);
+      // second arg
+      if (size(#)>1)
+      {
+        // exists 2nd arg
+        if ( typeof(#[2]) == "int" )
+        {
+          // the case: given ORD, given engine
+          eng = int(#[2]);
+        } 
+        else
+        {
+          eng = 0;
+        }
+      }
+      else
+      {
+        // no second arg
+        eng = 0;
+      }
+    }
+    else
+    {
+      if ( typeof(#[1]) == "int" )
+      {
+        // the case: default ORD, engine given
+        eng = int(#[1]);
+        // ORD = "(a(1..(P)..1),a(1..(P+P^2)..1),dp)";  //is already set
+      }
+      else
+      {
+        // incorr. 1st arg
+        ORD = string(st);
+      }
+    }
+  }
+  // size(#)=0, i.e. there is no elimination ordering and no engine given
+  // eng = 0; ORD = "(a(1..(P)..1),a(1..(P^2)..1),dp)";  //are already set
+  int ppl = printlevel-voice+2;
+  // returns a list with a ring and an ideal LD in it
+  // save, N, P and the indices are already defined 
+  int Nnew = 2*N+P+P^2;
+  list RL = ringlist(basering);
+  list L;
+  L[1] = RL[1];  //char
+  L[4] = RL[4];  //char, minpoly
+  // check whether vars have admissible names
+  list Name  = RL[2];
+  list RName;
+  // (i,j) <--> (i-1)*p+j
+  for (i=1; i<=P; i++)
+  {
+    RName[i] = "Dt("+string(i)+")";
+    for (j=1; j<=P; j++)
+    {
+      st = "s("+string(i)+")("+string(j)+")";
+      RName[P+(i-1)*P+j] = st;
+    }
+  }
+  for(i=1; i<=N; i++)
+  {
+    for(j=1; j<=size(RName); j++)
+    {
+      if (Name[i] == RName[j])
+      {
+        ERROR("Variable names should not include Dt(i), s(i)(j)");
+      }
+    }
+  }
+  // now, create the names for new vars
+  list DName;
+  for(i=1; i<=N; i++)
+  {
+    DName[i] = "D"+Name[i];  //concat
+  }
+  list NName = RName + Name + DName;
+  L[2] = NName;
+  // Name, Dname will be used further
+  kill NName;
+  //block ord (a(1..(P)..1),a(1..(P+P^2)..1),dp);
+  //export Nnew;
+  L[3] = ORDstr2list(ORD,Nnew);
+  // we are done with the list
+  def @R@ = ring(L);
+  setring @R@;
+  matrix @D[Nnew][Nnew];
+  // kronecker(i,j) equals (i==j)
+  // (i,j) <--> (i-1)*p+j
+  for (i=1; i<=P; i++)
+  {
+    for (j=1; j<=P; j++)
+    {
+      for (k=1; k<=P; k++)
+      {
+        //[sij,Dtk] = djk*Dti
+        @D[k,P+(i-1)*P+j] = (j==k)*Dt(i);
+        for (l=1; l<=P; l++)
+        {
+          if ( (i-k)*P < l-j )
+          {
+            //[sij,skl] = djk*sil - dil*skj
+            @D[P+(i-1)*P+j,P+(k-1)*P+l] = -(j==k)*s(i)(l) + (i==l)*s(k)(j);
+          }
+        }
+      }
+    }
+  }
+  for (i=1; i<=N; i++)
+  {
+    //[Dx,x]=1
+    @D[P+P^2+i,P+P^2+N+i] = 1;
+  }
+  def @R = nc_algebra(1,@D);
+  setring @R;
+  //@R@ will be used further
+  dbprint(ppl,"// -1-1- the ring @R(_Dt,_s,_x,_Dx) is ready");
+  dbprint(ppl-1, @R);
+  // create the ideal I
+  // (i,j) <--> (i-1)*p+j
+  ideal  F = imap(save,F);
+  ideal I;
+  for (i=1; i<=P; i++)
+  {
+    for (j=1; j<=P; j++)
+    {
+      I[(i-1)*P+j] = Dt(i)*F[j] + s(i)(j);
+    }
+  }
+  poly p,q;
+  for (i=1; i<=N; i++)
+  {
+    p=0;
+    for (j=1; j<=P; j++)
+    {
+      q = Dt(j);
+      q = q*diff(F[j],var(P+P^2+i));
+      p = p + q;
+    }
+    I = I, p + var(P+P^2+N+i);
+  }
+  // -------- the ideal I is ready ----------
+  dbprint(ppl,"// -1-2- starting the elimination of "+string(Dt(1..P))+" in @R");
+  dbprint(ppl-1, I);
+  ideal J = engine(I,eng);
+  ideal K = nselect(J,1..P);
+  kill I,J;
+  dbprint(ppl,"// -1-3- all Dt(i) are eliminated");
+  dbprint(ppl-1, K);  //K is without Dt(i)
+  // ----------- the ring @R2(_s,_x,_Dx) ------------
+  //come back to the ring save, recover L and remove all Dt(i)
+  //L[1],L[4] do not change
+  setring save;
+  list Lord, tmp;
+  // variables
+  tmp = L[2];
+  Lord = tmp[P+1..Nnew];
+  L[2] = Lord;
+  // ordering
+  // st = "(a(1..(P^2)..1),dp)";
+  st = "(a(" + string(1:P^2);
+  st = st + "),dp)";
+  tmp = ORDstr2list(st,Nnew-P);
+  L[3] = tmp;
+  def @R2@ = ring(L);
+  kill L;
+  // we are done with the list
+  setring @R2@;
+  matrix tmpM,LordM;
+  // non-commutative relations
+  intvec iv = P+1..Nnew;
+  tmpM = imap(@R@,@D);
+  kill @R@; 
+  LordM = submat(tmpM,iv,iv);
+  matrix @D2 = LordM;
+  def @R2 = nc_algebra(1,@D2);
+  setring @R2;
+  kill @R2@;
+  dbprint(ppl,"// -2-1- the ring @R(_s,_x,_Dx) is ready");
+  dbprint(ppl-1, @R2);
+  ideal K = imap(@R,K);
+  kill @R;
+  dbprint(ppl,"// -2-2- starting cosmetic Groebner basis computation");
+  dbprint(ppl-1, K);
+  K = engine(K,eng);
+  dbprint(ppl,"// -2-3- the cosmetic Groebner basis has been computed");
+  dbprint(ppl-1,K);
+  ideal LD = K;
+  attrib(LD,"isSB",1);
+  export LD;
+  return(@R2);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  ring R = 0,(x,y),Dp;
+  ideal F = x^3, y^5;
+  //ORD = "(a(1,1),a(1,1,1,1,1,1),dp)";
+  //eng = 0;
+  def A = SannfsVar(F);
+  setring A;
+  LD;
+}
+
+proc bfctVarAnn (ideal F, list #)
+"USAGE:  bfctVarAnn(F[,gid,eng]);  F an ideal, gid,eng optional ints
+RETURN:  list of ideal and intvec
+PURPOSE: computes the roots of the Bernstein-Sato polynomial and their
+@*       multiplicities for an affine algebraic variety defined by
+@*        F = F[1],..,F[r].
+ASSUME:  The basering is a commutative polynomial ring in char 0.
+BACKGROUND: In this proc, the annihilator of f^s in D[s] is computed and then a
+@*       system of linear equations is solved by linear reductions in order to
+@*       find the minimal polynomial of S = s(1)(1) + ... + s(P)(P), following
+@*       the ideas by Noro.
+NOTE:    In the output list, the ideal contains all the roots and the invec their
+@*        multiplicities.
+@*       If gid<>0, the ideal is used as given. Otherwise, and by default, a
+@*       heuristically better suited generating set is used.
+@*       If eng<>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.
+COMPUTATIONAL REMARK: The time of computation can be very different depending
+@*       on the chosen generators of F, although the result is always the same.
+EXAMPLE: example bfctVarAnn; shows examples
+"
+{
+  if (dmodvarAssumeViolation())
+  {
+    ERROR("Basering is inappropriate: characteristic>0 or qring present");
+  }
+  if (!isCommutative())
+  {
+    ERROR("Basering must be commutative");
+  }
+  int gid = 0; // default
+  int eng = 0; // default
+  if (size(#)>0)
+  {
+    if (typeof(#[1])=="int" || typeof(#[1])=="number")
+    {
+      gid = int(#[1]);
+    }
+    if (size(#)>1)
+    {
+      if (typeof(#[2])=="int" || typeof(#[2])=="number")
+      {
+	eng = int(#[2]);
+      }
+    }
+  }
+  def save = basering;
+  int ppl = printlevel - voice + 2;
+  printlevel = printlevel+1;
+  list L = smallGenCoDim(F,gid);
+  F = L[1];
+  int cd = L[2];
+  kill L;
+  def @R2 = SannfsVar(F,eng);
+  printlevel = printlevel-1;
+  setring @R2;
+  // we are in D[s] and LD is a std of SannfsVar(F)
+  ideal F = imap(save,F);
+  ideal GF = std(F);
+  ideal J = NF(LD,GF);
+  J = J, F;
+  dbprint(ppl,"// -3-1- starting Groebner basis of ann F^s + F ");
+  dbprint(ppl-1,J);
+  ideal K = engine(J,eng);
+  dbprint(ppl,"// -3-2- finished Groebner basis of ann F^s + F ");
+  dbprint(ppl-1,K);
+  poly S;
+  int i;
+  for (i=1; i<=size(F); i++)
+  {
+    S = S + s(i)(i);
+  }
+  dbprint(ppl,"// -4-1- computing the minimal polynomial of S");
+  //dbprint(ppl-1,"S = "+string(S));
+  module M = pIntersect(S,K);
+  dbprint(ppl,"// -4-2- the minimal polynomial has been computed");
+  //dbprint(ppl-1,M);
+  ring @R3 = 0,s,dp;
+  dbprint(ppl,"// -5-1- the ring @R3(s) is ready");
+  dbprint(ppl-1,@R3);
+  ideal M = imap(@R2,M);
+  //kill @R2;
+  poly p;
+  for (i=1; i<=size(M); i++)
+  {
+    p = p + M[i]*s^(i-1);
+  }
+  dbprint(ppl,"// -5-2- factorization of the minimal polynomial");
+  list P = factorize(p);          //with constants and multiplicities
+  dbprint(ppl-1,P);               //the Bernstein polynomial is monic,
+  ideal bs; intvec m;             //so we are not interested in constants
+  for (i=2; i<=ncols(P[1]); i++)   //and that is why we delete P[1][1] and P[2][1]
+  {
+    bs[i-1] = P[1][i];
+    m[i-1]  = P[2][i];
+  }
+  // convert factors to a list of their roots and multiplicities
+  bs =  normalize(bs);
+  bs = -subst(bs,s,0);
+  setring save;
+//   ideal GF = imap(@R2,GF);
+//   attrib(GF,"isSB",1);
+  kill @R2;
+  dbprint(ppl,"// -5-3- codimension of the variety");
+//   int cd = coDim(GF);
+  dbprint(ppl-1,cd);
+  ideal bs = imap(@R3,bs);
+  dbprint(ppl,"// -5-4- shifting BS(s)=minpoly(s-codim+1)");
+  for (i=1; i<=ncols(bs); i++)
+  {
+    bs[i] = bs[i] + cd - 1;
+  }
+  kill @R3;
+  list BS = bs,m;
+  return(BS);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  ring R = 0,(x,y,z),Dp;
+  ideal F = x^2+y^3, z;
+  bfctVarAnn(F);
+}
+
+proc makeIF (ideal F, list #)
+"USAGE:  makeIF(F [,ORD]);  F an ideal, ORD an optional string
+RETURN:  ring
+PURPOSE: create the ideal by Malgrange associated with F = F[1],..,F[P].
+NOTE:    activate this ring with the @code{setring} command. In this ring,
+@*       - the ideal IF is the ideal by Malgrange corresponding to F.
+@*       The value of ORD must be an arbitrary ordering in K<_t,_x,_Dt,_Dx>
+@*       written in the string form. By default ORD = 'dp'.
+@*       If printlevel=1, progress debug messages will be printed,
+@*       if printlevel>=2, all the debug messages will be printed.
+EXAMPLE: example makeIF; shows examples
+"
+{
+  string ORD = "dp";
+  if ( size(#)>0 )
+  {
+    if ( typeof(#[1]) == "string" )
+    {
+      ORD = string(#[1]);
+    }
+  }
+  int ppl = printlevel-voice+2;
+  def save = basering;
+  int N = nvars(save);
+  int P = ncols(F);
+  int Nnew = 2*P+2*N;
+  int i,j;
+  string st;
+  list RL = ringlist(save);
+  list L,Lord;
+  list tmp;
+  intvec iv;
+  L[1] = RL[1];
+  L[4] = RL[4];
+  //check whether vars have admissible names
+  list Name = RL[2];
+  list TName, DTName;
+  for (i=1; i<=P; i++)
+  {
+    TName[i] = "t("+string(i)+")";
+    DTName[i] = "Dt("+string(i)+")";
+  }
+  for (i=1; i<=N; i++)
+  {
+    for (j=1; j<=P; j++)
+    {
+      if (Name[i] == TName[j])
+      {
+        ERROR("Variable names should not include t(i)");
+      }
+    }
+  }
+  //now, create the names for new vars
+  list DName;
+  for (i=1; i<=N; i++)
+  {
+    DName[i] = "D"+Name[i];  //concat
+  }
+  list NName = TName + Name + DTName + DName;
+  L[2]   = NName;
+  // Name, Dname will be used further
+  kill NName, TName, Name, DTName, DName;
+  // ORD already set, default ord dp; 
+  L[3] = ORDstr2list(ORD,Nnew);
+  // we are done with the list
+  def @R@ = ring(L);
+  setring @R@;
+  matrix @D[Nnew][Nnew];
+  for (i=1; i<=N+P; i++)
+  {
+    @D[i,i+N+P]=1;
+  }
+  def @R = nc_algebra(1,@D);
+  setring @R;
+  kill @R@;
+  //dbprint(ppl,"// -1-1- the ring @R(_t,_x,_Dt,_Dx) is ready");
+  //dbprint(ppl-1, @R);
+  // create the ideal I
+  ideal  F = imap(save,F);
+  ideal I;
+  for (j=1; j<=P; j++)
+  {
+    I[j] = t(j) - F[j];
+  }
+  poly p,q;
+  for (i=1; i<=N; i++)
+  {
+    p=0;
+    for (j=1; j<=P; j++)
+    {
+      q = Dt(j);
+      q = diff(F[j],var(P+i))*q;
+      p = p + q;
+    }
+    I = I, p + var(2*P+N+i);
+  }
+  // -------- the ideal I is ready ----------
+  ideal IF = I;
+  export IF;
+  return(@R);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  ring R = 0,(x,y,z),Dp;
+  ideal I = x^2+y^3, z;
+  def W = makeIF(I);
+  setring W;
+  IF;
+}
+
+proc bfctVarIn (ideal I, list #)
+"USAGE:  bfctVarIn(I [,a,b,c]);  I an ideal, a,b,c optional ints
+RETURN:  list of ideal and intvec
+PURPOSE: computes the roots of the Bernstein-Sato polynomial and their
+@*       multiplicities for an affine algebraic variety defined by I.
+ASSUME:  The basering is commutative and of characteristic 0.
+@*       Varnames of the basering do not include t(1),...,t(r) and
+@*       Dt(1),...,Dt(r), where r is the number of entries of the input ideal.
+BACKGROUND:  In this proc, the initial ideal of the multivariate Malgrange ideal 
+@*       defined by I is computed and then a system of linear equations is solved
+@*       by linear reductions following the ideas by Noro.
+NOTE:    In the output list, say L,
+@*       - L[1] of type ideal contains all the rational roots of a b-function,
+@*       - L[2] of type intvec contains the multiplicities of above roots,
+@*       - optional L[3] of type string is the part of b-function without
+@*        rational roots.
+@*       Note, that a b-function of degree 0 is encoded via L[1][1]=0, L[2]=0 and
+@*       L[3] is 1 (for nonzero constant) or 0 (for zero b-function).
+@*       If a<>0, the ideal is used as given. Otherwise, and by default, a
+@*       heuristically better suited generating set is used to reduce computation
+@*       time.
+@*       If b<>0, @code{std} is used for GB computations in characteristic 0,
+@*       otherwise, and by default, @code{slimgb} is used. 
+@*       If c<>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 bfctVarIn; shows examples
+"
+{
+  if (dmodvarAssumeViolation())
+  {
+    ERROR("Basering is inappropriate: characteristic>0 or qring present");
+  }
+  if (!isCommutative())
+  {
+    ERROR("Basering must be commutative");
+  }
+  int ppl = printlevel - voice + 2;
+  int idealasgiven  = 0; // default
+  int whicheng      = 0; // default
+  int whichord      = 0; // default
+  if (size(#)>0)
+  {
+    if (typeof(#[1])=="int" || typeof(#[1])=="number")
+    {
+      idealasgiven = int(#[1]);
+    }
+    if (size(#)>1)
+    {
+      if (typeof(#[2])=="int" || typeof(#[2])=="number")
+      {
+	whicheng = int(#[2]);
+      }
+      if (size(#)>2)
+      {
+	if (typeof(#[3])=="int" || typeof(#[3])=="number")
+	{
+	  whichord = int(#[3]);
+	}
+      }
+    }
+  }
+  def save = basering;
+  int i;
+  int n = nvars(basering);
+  // step 0: get small generating set
+  I = simplify(I,2);
+  list L = smallGenCoDim(I,idealasgiven);
+  I = L[1];
+  int c = L[2];
+  kill L;
+  // step 1: setting up the multivariate Malgrange ideal
+  int r = size(I);
+  def D = makeIF(I);
+  setring D;
+  dbprint(ppl-1,"// Computing in " + string(n+r) + "-th Weyl algebra:", D);
+  dbprint(ppl-1,"// The Malgrange ideal: ", IF);
+  // step 2: compute the b-function of the Malgrange ideal w.r.t. approriate weights
+  intvec w = 1:r;
+  w[r+n] = 0;
+  dbprint(ppl,"// Computing the b-function of the Malgrange ideal...");
+  list L = bfctIdeal(IF,w,whicheng,whichord);
+  dbprint(ppl,"// ... done.");
+  dbprint(ppl-1,"// The b-function: ",L);
+  // step 3: shift the result
+  ring S = 0,s,dp;
+  list L = imap(D,L);
+  kill D;
+  if (size(L)==2)
+  {
+    ideal B = L[1];
+    for (i=1; i<=ncols(B); i++)
+    {
+      B[i] = -B[i]+c-r-1;
+    }
+    L[1] = B;
+  }
+  else // should never get here: BS poly has non-rational roots
+  {
+    string str = L[3];
+    L = delete(L,3);
+    str = "poly @b = (" + str + ")*(" + string(fl2poly(L,"s")) + ");";
+    execute(str);
+    poly b = subst(@b,s,-s+c-r-1);
+    L = bFactor(b);
+  }
+  setring save;
+  list L = imap(S,L);
+  return(L);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  ring R = 0,(x,y,z),dp;
+  ideal F = x^2+y^3, z;
+  list L = bfctVarIn(F);
+  L;
+}
+
+static proc smallGenCoDim (ideal I, int Iasgiven)
+{
+  // call from K[x]
+  // returns list L
+  // L[1]=I or L[1]=smaller generating set of I
+  // L[2]=codimension(I)
+  int ppl = printlevel - voice + 3;
+  int n = nvars(basering);
+  int c;
+  if (attrib(I,"isSB") == 1)
+  {
+    c = n - dim(I);
+    if (!Iasgiven)
+    {
+      list L = mstd(I);
+    }
+  }
+  else
+  {
+    def save = basering;
+    list RL = ringlist(save);
+    list @ord;
+    @ord[1] = list("dp", intvec(1:n));
+    @ord[2] = list("C", intvec(0));
+    RL[3] = @ord;
+    kill @ord;
+    if (size(RL)>4) // commutative G-algebra present
+    {
+      RL = RL[1..4];
+    }
+    def R = ring(RL);
+    kill RL;
+    setring R;
+    ideal I = imap(save,I);
+    if (!Iasgiven)
+    {
+      list L = mstd(I);
+      c = n - dim(L[1]);
+      setring save;
+      list L = imap(R,L);
+    }
+    else
+    {
+      I = std(I);
+      c = n - dim(I);
+      setring save;
+    }
+    kill R;
+  }
+  if (!Iasgiven)
+  {
+    if (size(L[2]) < size(I))
+    {
+      I = L[2];
+      dbprint(ppl,"// Found smaller generating set of the given variety: ", I);
+    }
+    else
+    {
+      dbprint(ppl,"// Have not found smaller generating set of the given variety.");
+    }
+  }
+  dbprint(ppl-1,"// The codim of the given variety is " + string(c) + ".");
+  if (!defined(L))
+  {
+    list L;
+  }
+  L[1] = I;
+  L[2] = c;
+  return(L);
+}
+
+
+// Some more examples
+
+static proc TXcups()
+{
+  "EXAMPLE:"; echo = 2;
+  //TX tangent space of X=V(x^2+y^3)
+  ring R = 0,(x0,x1,y0,y1),Dp;
+  ideal F = x0^2+y0^3, 2*x0*x1+3*y0^2*y1;
+  printlevel = 0;
+  //ORD = "(a(1,1),a(1,1,1,1,1,1),dp)";
+  //eng = 0;
+  def A = SannfsVar(F);
+  setring A;
+  LD;
+}
+
+static proc ex47()
+{
+  ring r7 = 0,(x0,x1,y0,y1),dp;
+  ideal I = x0^2+y0^3, 2*x0*x1+3*y0^2*y1;
+  bfctVarIn(I);
+  // second ex - too big
+  ideal J = x0^4+y0^5, 4*x0^3*x1+5*y0^4*y1;
+  bfctVarIn(J);
+}
+
+static proc ex48()
+{
+  ring r8 = 0,(x1,x2,x3),dp;
+  ideal I = x1^3-x2*x3, x2^2-x1*x3, x3^2-x1^2*x2;
+  bfctVarIn(I);
+}
+
+static proc ex49 ()
+{
+  ring r9 = 0,(z1,z2,z3,z4),dp;
+  ideal I = z3^2-z2*z4, z2^2*z3-z1*z4, z2^3-z1*z3;
+  bfctVarIn(I);
+}
+
+static proc ex410()
+{
+  LIB "toric.lib";
+  ring r = 0,(z(1..7)),dp;
+  intmat A[3][7];
+  A = 6,4,2,0,3,1,0,0,1,2,3,0,1,0,0,0,0,0,1,1,2;
+  ideal I = toric_ideal(A,"pt");
+  I = std(I);
+//ideal I = z(6)^2-z(3)*z(7), z(5)*z(6)-z(2)*z(7), z(5)^2-z(1)*z(7),
+//  z(4)*z(5)-z(3)*z(6), z(3)*z(5)-z(2)*z(6), z(2)*z(5)-z(1)*z(6), 
+//  z(3)^2-z(2)*z(4), z(2)*z(3)-z(1)*z(4), z(2)^2-z(1)*z(3);
+  bfctVarIn(I,1); // no result yet
+}
