source: git/Singular/LIB/bfct.lib @ a470eb

//////////////////////////////////////////////////////////////////////////////
version="$Id: bfct.lib,v 1.7 2008-12-09 16:50:21 levandov Exp $";
category="Noncommutative";
info="
LIBRARY: bfct.lib     M. Noro's Algorithm for 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.

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

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

SEE ALSO: dmod_lib, dmodapp_lib, gmssing_lib
";


LIB "qhmoduli.lib"; // for Max
LIB "dmod.lib"; // for SannfsBFCT etc
LIB "dmodapp.lib";  // for initialideal etc


proc testbfctlib ()
{
  // tests all procs for consistency
  "AUXILIARY PROCEDURES:";
  example ispositive;
  example isin;
  example scalarprod;
  example vec2poly;
  "MAIN PROCEDURES:";
  example bfct;
  example bfctsyz;
  example bfctann;
  example bfctonestep;
  example bfctideal;
  example pintersect;
  example pintersectsyz;
  example linreduce;
  example ncsolve;
}


//--------------- 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
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;
  LIB "bfct.lib";
  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 ispositive (intvec v)
"USAGE:  ispositive(v);  v an intvec
RETURN:  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
"
{
  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;
  ispositive(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
PURPOSE: check whether the second argument is a member of a list
EXAMPLE: example isin; shows an example
"
{
  int j;
  for (j=1; j<=size(l); j++)
  {
    if (l[j]==i)
    {
      return(j);
      break;
    }
  }
  return(0);
}
example
{
  "EXAMPLE:"; echo = 2;
  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
"
{
  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 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
"
{
  int ppl = printlevel - voice + 2;
  int remembercoeffs = 0; // default
  int redlm          = 0; // default
  if (size(#)>0)
  {
    if (typeof(#[1])=="int" || typeof(#[1])=="number")
    {
      remembercoeffs = #[1];
    }
    if (size(#)>1)
    {
      if (typeof(#[2])=="int" || typeof(#[2])=="number")
      {
        redlm = #[2];
      }
    }
  }
  int i,j,k;
  int sI = ncols(I);
  ideal lmI,lcI;
  for (i=1; i<=sI; i++)
  {
    lmI[i] = leadmonom(I[i]);
    lcI[i] = leadcoef(I[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++)
        {
          if (monomf[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);
            }
            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);
          }
        }
      }
    }
  }
  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);
}

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
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,
@*       if printlevel>=2, all the debug messages will be printed.
EXAMPLE: example ncsolve; 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;
    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, "ncsolve: 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");
      // 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;
  ncsolve(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
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,
@*       if printlevel>=2, all the debug messages will be printed.
EXAMPLE: example pintersect; shows examples
"
{
  // assume I is given in Groebner basis
  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))
    {
      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:  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.
@*       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,
@*       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;
  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,"ncsolve starts with: "+string(matrix(NI)));
    if (solveincharp<>0) // modular method
    {
      setring Rp;
      NI[i+1] = phi(newNF);
      v = ncsolve(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);
        if (v==0)
        {
          break;
        }
      }
      else // no modular solution
      {
        setring save;
        v = 0;
      }
    }
    else // non-modular method
    {
      v = ncsolve(NI,whichengine);
    }
    matrix MM[1][nrows(v)] = matrix(v);
    dbprint(ppl,"ncsolve ready  with: "+string(MM));
    kill MM;
    //  "ncsolve 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,"ncsolve: bad solution!");
      }
      else
      {
        dbprint(ppl,"ncsolve: 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:  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.
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 (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 (isin(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 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,
@*       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.
@*       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 && ispositive(#[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 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,
@*       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 printlevel>=2, all the debug messages will be printed.
EXAMPLE: example bfct; 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 && ispositive(#[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 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,
@*       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 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 = initialideal(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 = 3*x^2*Dy+2*y*Dx,2*x*Dx+3*y*Dy+6;
  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,
@*       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 printlevel>=2, all the debug messages will be printed.
EXAMPLE: example bfctonestep; shows examples
"
{
  int ppl = printlevel - voice +2;
  def save = basering;
  int n = nvars(save);
  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]);
      }
    }
  }
  def DDh = initialidealengine("bfctonestep", whichengine, methodord, f);
  setring DDh;
  dbprint(ppl, "the initial ideal:", string(matrix(inF)));
  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;
  setring save;
  ideal J; J[2] = var(1);
  map @m = DDh,J;
  ideal inF = @m(inF);
  poly p = inF[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;
  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,
@*       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 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);
}