//////////////////////////////////////////////////////////////////////////////
version="$Id$";
category="Noncommutative";
info="
LIBRARY: dmod.lib     Algorithms for algebraic D-modules
AUTHORS: Viktor Levandovskyy,     levandov@math.rwth-aachen.de
@*             Jorge Martin Morales,    jorge@unizar.es

OVERVIEW:
Let K be a field of characteristic 0. Given a polynomial ring
@*      R = K[x_1,...,x_n] and a polynomial F in R,
@*      one is interested in the R[1/F]-module of rank one, generated by
@*      the symbol F^s for a symbolic discrete variable s.
@* In fact, the module R[1/F]*F^s has a structure of a D(R)[s]-module, 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> and
@* D(R)[s] = D(R) tensored with K[s] over K.
@* Constructively, one needs to find a left ideal I = I(F^s) in D(R), such
@* that K[x_1,...,x_n,1/F]*F^s is isomorphic to D(R)/I as a D(R)-module.
@* We often write just D for D(R) and D[s] for D(R)[s].
@* One is interested in the following data:
@* - Ann F^s = I = I(F^s) in D(R)[s], denoted by LD in the output
@* - global Bernstein polynomial in K[s], 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
@* - Ann F^s0 = I(F^s0) in D(R), denoted by LD0 in the output
@*   (LD0 is a holonomic ideal in D(R))
@* - Ann^(1) F^s in D(R)[s], denoted by LD1 (logarithmic derivations)
@* - an operator in D(R)[s], denoted by PS, such that the functional equality
@*     PS*F^(s+1) = bs*F^s holds in K[x_1,...,x_n,1/F]*F^s.

REFERENCES:
@* We provide the following implementations of algorithms:
@* (OT) the classical Ann F^s algorithm from Oaku and Takayama (Journal of
@* Pure and Applied Math., 1999),
@* (LOT) Levandovskyy's modification of the Oaku-Takayama algorithm (ISSAC 2007)
@* (BM) the Ann F^s algorithm by Briancon and Maisonobe (Remarques sur
@*        l'ideal de Bernstein associe a des polynomes, preprint, 2002)
@* (LM08) V. Levandovskyy and J. Martin-Morales, ISSAC 2008
@* (C) Countinho, A Primer of Algebraic D-Modules,
@* (SST) Saito, Sturmfels, Takayama 'Groebner Deformations of Hypergeometric
@*         Differential Equations', Springer, 2000


Guide:
@* - Ann F^s = I(F^s) = LD in D(R)[s] can be computed by Sannfs [BM, OT, LOT]
@* - Ann^(1) F^s in D(R)[s] can be computed by Sannfslog
@* - global Bernstein polynomial bs in K[s] can be computed by bernsteinBM
@* - Ann F^s0 = I(F^s0) = LD0 in D(R) can be computed by annfs0, annfs, annfsBM,
@*    annfsOT, annfsLOT, annfs2, annfsRB etc.
@* - all the relevant data to F^s (LD, LD0, bs, PS) are computed by operatorBM
@* - operator PS can be computed via operatorModulo or operatorBM
@*
@* - annihilator of F^{s1} for a number s1 is computed with annfspecial
@* - annihilator of F_1^s_1 * ... * F_p^s_p is computed with annfsBMI
@* - computing the multiplicity of a rational number r in the Bernstein poly
@*   of a given ideal goes with checkRoot
@* - check, whether a given univariate polynomial divides the Bernstein poly
@*   goes with checkFactor


PROCEDURES:

annfs(F[,S,eng]);       compute Ann F^s0 in D and Bernstein polynomial for a poly F
annfspecial(I, F, m, n);  compute Ann F^n from Ann F^s for a polynomial F and a number n
Sannfs(F[,S,eng]);      compute Ann F^s in D[s] for a polynomial F
Sannfslog(F[,eng]);     compute Ann^(1) F^s in D[s] for a polynomial F
bernsteinBM(F[,eng]);   compute global Bernstein polynomial for a polynomial F (algorithm of Briancon-Maisonobe)
bernsteinLift(I,F [,eng]);  compute a possible multiple of Bernstein polynomial via lift-like procedure
operatorBM(F[,eng]);    compute Ann F^s, Ann F^s0, BS and PS for a polynomial F (algorithm of Briancon-Maisonobe)
operatorModulo(F, I, b); compute PS via the modulo approach
annfsParamBM(F[,eng]);  compute the generic Ann F^s (algorithm by Briancon and Maisonobe) and exceptional parametric constellations for a polynomial F with parametric coefficients
annfsBMI(F[,eng]);      compute Ann F^s and Bernstein ideal for a polynomial F=f1*..*fP (multivariate algorithm of Briancon-Maisonobe)
checkRoot(F,a[,S,eng]); check if a given rational is a root of the global Bernstein polynomial of F and compute its multiplicity
SannfsBFCT(F[,eng]);      compute Ann F^s in D[s] for a polynomial F (algorithm of Briancon-Maisonobe, other output ordering)
annfs0(I,F [,eng]);          compute Ann F^s0 in D and Bernstein polynomial from the known Ann F^s in D[s]
annfs2(I,F [,eng]);           compute Ann F^s0 in D and Bernstein polynomial from the known Ann F^s in D[s] by using a trick of Noro
annfsRB(I,F [,eng]);          compute Ann F^s0 in D and Bernstein polynomial from the known Ann F^s in D[s] by using Jacobian ideal
checkFactor(I,F,q[,eng]); check whether a polynomial q in K[s] is a factor of the global Bernstein polynomial of F from the known Ann F^s in D[s]


arrange(p);           create a poly, describing a full hyperplane arrangement
reiffen(p,q);         create a poly, describing a Reiffen curve
isHolonomic(M);   check whether a module is holonomic
convloc(L);         replace global orderings with local in the ringlist L
minIntRoot(P,fact); minimal integer root among the roots of a maximal ideal P
varNum(s);          the number of the variable with the name s
isRational(n);  check whether n is a rational number

SEE ALSO: bfun_lib, dmodapp_lib, dmodvar_lib, gmssing_lib

KEYWORDS: D-module; D-module structure; left annihilator ideal; Bernstein-Sato polynomial; global Bernstein-Sato polynomial;
Weyl algebra; Bernstein operator; logarithmic annihilator ideal; parametric annihilator; root of Bernstein-Sato polynomial;
hyperplane arrangement; Oaku-Takayama algorithm; Briancon-Maisonobe algorithm; LOT algorithm
";

// reworked by JM+VL on 9.3.2010: Sannfslog
// added by VL on 2.3.2010: bernsteinLift
// ****** commented out for better readability by VL on 2.3.2010
// annfsBM(F[,eng]);          compute Ann F^s0 in D and Bernstein polynomial for a polynomial F (algorithm of Briancon-Maisonobe)
// annfsLOT(F[,eng]);         compute Ann F^s0 in D and Bernstein polynomial for a polynomial F (Levandovskyy modification of the Oaku-Takayama algorithm)
// annfsOT(F[,eng]);          compute Ann F^s0 in D and Bernstein polynomial for a polynomial F (algorithm of Oaku-Takayama)
// SannfsBM(F[,eng]);         compute Ann F^s in D[s] for a polynomial F (algorithm of Briancon-Maisonobe)
// SannfsLOT(F[,eng]);        compute Ann F^s in D[s] for a polynomial F (Levandovskyy modification of the Oaku-Takayama algorithm)
// SannfsOT(F[,eng]);         compute Ann F^s in D[s] for a polynomial F (algorithm of Oaku-Takayama)
// checkRoot1(I,F,a[,eng]);   check whether a rational is a root of the global Bernstein polynomial of F from the known Ann F^s in D[s]
// checkRoot2(I,F,a[,eng]);   check whether a rational is a root of the global Bernstein polynomial of F and compute its multiplicity from the known Ann F^s in D[s]

LIB "matrix.lib"; // for submat
LIB "nctools.lib"; // makeModElimRing etc.
LIB "elim.lib"; // for nselect
LIB "qhmoduli.lib"; // for Max
LIB "gkdim.lib";
LIB "gmssing.lib";
LIB "control.lib";  // for genericity
LIB "dmodapp.lib"; // for e.g. engine
LIB "poly.lib";

proc testdmodlib()
{
  /* tests all procs for consistency */
  /* adding the new proc, add it here */

  "MAIN PROCEDURES:";
  example annfs;
  example Sannfs;
  example Sannfslog;
  example bernsteinBM;
  example operatorBM;
  example annfspecial;
  example annfsParamBM;
  example annfsBMI;
  example checkRoot;
  example annfs0;
  example annfs2;
  example annfsRB;
  example annfs2;
  example operatorModulo;
  "SECONDARY D-MOD PROCEDURES:";
  example annfsBM;
  example annfsLOT;
  example annfsOT;
  example SannfsBM;
  example SannfsLOT;
  example SannfsOT;
  example SannfsBFCT;
  example checkRoot1;
  example checkRoot2;
  example checkFactor;
}



// new top-level procedures
proc annfs(poly F, list #)
"USAGE:  annfs(f [,S,eng]);  f a poly, S a string, eng an optional int
RETURN:  ring
PURPOSE: compute the D-module structure of basering[1/f]*f^s with the algorithm
@*  given in S and with the Groebner basis engine given in ''eng''
NOTE:  activate the output ring with the @code{setring} command.
@*    String S; S can be one of the following:
@*    'bm' (default) - for the algorithm of Briancon and Maisonobe,
@*    'ot'  - for the algorithm of Oaku and Takayama,
@*    'lot' - for the Levandovskyy's modification of the algorithm of OT.
@*  If eng <>0, @code{std} is used for Groebner basis computations,
@*  otherwise and by default @code{slimgb} is used.
@*  In the output ring:
@*  - the ideal LD (which is a Groebner basis) is the needed D-module structure,
@*  - the list  BS contains roots and multiplicities of a BS-polynomial of f.
DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
@*          if @code{printlevel}>=2, all the debug messages will be printed.
EXAMPLE: example annfs; shows examples
"
{
  int eng = 0;
  int chs = 0; // choice
  string algo = "bm";
  string st;
  if ( size(#)>0 )
  {
   if ( typeof(#[1]) == "string" )
   {
     st = string(#[1]);
     if ( (st == "BM") || (st == "Bm") || (st == "bM") ||(st == "bm"))
     {
       algo = "bm";
       chs  = 1;
     }
     if ( (st == "OT") || (st == "Ot") || (st == "oT") || (st == "ot"))
     {
       algo = "ot";
       chs  = 1;
     }
     if ( (st == "LOT") || (st == "lOT") || (st == "Lot") || (st == "lot"))
     {
       algo = "lot";
       chs  = 1;
     }
     if (chs != 1)
     {
       // incorrect value of S
       print("Incorrect algorithm given, proceed with the default BM");
       algo = "bm";
     }
     // second arg
     if (size(#)>1)
     {
       // exists 2nd arg
       if ( typeof(#[2]) == "int" )
       {
         // the case: given alg, given engine
         eng = int(#[2]);
       }
       else
       {
         eng = 0;
       }
     }
     else
     {
       // no second arg
       eng = 0;
     }
   }
   else
   {
     if ( typeof(#[1]) == "int" )
     {
     // the case: default alg, engine
     eng = int(#[1]);
     //      algo = "bm";  //is already set
     }
     else
     {
       // incorr. 1st arg
       algo = "bm";
     }
   }
  }

   // size(#)=0, i.e. there is no algorithm and no engine given
   //  eng = 0; algo = "bm";  //are already set
   // int ppl = printlevel-voice+2;

  int old_printlevel = printlevel;
  printlevel=printlevel+1;
  def save = basering;
  if ( algo=="ot")
  {
    def @A = annfsOT(F,eng);
  }
  else
  {
    if ( algo=="lot")
    {
      def @A = annfsLOT(F,eng);
    }
    else
    {
      // bm = default
      def @A = annfsBM(F,eng);
    }
  }
  printlevel = old_printlevel;
  return(@A);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),Dp;
  poly F = z*x^2+y^3;
  def A  = annfs(F); // here, the default BM algorithm will be used
  setring A; // the Weyl algebra in (x,y,z,Dx,Dy,Dz)
  LD; //the annihilator of F^{-1} over A
  BS; // roots with multiplicities of BS polynomial
}


proc Sannfs(poly F, list #)
"USAGE:  Sannfs(f [,S,eng]);  f a poly, S a string, eng an optional int
RETURN:  ring
PURPOSE: compute the D-module structure of basering[f^s] with the algorithm
@*  given in S and with the Groebner basis engine given in eng
NOTE:    activate the output ring with the @code{setring} command.
@*       The value of a string S can be
@*       'bm' (default) - for the algorithm of  Briancon and Maisonobe,
@*       'lot' - for the Levandovskyy's modification of the algorithm of OT,
@*       'ot'  - for the algorithm of Oaku and Takayama.
@*       If eng <>0, @code{std} is used for Groebner basis computations,
@*       otherwise, and by default @code{slimgb} is used.
@*       In the output ring:
@*       - the ideal LD is the needed D-module structure.
DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
@*          if @code{printlevel}>=2, all the debug messages will be printed.
EXAMPLE: example Sannfs; shows examples
"
{
  int eng = 0;
  int chs = 0; // choice
  string algo = "bm";
  string st;
  if ( size(#)>0 )
  {
   if ( typeof(#[1]) == "string" )
   {
     st = string(#[1]);
     if ( (st == "BM") || (st == "Bm") || (st == "bM") ||(st == "bm"))
     {
       algo = "bm";
       chs  = 1;
     }
     if ( (st == "OT") || (st == "Ot") || (st == "oT") || (st == "ot"))
     {
       algo = "ot";
       chs  = 1;
     }
     if ( (st == "LOT") || (st == "lOT") || (st == "Lot") || (st == "lot"))
     {
       algo = "lot";
       chs  = 1;
     }
     if (chs != 1)
     {
       // incorrect value of S
       print("Incorrect algorithm given, proceed with the default BM");
       algo = "bm";
     }
     // second arg
     if (size(#)>1)
     {
       // exists 2nd arg
       if ( typeof(#[2]) == "int" )
       {
         // the case: given alg, given engine
         eng = int(#[2]);
       }
       else
       {
         eng = 0;
       }
     }
     else
     {
       // no second arg
       eng = 0;
     }
   }
   else
   {
     if ( typeof(#[1]) == "int" )
     {
     // the case: default alg, engine
     eng = int(#[1]);
     //      algo = "bm";  //is already set
     }
     else
     {
       // incorr. 1st arg
       algo = "bm";
     }
   }
  }
  // size(#)=0, i.e. there is no algorithm and no engine given
  //  eng = 0; algo = "bm";  //are already set
  // int ppl = printlevel-voice+2;
  printlevel=printlevel+1;
  def save = basering;
  if ( algo=="ot")
  {
    def @A = SannfsOT(F,eng);
  }
  else
  {
    if ( algo=="lot")
    {
      def @A = SannfsLOT(F,eng);
    }
    else
    {
      // bm = default
      def @A = SannfsBM(F,eng);
    }
  }
  printlevel=printlevel-1;
  return(@A);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),Dp;
  poly F = x^3+y^3+z^3;
  printlevel = 0;
  def A  = Sannfs(F); // here, the default BM algorithm will be used
  setring A;
  LD;
}

proc Sannfslog (poly F, list #)
"USAGE:  Sannfslog(f [,eng]);  f a poly, eng an optional int
RETURN:  ring
PURPOSE: compute the D-module structure of basering[1/f]*f^s
NOTE:    activate the output ring with the @code{setring} command.
@*   In the output ring D[s], the ideal LD1 is generated by the elements
@*   in Ann F^s in D[s], coming from logarithmic derivations.
@*       If eng <>0, @code{std} is used for Groebner basis computations,
@*       otherwise, and by default @code{slimgb} is used.
DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
@*          if @code{printlevel}>=2, all the debug messages will be printed.
EXAMPLE: example Sannfslog; shows examples
"
{
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  int ppl = printlevel-voice+2;
  def save = basering;
  int N = nvars(basering);
  int Nnew = 2*N+1;
  int i;
  string s;
  list RL = ringlist(basering);
  list L, Lord;
  list tmp;
  intvec iv;
  L[1] = RL[1]; // char
  L[4] = RL[4]; // char, minpoly
  // check whether vars have admissible names
  list Name = RL[2];
  for (i=1; i<=N; i++)
  {
    if (Name[i] == "s")
    {
      ERROR("Variable names should not include s");
    }
  }
  // the ideal I
  ideal I = -F, jacob(F);
  dbprint(ppl,"// -1-1- starting the computation of syz(-F,_Dx(F))");
  dbprint(ppl-1, I);
  matrix M = syz(I);
  M = transpose(M);  // it is more usefull working with columns
  dbprint(ppl,"// -1-2- the module syz(-F,_Dx(F)) has been computed");
  dbprint(ppl-1, M);
  // ------------ the ring @R ------------
  // _x, _Dx, s;  elim.ord for _x,_Dx.
  // now, create the names for new vars
  list DName;
  for (i=1; i<=N; i++)
  {
    DName[i] = "D"+Name[i]; // concat
  }
  tmp[1] = "s";
  list NName;
  NName = Name + DName + tmp;
  L[2] = NName;
  tmp = 0;
  // block ord (dp(N),dp);
  s = "iv=";
  for (i=1; i<=Nnew-1; i++)
  {
    s = s+"1,";
  }
  s[size(s)]=";";
  execute(s);
  tmp[1] = "dp";  // string
  tmp[2] = iv;   // intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1] = "dp";  // string
  s[size(s)] = ",";
  s = s+"1;";
  execute(s);
  kill s;
  kill NName;
  tmp[2]      = iv;
  Lord[2]     = tmp;
  tmp[1]      = "C";  iv  = 0;  tmp[2]=iv;
  Lord[3]     = tmp;  tmp = 0;
  L[3]        = Lord;
  // we are done with the list. Now add a Plural part
  def @R@ = ring(L);
  setring @R@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R = nc_algebra(1,@D);
  setring @R;
  kill @R@;
  dbprint(ppl,"// -2-1- the ring @R(_x,_Dx,s) is ready");
  dbprint(ppl-1, @R);
  matrix M = imap(save,M);
  // now, create the vector [-s,_Dx]
  vector v = [-s];  // now s is a variable
  for (i=1; i<=N; i++)
  {
    v = v + var(i+N)*gen(i+1);
  }
  ideal J = ideal(M*v);
  // make leadcoeffs positive
  for (i=1; i<= ncols(J); i++)
  {
    if ( leadcoef(J[i])<0 )
    {
      J[i] = -J[i];
    }
  }
  ideal LD1 = J;
  kill J;
  export LD1;
  return(@R);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),Dp;
  poly F = x4+y5+x*y4;
  printlevel = 0;
  def A  = Sannfslog(F);
  setring A;
  LD1;
}

// JM+VL: output ring restructured into "normal"

// proc Sannfslog (poly F, list #)
// "USAGE:  Sannfslog(f [,eng]);  f a poly, eng an optional int
// RETURN:  ring
// PURPOSE: compute the D-module structure of basering[1/f]*f^s
// NOTE:    activate the output ring with the @code{setring} command.
// @*   In the output ring D[s], the ideal LD1 is generated by the elements
// @*   in Ann F^s in D[s], coming from logarithmic derivations.
// @*       If eng <>0, @code{std} is used for Groebner basis computations,
// @*       otherwise, and by default @code{slimgb} is used.
// DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
// @*          if @code{printlevel}>=2, all the debug messages will be printed.
// EXAMPLE: example Sannfslog; shows examples
// "
// {
//   int eng = 0;
//   if ( size(#)>0 )
//   {
//     if ( typeof(#[1]) == "int" )
//     {
//       eng = int(#[1]);
//     }
//   }
//   int ppl = printlevel-voice+2;
//   def save = basering;
//   int N = nvars(basering);
//   int Nnew = 2*N+1;
//   int i;
//   string s;
//   list RL = ringlist(basering);
//   list L, Lord;
//   list tmp;
//   intvec iv;
//   L[1] = RL[1]; // char
//   L[4] = RL[4]; // char, minpoly
//   // check whether vars have admissible names
//   list Name = RL[2];
//   for (i=1; i<=N; i++)
//   {
//     if (Name[i] == "s")
//     {
//       ERROR("Variable names should not include s");
//     }
//   }
//   // the ideal I
//   ideal I = -F, jacob(F);
//   dbprint(ppl,"// -1-1- starting the computation of syz(-F,_Dx(F))");
//   dbprint(ppl-1, I);
//   matrix M = syz(I);
//   M = transpose(M);  // it is more usefull working with columns
//   dbprint(ppl,"// -1-2- the module syz(-F,_Dx(F)) has been computed");
//   dbprint(ppl-1, M);
//   // ------------ the ring @R ------------
//   // _x, _Dx, s;  elim.ord for _x,_Dx.
//   // now, create the names for new vars
//   list DName;
//   for (i=1; i<=N; i++)
//   {
//     DName[i] = "D"+Name[i]; // concat
//   }
//   tmp[1] = "s";
//   list NName;
//   for (i=1; i<=N; i++)
//   {
//     NName[2*i-1] = Name[i];
//     NName[2*i] = DName[i];
//     //NName[2*i-1] = DName[i];
//     //NName[2*i] = Name[i];
//   }
//   NName[Nnew] = tmp[1];
//   L[2] = NName;
//   tmp = 0;
//   // block ord (a(1,1),a(0,0,1,1),...,dp);
//   //list("a",intvec(1,1)), list("a",intvec(0,0,1,1)), ...
//   tmp[1] = "a";  // string
//   for (i=1; i<=N; i++)
//   {
//     iv[2*i-1] = 1;
//     iv[2*i]   = 1;
//     tmp[2]    = iv;  iv = 0;  // intvec
//     Lord[i]   = tmp;
//   }
//   //list("dp",intvec(1,1,1,1,1,...))
//   s = "iv=";
//   for (i=1; i<=Nnew; i++)
//   {
//     s = s+"1,";
//   }
//   s[size(s)]=";";
//   execute(s);
//   kill s;
//   tmp[1] = "dp";  // string
//   tmp[2] = iv;    // intvec
//   Lord[N+1] = tmp;
//   //list("C",intvec(0))
//   tmp[1] = "C";  // string
//   iv = 0;
//   tmp[2] = iv;   // intvec
//   Lord[N+2] = tmp;
//   tmp = 0;
//   L[3]    = Lord;
//   // we are done with the list. Now add a Plural part
//   def @R@ = ring(L);
//   setring @R@;
//   matrix @D[Nnew][Nnew];
//   for (i=1; i<=N; i++)
//   {
//     @D[2*i-1,2*i]=1;
//     //@D[2*i-1,2*i]=-1;
//   }
//   def @R = nc_algebra(1,@D);
//   setring @R;
//   kill @R@;
//   dbprint(ppl,"// -2-1- the ring @R(_x,_Dx,s) is ready");
//   dbprint(ppl-1, @R);
//   matrix M = imap(save,M);
//   // now, create the vector [-s,_Dx]
//   vector v = [-s];  // now s is a variable
//   for (i=1; i<=N; i++)
//   {
//     v = v + var(2*i)*gen(i+1);
//     //v = v + var(2*i-1)*gen(i+1);
//   }
//   ideal J = ideal(M*v);
//   // make leadcoeffs positive
//   for (i=1; i<= ncols(J); i++)
//   {
//     if ( leadcoef(J[i])<0 )
//     {
//       J[i] = -J[i];
//     }
//   }
//   ideal LD1 = J;
//   kill J;
//   export LD1;
//   return(@R);
// }
// example
// {
//   "EXAMPLE:"; echo = 2;
//   ring r = 0,(x,y),Dp;
//   poly F = x^4+y^5+x*y^4;
//   printlevel = 0;
//   def A  = Sannfslog(F);
//   setring A;
//   LD1;
// }


// alternative code for SannfsBM, renamed from annfsBM to ALTannfsBM
// is superfluos - will not be included in the official documentation
static proc ALTannfsBM (poly F, list #)
"USAGE:  ALTannfsBM(f [,eng]);  f a poly, eng an optional int
RETURN:  ring
PURPOSE: compute the annihilator ideal of f^s in D[s], where D is the Weyl
@*   algebra, according to the algorithm by Briancon and Maisonobe
NOTE:  activate the output ring with the @code{setring} command. In this ring,
@*       - the ideal LD is the annihilator of f^s.
@*       If eng <>0, @code{std} is used for Groebner basis computations,
@*       otherwise, and by default @code{slimgb} is used.
DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
@*          if @code{printlevel}>=2, all the debug messages will be printed.
EXAMPLE: example ALTannfsBM; shows examples
"
{
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  // returns a list with a ring and an ideal LD in it
  int ppl = printlevel-voice+2;
  //  printf("plevel :%s, voice: %s",printlevel,voice);
  def save = basering;
  int N = nvars(basering);
  int Nnew = 2*N+2;
  int i,j;
  string s;
  list RL = ringlist(basering);
  list L, Lord;
  list tmp;
  intvec iv;
  L[1] = RL[1];  //char
  L[4] = RL[4];  //char, minpoly
  // check whether vars have admissible names
  list Name  = RL[2];
  list RName;
  RName[1] = "t";
  RName[2] = "s";
  for (i=1; i<=N; i++)
  {
    for(j=1; j<=size(RName); j++)
    {
      if (Name[i] == RName[j])
      {
        ERROR("Variable names should not include t,s");
      }
    }
  }
  // now, create the names for new vars
  list DName;
  for (i=1; i<=N; i++)
  {
    DName[i] = "D"+Name[i];  //concat
  }
  tmp[1] = "t";
  tmp[2] = "s";
  list NName = tmp + Name + DName;
  L[2]   = NName;
  // Name, Dname will be used further
  kill NName;
  // block ord (lp(2),dp);
  tmp[1]  = "lp"; // string
  iv      = 1,1;
  tmp[2]  = iv; //intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1]  = "dp"; // string
  s       = "iv=";
  for (i=1; i<=Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)]= ";";
  execute(s);
  kill s;
  tmp[2]    = iv;
  Lord[2]   = tmp;
  tmp[1]    = "C";
  iv        = 0;
  tmp[2]    = iv;
  Lord[3]   = tmp;
  tmp       = 0;
  L[3]      = Lord;
  // we are done with the list
  def @R@ = ring(L);
  setring @R@;
  matrix @D[Nnew][Nnew];
  @D[1,2]=t;
  for(i=1; i<=N; i++)
  {
    @D[2+i,N+2+i]=1;
  }
  // L[5] = matrix(UpOneMatrix(Nnew));
  // L[6] = @D;
  def @R = nc_algebra(1,@D);
  setring @R;
  kill @R@;
  dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready");
  dbprint(ppl-1, @R);
  // create the ideal I
  poly  F = imap(save,F);
  ideal I = t*F+s;
  poly p;
  for(i=1; i<=N; i++)
  {
    p = t;  //t
    p = diff(F,var(2+i))*p;
    I = I, var(N+2+i) + p;
  }
  // -------- the ideal I is ready ----------
  dbprint(ppl,"// -1-2- starting the elimination of t in @R");
  dbprint(ppl-1, I);
  ideal J = engine(I,eng);
  ideal K = nselect(J,1);
  kill I,J;
  dbprint(ppl,"// -1-3- t is eliminated");
  dbprint(ppl-1, K);  //K is without t
  // create Dn[s], where Dn is the ordinary Weyl Algebra, and put the result into it,
  // thus creating the ring @R2
  // keep: N, i,j,s, tmp, RL
  setring save;
  Nnew = 2*N+1;
  // list RL = ringlist(save);  //is defined earlier
  kill Lord, tmp, iv;
  L = 0;
  list Lord, tmp;
  intvec iv;
  L[1] = RL[1];
  L[4] = RL[4];  //char, minpoly
  // check whether vars have admissible names -> done earlier
  // list Name = RL[2]
  // DName is defined earlier
  tmp[1] = "s";
  list NName = Name + DName + tmp;
  L[2] = NName;
  // dp ordering;
  string s = "iv=";
  for (i=1; i<=Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)] = ";";
  execute(s);
  kill s;
  tmp     = 0;
  tmp[1]  = "dp";  //string
  tmp[2]  = iv;    //intvec
  Lord[1] = tmp;
  tmp[1]  = "C";
  iv      = 0;
  tmp[2]  = iv;
  Lord[2] = tmp;
  tmp     = 0;
  L[3]    = Lord;
  // we are done with the list
  // Add: Plural part
  def @R2@ = ring(L);
  setring @R2@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R2 = nc_algebra(1,@D);
  setring @R2;
  kill @R2@;
  dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,s) is ready");
  dbprint(ppl-1, @R2);
  ideal K = imap(@R,K);
  option(redSB);
  //dbprint(ppl,"// -2-2- the final cosmetic std");
  //K = engine(K,eng);  //std does the job too
  // total cleanup
  kill @R;
  ideal LD = K;
  export LD;
  return(@R2);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z,w),Dp;
  poly F = x^3+y^3+z^2*w;
  printlevel = 0;
  def A = ALTannfsBM(F);
  setring A;
  LD;
}

proc bernsteinBM(poly F, list #)
"USAGE:  bernsteinBM(f [,eng]);  f a poly, eng an optional int
RETURN:  list (of roots of the Bernstein polynomial b and their multiplicies)
PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface,
@* defined by f, according to the algorithm by Briancon and Maisonobe
NOTE:    If eng <>0, @code{std} is used for Groebner basis computations,
@*       otherwise, and by default @code{slimgb} is used.
DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
@*          if @code{printlevel}>=2, all the debug messages will be printed.
EXAMPLE: example bernsteinBM; shows examples
"
{
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  // returns a list with a ring and an ideal LD in it
  int ppl = printlevel-voice+2;
  //  printf("plevel :%s, voice: %s",printlevel,voice);
  def save = basering;
  int N = nvars(basering);
  int Nnew = 2*N+2;
  int i,j;
  string s;
  list RL = ringlist(basering);
  list L, Lord;
  list tmp;
  intvec iv;
  L[1] = RL[1];  //char
  L[4] = RL[4];  //char, minpoly
  // check whether vars have admissible names
  list Name  = RL[2];
  list RName;
  RName[1] = "t";
  RName[2] = "s";
  for (i=1; i<=N; i++)
  {
    for(j=1; j<=size(RName); j++)
    {
      if (Name[i] == RName[j])
      {
        ERROR("Variable names should not include t,s");
      }
    }
  }
  // now, create the names for new vars
  list DName;
  for (i=1; i<=N; i++)
  {
    DName[i] = "D"+Name[i];  //concat
  }
  tmp[1] = "t";
  tmp[2] = "s";
  list NName = tmp + Name + DName;
  L[2]   = NName;
  // Name, Dname will be used further
  kill NName;
  // block ord (lp(2),dp);
  tmp[1]  = "lp"; // string
  iv      = 1,1;
  tmp[2]  = iv; //intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1]  = "dp"; // string
  s       = "iv=";
  for (i=1; i<=Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)]= ";";
  execute(s);
  kill s;
  tmp[2]    = iv;
  Lord[2]   = tmp;
  tmp[1]    = "C";
  iv        = 0;
  tmp[2]    = iv;
  Lord[3]   = tmp;
  tmp       = 0;
  L[3]      = Lord;
  // we are done with the list
  def @R@ = ring(L);
  setring @R@;
  matrix @D[Nnew][Nnew];
  @D[1,2]=t;
  for(i=1; i<=N; i++)
  {
    @D[2+i,N+2+i]=1;
  }
  // L[5] = matrix(UpOneMatrix(Nnew));
  // L[6] = @D;
  def @R = nc_algebra(1,@D);
  setring @R;
  kill @R@;
  dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready");
  dbprint(ppl-1, @R);
  // create the ideal I
  poly  F = imap(save,F);
  ideal I = t*F+s;
  poly p;
  for(i=1; i<=N; i++)
  {
    p = t;  //t
    p = diff(F,var(2+i))*p;
    I = I, var(N+2+i) + p;
  }
  // -------- the ideal I is ready ----------
  dbprint(ppl,"// -1-2- starting the elimination of t in @R");
  dbprint(ppl-1, I);
  ideal J = engine(I,eng);
  ideal K = nselect(J,1);
  kill I,J;
  dbprint(ppl,"// -1-3- t is eliminated");
  dbprint(ppl-1, K);  //K is without t
  // ----------- the ring @R2 ------------
  // _x, _Dx,s;  elim.ord for _x,_Dx.
  // keep: N, i,j,s, tmp, RL
  setring save;
  Nnew = 2*N+1;
  kill Lord, tmp, iv, RName;
  list Lord, tmp;
  intvec iv;
  L[1] = RL[1];
  L[4] = RL[4];  //char, minpoly
  // check whether vars hava admissible names -> done earlier
  // now, create the names for new var
  tmp[1] = "s";
  // DName is defined earlier
  list NName = Name + DName + tmp;
  L[2] = NName;
  tmp = 0;
  // block ord (dp(N),dp);
  string s = "iv=";
  for (i=1; i<=Nnew-1; i++)
  {
    s = s+"1,";
  }
  s[size(s)]=";";
  execute(s);
  tmp[1] = "dp";  //string
  tmp[2] = iv;    //intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1] = "dp";  //string
  s[size(s)] = ",";
  s = s+"1;";
  execute(s);
  kill s;
  kill NName;
  tmp[2]  = iv;
  Lord[2] = tmp;
  tmp[1]  = "C";
  iv      = 0;
  tmp[2]  = iv;
  Lord[3] = tmp;
  tmp     = 0;
  L[3]    = Lord;
  // we are done with the list. Now add a Plural part
  def @R2@ = ring(L);
  setring @R2@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R2 = nc_algebra(1,@D);
  setring @R2;
  kill @R2@;
  dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,s) is ready");
  dbprint(ppl-1, @R2);
  ideal MM = maxideal(1);
  MM = 0,s,MM;
  map R01 = @R, MM;
  ideal K = R01(K);
  kill @R, R01;
  poly F = imap(save,F);
  K = K,F;
  dbprint(ppl,"// -2-2- starting the elimination of _x,_Dx in @R2");
  dbprint(ppl-1, K);
  ideal M = engine(K,eng);
  ideal K2 = nselect(M,1..Nnew-1);
  kill K,M;
  dbprint(ppl,"// -2-3- _x,_Dx are eliminated in @R2");
  dbprint(ppl-1, K2);
  // the ring @R3 and the search for minimal negative int s
  ring @R3 = 0,s,dp;
  dbprint(ppl,"// -3-1- the ring @R3(s) is ready");
  ideal K3 = imap(@R2,K2);
  kill @R2;
  poly p = K3[1];
  dbprint(ppl,"// -3-2- factorization");
  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];
  }
  // "--------- b-function factorizes into ---------";  P;
  //int sP = minIntRoot(bs,1);
  //dbprint(ppl,"// -3-3- minimal integer root found");
  //dbprint(ppl-1, sP);
  // convert factors to a list of their roots and multiplicities
  bs =  normalize(bs);
  bs = -subst(bs,s,0);
  setring save;
  ideal bs = imap(@R3,bs);
  kill @R3;
  list BS = bs,m;
  return(BS);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z,w),Dp;
  poly F = x^3+y^3+z^2*w;
  printlevel = 0;
  bernsteinBM(F);
}

// some changes
proc annfsBM (poly F, list #)
"USAGE:  annfsBM(f [,eng]);  f a poly, eng an optional int
RETURN:  ring
PURPOSE: compute the D-module structure of basering[1/f]*f^s, according
@* to the algorithm by Briancon and Maisonobe
NOTE:    activate the output ring with the @code{setring} command. In this ring,
@*       - the ideal LD (which is a Groebner basis) is the needed D-module structure,
@*         which is obtained by substituting the minimal integer root of a Bernstein
@*         polynomial into the s-parametric ideal;
@*       - the list BS is the list of roots and multiplicities of a Bernstein polynomial of f.
@*       If eng <>0, @code{std} is used for Groebner basis computations,
@*       otherwise, and by default @code{slimgb} is used.
DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
@*          if @code{printlevel}>=2, all the debug messages will be printed.
EXAMPLE: example annfsBM; shows examples
"
{
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  // returns a list with a ring and an ideal LD in it
  int ppl = printlevel-voice+2;
  //  printf("plevel :%s, voice: %s",printlevel,voice);
  def save = basering;
  int N = nvars(basering);
  int Nnew = 2*N+2;
  int i,j;
  string s;
  list RL = ringlist(basering);
  list L, Lord;
  list tmp;
  intvec iv;
  L[1] = RL[1];  //char
  L[4] = RL[4];  //char, minpoly
  // check whether vars have admissible names
  list Name  = RL[2];
  list RName;
  RName[1] = "t";
  RName[2] = "s";
  for (i=1; i<=N; i++)
  {
    for(j=1; j<=size(RName); j++)
    {
      if (Name[i] == RName[j])
      {
        ERROR("Variable names should not include t,s");
      }
    }
  }
  // now, create the names for new vars
  list DName;
  for (i=1; i<=N; i++)
  {
    DName[i] = "D"+Name[i];  //concat
  }
  tmp[1] = "t";
  tmp[2] = "s";
  list NName = tmp + Name + DName;
  L[2]   = NName;
  // Name, Dname will be used further
  kill NName;
  // block ord (lp(2),dp);
  tmp[1]  = "lp"; // string
  iv      = 1,1;
  tmp[2]  = iv; //intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1]  = "dp"; // string
  s       = "iv=";
  for (i=1; i<=Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)]= ";";
  execute(s);
  kill s;
  tmp[2]    = iv;
  Lord[2]   = tmp;
  tmp[1]    = "C";
  iv        = 0;
  tmp[2]    = iv;
  Lord[3]   = tmp;
  tmp       = 0;
  L[3]      = Lord;
  // we are done with the list
  def @R@ = ring(L);
  setring @R@;
  matrix @D[Nnew][Nnew];
  @D[1,2]=t;
  for(i=1; i<=N; i++)
  {
    @D[2+i,N+2+i]=1;
  }
  // L[5] = matrix(UpOneMatrix(Nnew));
  // L[6] = @D;
  def @R = nc_algebra(1,@D);
  setring @R;
  kill @R@;
  dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready");
  dbprint(ppl-1, @R);
  // create the ideal I
  poly  F = imap(save,F);
  ideal I = t*F+s;
  poly p;
  for(i=1; i<=N; i++)
  {
    p = t;  //t
    p = diff(F,var(2+i))*p;
    I = I, var(N+2+i) + p;
  }
  // -------- the ideal I is ready ----------
  dbprint(ppl,"// -1-2- starting the elimination of t in @R");
  dbprint(ppl-1, I);
  ideal J = engine(I,eng);
  ideal K = nselect(J,1);
  kill I,J;
  dbprint(ppl,"// -1-3- t is eliminated");
  dbprint(ppl-1, K);  //K is without t
  setring save;
  // ----------- the ring @R2 ------------
  // _x, _Dx,s;  elim.ord for _x,_Dx.
  // keep: N, i,j,s, tmp, RL
  Nnew = 2*N+1;
  kill Lord, tmp, iv, RName;
  list Lord, tmp;
  intvec iv;
  L[1] = RL[1];
  L[4] = RL[4];  //char, minpoly
  // check whether vars hava admissible names -> done earlier
  // now, create the names for new var
  tmp[1] = "s";
  // DName is defined earlier
  list NName = Name + DName + tmp;
  L[2] = NName;
  tmp = 0;
  // block ord (dp(N),dp);
  string s = "iv=";
  for (i=1; i<=Nnew-1; i++)
  {
    s = s+"1,";
  }
  s[size(s)]=";";
  execute(s);
  tmp[1] = "dp";  //string
  tmp[2] = iv;    //intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1] = "dp";  //string
  s[size(s)] = ",";
  s = s+"1;";
  execute(s);
  kill s;
  kill NName;
  tmp[2]  = iv;
  Lord[2] = tmp;
  tmp[1]  = "C";
  iv      = 0;
  tmp[2]  = iv;
  Lord[3] = tmp;
  tmp     = 0;
  L[3]    = Lord;
  // we are done with the list. Now add a Plural part
  def @R2@ = ring(L);
  setring @R2@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R2 = nc_algebra(1,@D);
  setring @R2;
  kill @R2@;
  dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,s) is ready");
  dbprint(ppl-1, @R2);
  ideal MM = maxideal(1);
  MM = 0,s,MM;
  map R01 = @R, MM;
  ideal K = R01(K);
  poly F = imap(save,F);
  K = K,F;
  dbprint(ppl,"// -2-2- starting the elimination of _x,_Dx in @R2");
  dbprint(ppl-1, K);
  ideal M = engine(K,eng);
  ideal K2 = nselect(M,1..Nnew-1);
  kill K,M;
  dbprint(ppl,"// -2-3- _x,_Dx are eliminated in @R2");
  dbprint(ppl-1, K2);
  // the ring @R3 and the search for minimal negative int s
  ring @R3 = 0,s,dp;
  dbprint(ppl,"// -3-1- the ring @R3(s) is ready");
  ideal K3 = imap(@R2,K2);
  poly p = K3[1];
  dbprint(ppl,"// -3-2- factorization");
  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 ignore P[1][1] (constant) and P[2][1] (its mult.)
  {
    bs[i-1] = P[1][i];
    m[i-1]  = P[2][i];
  }
  // "--------- b-function factorizes into ---------";  P;
  int sP = minIntRoot(bs,1);
  dbprint(ppl,"// -3-3- minimal integer root found");
  dbprint(ppl-1, sP);
  // convert factors to a list of their roots
  bs = normalize(bs);
  bs = -subst(bs,s,0);
  list BS = bs,m;
  //TODO: sort BS!
  // --------- substitute s found in the ideal ---------
  // --------- going back to @R and substitute ---------
  setring @R;
  ideal K2 = subst(K,s,sP);
  kill K;
  // create the ordinary Weyl algebra and put the result into it,
  // thus creating the ring @R4
  // keep: N, i,j,s, tmp, RL
  setring save;
  Nnew = 2*N;
  // list RL = ringlist(save);  //is defined earlier
  kill Lord, tmp, iv;
  L = 0;
  list Lord, tmp;
  intvec iv;
  L[1] = RL[1];
  L[4] = RL[4];  //char, minpoly
  // check whether vars have admissible names -> done earlier
  // list Name = RL[2]M
  // DName is defined earlier
  list NName = Name + DName;
  L[2] = NName;
  // dp ordering;
  string s = "iv=";
  for (i=1; i<=Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)] = ";";
  execute(s);
  kill s;
  tmp     = 0;
  tmp[1]  = "dp";  //string
  tmp[2]  = iv;    //intvec
  Lord[1] = tmp;
  tmp[1]  = "C";
  iv      = 0;
  tmp[2]  = iv;
  Lord[2] = tmp;
  tmp     = 0;
  L[3]    = Lord;
  // we are done with the list
  // Add: Plural part
  def @R4@ = ring(L);
  setring @R4@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R4 = nc_algebra(1,@D);
  setring @R4;
  kill @R4@;
  dbprint(ppl,"// -4-1- the ring @R4(_x,_Dx) is ready");
  dbprint(ppl-1, @R4);
  ideal K4 = imap(@R,K2);
  intvec saveopt=option(get);
  option(redSB);
  dbprint(ppl,"// -4-2- the final cosmetic std");
  K4 = engine(K4,eng);  //std does the job too
  // total cleanup
  kill @R;
  kill @R2;
  list BS = imap(@R3,BS);
  export BS;
  kill @R3;
  ideal LD = K4;
  export LD;
  option(set,saveopt);
  return(@R4);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),Dp;
  poly F = z*x^2+y^3;
  printlevel = 0;
  def A = annfsBM(F);
  setring A;
  LD;
  BS;
}


// replacing s with -s-1 => data is shorter
// analogue of annfs0
proc annfs2(ideal I, poly F, list #)
"USAGE:  annfs2(I, F [,eng]);  I an ideal, F a poly, eng an optional int
RETURN:  ring
PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra,
@*       based on the output of Sannfs-like procedure
@*       annfs2 uses shorter expressions in the variable s (the idea of Noro).
NOTE: activate the output ring with the @code{setring} command. In this ring,
@*      - the ideal LD (which is a Groebner basis) is the annihilator of f^s,
@*      - the list BS contains the roots with multiplicities of the BS polynomial.
@*    If eng <>0, @code{std} is used for Groebner basis computations,
@*    otherwise and by default @code{slimgb} is used.
DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
@*       if @code{printlevel}>=2, all the debug messages will be printed.
EXAMPLE: example annfs2; shows examples
"
{
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  def @R2 = basering;
  // we're in D_n[s], where the elim ord for s is set
  ideal J = NF(I,std(F));
  // make leadcoeffs positive
  int i;
  J = subst(J,s,-s-1);
  for (i=1; i<= ncols(J); i++)
  {
    if (leadcoef(J[i]) <0 )
    {
      J[i] = -J[i];
    }
  }
  J = J,F;
  ideal M = engine(J,eng);
  int Nnew = nvars(@R2);
  ideal K2 = nselect(M,1..Nnew-1);
  int ppl = printlevel-voice+2;
  dbprint(ppl,"// -1-1- _x,_Dx are eliminated in basering");
  dbprint(ppl-1, K2);
  // the ring @R3 and the search for minimal negative int s
  ring @R3 = 0,s,dp;
  dbprint(ppl,"// -2-1- the ring @R3 i.e. K[s] is ready");
  ideal K3 = imap(@R2,K2);
  poly p = K3[1];
  dbprint(ppl,"// -2-2- factorization");
  //  ideal P = factorize(p,1);  //without constants and multiplicities
  //  "--------- b-function factorizes into ---------";   P;
  // convert factors to the list of their roots with mults
  // assume all factors are linear
  //  ideal BS = normalize(P);
  //  BS = subst(BS,s,0);
  //  BS = -BS;
  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]; bs[i-1] = subst(bs[i-1],s,-s-1);
    m[i-1]  = P[2][i];
  }
  int sP = minIntRoot(bs,1);
  bs =  normalize(bs);
  bs = -subst(bs,s,0);
  dbprint(ppl,"// -2-3- minimal integer root found");
  dbprint(ppl-1, sP);
 //TODO: sort BS!
  // --------- substitute s found in the ideal ---------
  // --------- going back to @R and substitute ---------
  setring @R2;
  K2 = subst(I,s,sP);
  // create the ordinary Weyl algebra and put the result into it,
  // thus creating the ring @R5
  // keep: N, i,j,s, tmp, RL
  Nnew = Nnew - 1; // former 2*N;
  // list RL = ringlist(save);  // is defined earlier
  //  kill Lord, tmp, iv;
  list L = 0;
  list Lord, tmp;
  intvec iv;
  list RL = ringlist(basering);
  L[1] = RL[1];
  L[4] = RL[4];  //char, minpoly
  // check whether vars have admissible names -> done earlier
  // list Name = RL[2]M
  // DName is defined earlier
  list NName; // = RL[2]; // skip the last var 's'
  for (i=1; i<=Nnew; i++)
  {
    NName[i] =  RL[2][i];
  }
  L[2] = NName;
  // dp ordering;
  string s = "iv=";
  for (i=1; i<=Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)] = ";";
  execute(s);
  tmp     = 0;
  tmp[1]  = "dp";  // string
  tmp[2]  = iv;  // intvec
  Lord[1] = tmp;
  kill s;
  tmp[1]  = "C";
  iv = 0;
  tmp[2]  = iv;
  Lord[2] = tmp;
  tmp     = 0;
  L[3]    = Lord;
  // we are done with the list
  // Add: Plural part
  def @R4@ = ring(L);
  setring @R4@;
  int N = Nnew div 2;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R4 = nc_algebra(1,@D);
  setring @R4;
  kill @R4@;
  dbprint(ppl,"// -3-1- the ring @R4 is ready");
  dbprint(ppl-1, @R4);
  ideal K4 = imap(@R2,K2);
  option(redSB);
  dbprint(ppl,"// -3-2- the final cosmetic std");
  K4 = engine(K4,eng);  // std does the job too
  // total cleanup
  ideal bs = imap(@R3,bs);
  kill @R3;
  list BS = bs,m;
  export BS;
  ideal LD = K4;
  export LD;
  return(@R4);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),Dp;
  poly F = x^3+y^3+z^3;
  printlevel = 0;
  def A = SannfsBM(F);
  setring A;
  LD;
  poly F = imap(r,F);
  def B  = annfs2(LD,F);
  setring B;
  LD;
  BS;
}

// try to replace s with -s-1 => data is shorter as in annfs2
// and use what Macaulay2 people call reduceB strategy, that is add
// not F but Tjurina ideal <F,dF/dx1,...,dF/dxN>; the resulting B-function
// has to be multiplied with (s+1) at the very end
proc annfsRB(ideal I, poly F, list #)
"USAGE:  annfsRB(I, F [,eng]);  I an ideal, F a poly, eng an optional int
RETURN:  ring
PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra,
@* based on the output of Sannfs like procedure
NOTE:    activate the output ring with the @code{setring} command. In this ring,
@*       - the ideal LD (which is a Groebner basis) is the annihilator of f^s,
@*       - the list BS contains the roots with multiplicities of a Bernstein polynomial of f.
@*       If eng <>0, @code{std} is used for Groebner basis computations,
@*       otherwise and by default @code{slimgb} is used.
@*       This procedure uses in addition to F its Jacobian ideal.
DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
@*       if @code{printlevel}>=2, all the debug messages will be printed.
EXAMPLE: example annfsRB; shows examples
"
{
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  def @R2 = basering;
  int ppl = printlevel-voice+2;
  // we're in D_n[s], where the elim ord for s is set
  // switch to comm. ring in X's and compute the GB of Tjurina ideal
  dbprint(ppl,"// -1-0- creating K[x] and Tjurina ideal");
  list nL = ringlist(@R2);
  list temp,t2;
  temp[1] = nL[1];
  temp[4] = nL[4];
  int @n = int((nvars(@R2)-1) div 2); // # of x's
  int i;
  for (i=1; i<=@n; i++)
  {
    t2[i] = nL[2][i];
  }
  temp[2] = t2;
  t2 = 0;
  t2[1] = nL[3][1]; // more weights than vars?
  t2[2] = nL[3][3];
  temp[3] = t2;
  def @R22 = ring(temp);
  setring @R22;
  poly F = imap(@R2,F);
  ideal J = F;
  for (i=1; i<=@n; i++)
  {
    J = J, diff(F,var(i));
  }
  J = engine(J,eng);
  dbprint(ppl,"// -1-1- finished computing the GB of Tjurina ideal");
  dbprint(ppl-1, J);
  setring @R2;
  ideal JF = imap(@R22,J);
  kill @R22;
  attrib(JF,"isSB",1); // embedded comm ring is used
  ideal J = NF(I,JF);
  dbprint(ppl,"// -1-2- finished computing the NF of I w.r.t. Tjurina ideal");
  dbprint(ppl-1, J2);
  // make leadcoeffs positive
  J = subst(J,s,-s-1);
  for (i=1; i<= ncols(J); i++)
  {
    if (leadcoef(J[i]) <0 )
    {
      J[i] = -J[i];
    }
  }
  J = J,JF;
  ideal M = engine(J,eng);
  int Nnew = nvars(@R2);
  ideal K2 = nselect(M,1..Nnew-1);
  dbprint(ppl,"// -2-0- _x,_Dx are eliminated in basering");
  dbprint(ppl-1, K2);
  // the ring @R3 and the search for minimal negative int s
  ring @R3 = 0,s,dp;
  dbprint(ppl,"// -2-1- the ring @R3 i.e. K[s] is ready");
  ideal K3 = imap(@R2,K2);
  poly p = K3[1];
  p = s*p; // mult with the lost (s+1) factor
  dbprint(ppl,"// -2-2- factorization");
  //  ideal P = factorize(p,1);  //without constants and multiplicities
  //  "--------- b-function factorizes into ---------";   P;
  // convert factors to the list of their roots with mults
  // assume all factors are linear
  //  ideal BS = normalize(P);
  //  BS = subst(BS,s,0);
  //  BS = -BS;
  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]; bs[i-1] = subst(bs[i-1],s,-s-1);
    m[i-1]  = P[2][i];
  }
  int sP = minIntRoot(bs,1);
  bs =  normalize(bs);
  bs = -subst(bs,s,0);
  dbprint(ppl,"// -2-3- minimal integer root found");
  dbprint(ppl-1, sP);
 //TODO: sort BS!
  // --------- substitute s found in the ideal ---------
  // --------- going back to @R and substitute ---------
  setring @R2;
  K2 = subst(I,s,sP);
  // create the ordinary Weyl algebra and put the result into it,
  // thus creating the ring @R5
  // keep: N, i,j,s, tmp, RL
  Nnew = Nnew - 1; // former 2*N;
  // list RL = ringlist(save);  // is defined earlier
  //  kill Lord, tmp, iv;
  list L = 0;
  list Lord, tmp;
  intvec iv;
  list RL = ringlist(basering);
  L[1] = RL[1];
  L[4] = RL[4];  //char, minpoly
  // check whether vars have admissible names -> done earlier
  // list Name = RL[2]M
  // DName is defined earlier
  list NName; // = RL[2]; // skip the last var 's'
  for (i=1; i<=Nnew; i++)
  {
    NName[i] =  RL[2][i];
  }
  L[2] = NName;
  // dp ordering;
  string s = "iv=";
  for (i=1; i<=Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)] = ";";
  execute(s);
  tmp     = 0;
  tmp[1]  = "dp";  // string
  tmp[2]  = iv;  // intvec
  Lord[1] = tmp;
  kill s;
  tmp[1]  = "C";
  iv = 0;
  tmp[2]  = iv;
  Lord[2] = tmp;
  tmp     = 0;
  L[3]    = Lord;
  // we are done with the list
  // Add: Plural part
  def @R4@ = ring(L);
  setring @R4@;
  int N = Nnew div 2;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R4 = nc_algebra(1,@D);
  setring @R4;
  kill @R4@;
  dbprint(ppl,"// -3-1- the ring @R4 is ready");
  dbprint(ppl-1, @R4);
  ideal K4 = imap(@R2,K2);
  intvec saveopt=option(get);
  option(redSB);
  dbprint(ppl,"// -3-2- the final cosmetic std");
  K4 = engine(K4,eng);  // std does the job too
  // total cleanup
  ideal bs = imap(@R3,bs);
  kill @R3;
  list BS = bs,m;
  export BS;
  ideal LD = K4;
  export LD;
  option(set,saveopt);
  return(@R4);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),Dp;
  poly F = x^3+y^3+z^3;
  printlevel = 0;
  def A = SannfsBM(F); setring A;
  LD; // s-parametric ahhinilator
  poly F = imap(r,F);
  def B  = annfsRB(LD,F); setring B;
  LD;
  BS;
}

proc operatorBM(poly F, list #)
"USAGE:  operatorBM(f [,eng]);  f a poly, eng an optional int
RETURN:  ring
PURPOSE: compute the B-operator and other relevant data for Ann F^s,
@*  using e.g. algorithm by Briancon and Maisonobe for Ann F^s and BS.
NOTE:    activate the output ring with the @code{setring} command. In the output ring D[s]
@*       - the polynomial F is the same as the input,
@*       - the ideal LD is the annihilator of f^s in Dn[s],
@*       - the ideal LD0 is the needed D-mod structure, where LD0 = LD|s=s0,
@*       - the polynomial bs is the global Bernstein polynomial of f in the variable s,
@*       - the list BS contains all the roots with multiplicities of the global Bernstein polynomial of f,
@*       - the polynomial PS is an operator in Dn[s] such that PS*f^(s+1) = bs*f^s.
@*       If eng <>0, @code{std} is used for Groebner basis computations,
@*       otherwise and by default @code{slimgb} is used.
DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
@*       if @code{printlevel}>=2, all the debug messages will be printed.
EXAMPLE: example operatorBM; shows examples
"
{
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  // returns a list with a ring and an ideal LD in it
  int ppl = printlevel-voice+2;
  //  printf("plevel :%s, voice: %s",printlevel,voice);
  def save = basering;
  int N = nvars(basering);
  int Nnew = 2*N+2;
  int i,j;
  string s;
  list RL = ringlist(basering);
  list L, Lord;
  list tmp;
  intvec iv;
  L[1] = RL[1];  //char
  L[4] = RL[4];  //char, minpoly
  // check whether vars have admissible names
  list Name  = RL[2];
  list RName;
  RName[1] = "t";
  RName[2] = "s";
  for (i=1; i<=N; i++)
  {
    for(j=1; j<=size(RName); j++)
    {
      if (Name[i] == RName[j])
      {
        ERROR("Variable names should not include t,s");
      }
    }
  }
  // now, create the names for new vars
  list DName;
  for (i=1; i<=N; i++)
  {
    DName[i] = "D"+Name[i];  //concat
  }
  tmp[1] = "t";
  tmp[2] = "s";
  list NName = tmp + Name + DName;
  L[2]   = NName;
  // Name, Dname will be used further
  kill NName;
  // block ord (lp(2),dp);
  tmp[1]  = "lp"; // string
  iv      = 1,1;
  tmp[2]  = iv; //intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1]  = "dp"; // string
  s       = "iv=";
  for (i=1; i<=Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)]= ";";
  execute(s);
  kill s;
  tmp[2]    = iv;
  Lord[2]   = tmp;
  tmp[1]    = "C";
  iv        = 0;
  tmp[2]    = iv;
  Lord[3]   = tmp;
  tmp       = 0;
  L[3]      = Lord;
  // we are done with the list
  def @R@ = ring(L);
  setring @R@;
  matrix @D[Nnew][Nnew];
  @D[1,2]=t;
  for(i=1; i<=N; i++)
  {
    @D[2+i,N+2+i]=1;
  }
  // L[5] = matrix(UpOneMatrix(Nnew));
  // L[6] = @D;
  def @R = nc_algebra(1,@D);
  setring @R;
  kill @R@;
  dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready");
  dbprint(ppl-1, @R);
  // create the ideal I
  poly  F = imap(save,F);
  ideal I = t*F+s;
  poly p;
  for(i=1; i<=N; i++)
  {
    p = t;  //t
    p = diff(F,var(2+i))*p;
    I = I, var(N+2+i) + p;
  }
  // -------- the ideal I is ready ----------
  dbprint(ppl,"// -1-2- starting the elimination of t in @R");
  dbprint(ppl-1, I);
  ideal J = engine(I,eng);
  ideal K = nselect(J,1);
  kill I,J;
  dbprint(ppl,"// -1-3- t is eliminated");
  dbprint(ppl-1, K);  //K is without t
  setring save;
  // ----------- the ring @R2 ------------
  // _x, _Dx,s;  elim.ord for _x,_Dx.
  // keep: N, i,j,s, tmp, RL
  Nnew = 2*N+1;
  kill Lord, tmp, iv, RName;
  list Lord, tmp;
  intvec iv;
  L[1] = RL[1];
  L[4] = RL[4];  //char, minpoly
  // check whether vars hava admissible names -> done earlier
  // now, create the names for new var
  tmp[1] = "s";
  // DName is defined earlier
  list NName = Name + DName + tmp;
  L[2] = NName;
  tmp = 0;
  // block ord (dp(N),dp);
  string s = "iv=";
  for (i=1; i<=Nnew-1; i++)
  {
    s = s+"1,";
  }
  s[size(s)]=";";
  execute(s);
  tmp[1] = "dp";  //string
  tmp[2] = iv;    //intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1] = "dp";  //string
  s[size(s)] = ",";
  s = s+"1;";
  execute(s);
  kill s;
  kill NName;
  tmp[2]  = iv;
  Lord[2] = tmp;
  tmp[1]  = "C";
  iv      = 0;
  tmp[2]  = iv;
  Lord[3] = tmp;
  tmp     = 0;
  L[3]    = Lord;
  // we are done with the list. Now add a Plural part
  def @R2@ = ring(L);
  setring @R2@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R2 = nc_algebra(1,@D);
  setring @R2;
  kill @R2@;
  dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,s) is ready");
  dbprint(ppl-1, @R2);
  ideal MM = maxideal(1);
  MM = 0,s,MM;
  map R01 = @R, MM;
  ideal K = R01(K);
  poly F = imap(save,F);
  K = K,F;
  dbprint(ppl,"// -2-2- starting the elimination of _x,_Dx in @R2");
  dbprint(ppl-1, K);
  ideal M = engine(K,eng);
  ideal K2 = nselect(M,1..Nnew-1);
  kill K,M;
  dbprint(ppl,"// -2-3- _x,_Dx are eliminated in @R2");
  dbprint(ppl-1, K2);
  // the ring @R3 and the search for minimal negative int s
  ring @R3 = 0,s,dp;
  dbprint(ppl,"// -3-1- the ring @R3(s) is ready");
  ideal K3 = imap(@R2,K2);
  kill @R2;
  poly p = K3[1];
  dbprint(ppl,"// -3-2- factorization");
  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];
  }
  // "--------- b-function factorizes into ---------";  P;
  int sP = minIntRoot(bs,1);
  dbprint(ppl,"// -3-3- minimal integer root found");
  dbprint(ppl-1, sP);
  // convert factors to a list of their roots with multiplicities
  bs = normalize(bs);
  bs = -subst(bs,s,0);
  list BS = bs,m;
  //TODO: sort BS!
  // --------- substitute s found in the ideal ---------
  // --------- going back to @R and substitute ---------
  setring @R;
  ideal K2 = subst(K,s,sP);
  // create Dn[s], where Dn is the ordinary Weyl algebra, and put the result into it,
  // thus creating the ring @R4
  // keep: N, i,j,s, tmp, RL
  setring save;
  Nnew = 2*N+1;
  // list RL = ringlist(save);  //is defined earlier
  kill Lord, tmp, iv;
  L = 0;
  list Lord, tmp;
  intvec iv;
  L[1] = RL[1];
  L[4] = RL[4];  //char, minpoly
  // check whether vars have admissible names -> done earlier
  // list Name = RL[2]
  // DName is defined earlier
  tmp[1] = "s";
  list NName = Name + DName + tmp;
  L[2] = NName;
  // dp ordering;
  string s = "iv=";
  for (i=1; i<=Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)] = ";";
  execute(s);
  kill s;
  tmp     = 0;
  tmp[1]  = "dp";  //string
  tmp[2]  = iv;    //intvec
  Lord[1] = tmp;
  tmp[1]  = "C";
  iv      = 0;
  tmp[2]  = iv;
  Lord[2] = tmp;
  tmp     = 0;
  L[3]    = Lord;
  // we are done with the list
  // Add: Plural part
  def @R4@ = ring(L);
  setring @R4@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R4 = nc_algebra(1,@D);
  setring @R4;
  kill @R4@;
  dbprint(ppl,"// -4-1- the ring @R4(_x,_Dx,s) is ready");
  dbprint(ppl-1, @R4);
  ideal LD0 = imap(@R,K2);
  ideal LD  = imap(@R,K);
  kill @R;
  poly bs = imap(@R3,p);
  list BS = imap(@R3,BS);
  kill @R3;
  bs = normalize(bs);
  poly F = imap(save,F);
  dbprint(ppl,"// -4-2- starting the computation of PS via lift");
//better liftstd, I didn't knot it works also for Plural, liftslimgb?
// liftstd may give extra coeffs in the resulting ideal
  matrix T = lift(F+LD,bs);
  poly PS = T[1,1];
  dbprint(ppl,"// -4-3- an operator PS found, PS*f^(s+1) = b(s)*f^s");
  dbprint(ppl-1,PS);
  intvec saveopt=option(get);
  option(redSB);
  dbprint(ppl,"// -4-4- the final cosmetic std");
  LD0 = engine(LD0,eng);  //std does the job too
  LD  = engine(LD,eng);
  export F,LD,LD0,bs,BS,PS;
  option(set,saveopt);
  return(@R4);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),Dp;
  poly F = x^3+y^3+z^3;
  printlevel = 0;
  def A = operatorBM(F);
  setring A;
  F; // the original polynomial itself
  LD; // generic annihilator
  LD0; // annihilator
  bs; // normalized Bernstein poly
  BS; // roots and multiplicities of the Bernstein poly
  PS; // the operator, s.t. PS*F^{s+1} = bs*F^s mod LD
  reduce(PS*F-bs,LD); // check the property of PS
}

// more interesting:
//  ring r = 0,(x,y,z,w),Dp;
//  poly F = x^3+y^3+z^2*w;

// need: (c,<) ordering for such comp's

proc operatorModulo(poly F, ideal I, poly b)
"USAGE:  operatorModulo(f,I,b);  f a poly, I an ideal, b a poly
RETURN:  poly
PURPOSE: compute the B-operator from the polynomial f,
@*           ideal I = Ann f^s and Bernstein-Sato polynomial b
@*           using modulo i.e. kernel of module homomorphism
NOTE:  The computations take place in the ring, similar to the one
@*       returned by Sannfs procedure.
@*     Note, that operator is not completely reduced wrt Ann f^{s+1}.
@*     If printlevel=1, progress debug messages will be printed,
@*     if printlevel>=2, all the debug messages will be printed.
EXAMPLE: example operatorModulo; shows examples
"
{
  int ppl = printlevel-voice+2;
  def save = basering;
  // change the ordering on the currRing
  def mering = makeModElimRing(save);
  setring mering;
  poly b = imap(save, b);
  poly F = imap(save, F);
  ideal I = imap(save, I);
  matrix N = matrix(I); // ann f^s
  //  matrix K = hom_kernel(AA,M,N);
  // option(noreturnSB)?
  /// matrix K = modulo(AA,N); // too slow: do it with slim!
    module M = b,-F;
    dbprint(ppl,"starting modulo computation");
    module K = moduloSlim(M,N);
    dbprint(ppl,"finished modulo computation");
    //    K = transpose(K);
//     matrix M[3][s+2] = F,-b,I[1..s], 1,0:(s+1),0,1,0:(s);
//     module GM = slimgb(M);
//     module GMT = transpose(G);
//     GMT = GMT[2],GMT[3]; // modulo matrix
//     module K = GMT[2];
//     GMT = transpose(GMT);
//     K = transpose(K);
//     matrix K = GMT;
//////////////////////////////////////////////////
//  now select those elts whose 0's entry is nonzero
// if there is constant => done
// if not => compute GB and get it
    module L;
    ideal J;
    int i;
    poly t; number n;
    for(i=1; i<=ncols(K); i++)
    {
      if (K[1,i]!=0)
      {
        L = L,K[i];
        if ( leadmonom(K[1,i]) == 1)
        {
          t = K[2,i];
          n = leadcoef(K[1,i]);
          t = t/n;
          break;
          //          return(t);
        }
      }
    }
    if (n!=0)
    {
      // constant found
      setring save; poly t = imap(mering,t); kill mering;
      return(t);
    }
    dbprint(ppl,"no explicit constant. Start one more GB computation");
    // else: compute GB and do the same
    L = L[2..ncols(L)];
    K  = slimgb(L);
    dbprint(ppl,"finished GB computation");
    for(i=1; i<=ncols(K); i++)
    {
      if (K[1,i]!=0)
      {
        if ( leadmonom(K[1,i]) == 1)
        {
          t = K[2,i];
          n = leadcoef(K[1,i]);
          t = t/n;
//          break;
          setring save; poly t = imap(mering,t); kill mering;
          return(t);
        }
      }
    }

    // we are here if no constant found
    "ERROR: must never get here!";
    return(0);
//     for(i=1; i<=nrows(K); i++)
//     {
//       if (K[i,2]!=0)
//       {
//         if ( leadmonom(K[i,2]) == 1)
//         {
//           t = K[i,1];
//           n = leadcoef(K[i,2]);
//           t = t/n;
// //        J = J, K[i][2];
//           break;
//         }
//      }
//   }
//   ideal J = groebner(subst(I,s,s+1)); // for NF
//   t = NF(t,J);
//   "candidate:"; t;
//   J = subst(J,s,s-1);
//   // test:
//   if ( NF(t*F-b,J) !=0)
//   {
//     "Problem: PS does not work on F";
//   }
//   return(t);
}
example
{
  "EXAMPLE:"; echo = 2;
  //  LIB "dmod.lib"; option(prot); option(mem);
  ring r = 0,(x,y),Dp;
  poly F = x^3+y^3+x*y^3;
  def A = Sannfs(F); // here we get LD = ann f^s
  setring A;
  poly F = imap(r,F);
  def B = annfs0(LD,F); // to obtain BS polynomial
  list BS = imap(B,BS);   poly bs = fl2poly(BS,"s");
  poly PS = operatorModulo(F,LD,bs);
  LD = groebner(LD);
  PS = NF(PS,subst(LD,s,s+1));; // reduction modulo Ann s^{s+1}
  size(PS);
  lead(PS);
  reduce(PS*F-bs,LD); // check the defining property of PS
}

proc annfsParamBM (poly F, list #)
"USAGE:  annfsParamBM(f [,eng]);  f a poly, eng an optional int
RETURN:  ring
PURPOSE: compute the generic Ann F^s and exceptional parametric constellations
@* of a polynomial with parametric coefficients with the BM algorithm
NOTE:    activate the output ring with the @code{setring} command. In this ring,
@*       - the ideal LD is the D-module structure oa Ann F^s
@*       - the ideal Param contains special parameters as entries
@*       If eng <>0, @code{std} is used for Groebner basis computations,
@*       otherwise, and by default @code{slimgb} is used.
DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
@*          if @code{printlevel}>=2, all the debug messages will be printed.
EXAMPLE: example annfsParamBM; shows examples
"
{
  //PURPOSE: compute the list of all possible Bernstein-Sato polynomials for a polynomial with parametric coefficients, according to the algorithm by Briancon and Maisonobe
  // @*       - the list BS is the list of roots and multiplicities of a Bernstein polynomial of f.
  // ***** not implented yet ****
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  // returns a list with a ring and an ideal LD in it
  int ppl = printlevel-voice+2;
  //  printf("plevel :%s, voice: %s",printlevel,voice);
  def save = basering;
  int N = nvars(basering);
  int Nnew = 2*N+2;
  int i,j;
  string s;
  list RL = ringlist(basering);
  list L, Lord;
  list tmp;
  intvec iv;
  L[1] = RL[1];  //char
  L[4] = RL[4];  //char, minpoly
  // check whether vars have admissible names
  list Name  = RL[2];
  list RName;
  RName[1] = "t";
  RName[2] = "s";
  for (i=1; i<=N; i++)
  {
    for(j=1; j<=size(RName); j++)
    {
      if (Name[i] == RName[j])
      {
        ERROR("Variable names should not include t,s");
      }
    }
  }
  // now, create the names for new vars
  list DName;
  for (i=1; i<=N; i++)
  {
    DName[i] = "D"+Name[i];  //concat
  }
  tmp[1] = "t";
  tmp[2] = "s";
  list NName = tmp + Name + DName;
  L[2]   = NName;
  // Name, Dname will be used further
  kill NName;
  // block ord (lp(2),dp);
  tmp[1]  = "lp"; // string
  iv      = 1,1;
  tmp[2]  = iv; //intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1]  = "dp"; // string
  s       = "iv=";
  for (i=1; i<=Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)]= ";";
  execute(s);
  kill s;
  tmp[2]    = iv;
  Lord[2]   = tmp;
  tmp[1]    = "C";
  iv        = 0;
  tmp[2]    = iv;
  Lord[3]   = tmp;
  tmp       = 0;
  L[3]      = Lord;
  // we are done with the list
  def @R@ = ring(L);
  setring @R@;
  matrix @D[Nnew][Nnew];
  @D[1,2]=t;
  for(i=1; i<=N; i++)
  {
    @D[2+i,N+2+i]=1;
  }
  // L[5] = matrix(UpOneMatrix(Nnew));
  // L[6] = @D;
  def @R = nc_algebra(1,@D);
  setring @R;
  kill @R@;
  dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready");
  dbprint(ppl-1, @R);
  // create the ideal I
  poly  F = imap(save,F);
  ideal I = t*F+s;
  poly p;
  for(i=1; i<=N; i++)
  {
    p = t;  //t
    p = diff(F,var(2+i))*p;
    I = I, var(N+2+i) + p;
  }
  // -------- the ideal I is ready ----------
  dbprint(ppl,"// -1-2- starting the elimination of t in @R");
  dbprint(ppl-1, I);
  ideal J = engine(I,eng);
  ideal K = nselect(J,1);
  dbprint(ppl,"// -1-3- t is eliminated");
  dbprint(ppl-1, K);  //K is without t
  // ----- looking for special parameters -----
  dbprint(ppl,"// -2-1- starting the computation of the transformation matrix (via lift)");
  J = normalize(J);
  matrix T = lift(I,J);  //try also with liftstd
  kill I,J;
  dbprint(ppl,"// -2-2- the transformation matrix has been computed");
  dbprint(ppl-1, T);  //T is the transformation matrix
  dbprint(ppl,"// -2-3- genericity does the job");
  list lParam = genericity(T);
  int ip = size(lParam);
  int cip;
  string sParam;
  if (sParam[1]=="-") { sParam=""; } //genericity returns "-"
  // if no parameters exist in a basering
  for (cip=1; cip <= ip; cip++)
  {
    sParam = sParam + "," +lParam[cip];
  }
  if (size(sParam) >=2)
  {
    sParam = sParam[2..size(sParam)]; // removes the 1st colon
  }
  export sParam;
  kill T;
  dbprint(ppl,"// -2-4- the special parameters has been computed");
  dbprint(ppl, sParam);
  // create Dn[s], where Dn is the ordinary Weyl Algebra, and put the result into it,
  // thus creating the ring @R2
  // keep: N, i,j,s, tmp, RL
  setring save;
  Nnew = 2*N+1;
  // list RL = ringlist(save);  //is defined earlier
  kill Lord, tmp, iv;
  L = 0;
  list Lord, tmp;
  intvec iv;
  L[1] = RL[1];
  L[4] = RL[4];  //char, minpoly
  // check whether vars have admissible names -> done earlier
  // list Name = RL[2]M
  // DName is defined earlier
  tmp[1] = "s";
  list NName = Name + DName + tmp;
  L[2] = NName;
  // dp ordering;
  string s = "iv=";
  for (i=1; i<=Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)] = ";";
  execute(s);
  kill s;
  tmp     = 0;
  tmp[1]  = "dp";  //string
  tmp[2]  = iv;    //intvec
  Lord[1] = tmp;
  tmp[1]  = "C";
  iv      = 0;
  tmp[2]  = iv;
  Lord[2] = tmp;
  tmp     = 0;
  L[3]    = Lord;
  // we are done with the list
  // Add: Plural part
  def @R2@ = ring(L);
  setring @R2@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R2 = nc_algebra(1,@D);
  setring @R2;
  kill @R2@;
  dbprint(ppl,"// -3-1- the ring @R2(_x,_Dx,s) is ready");
  dbprint(ppl-1, @R2);
  ideal K = imap(@R,K);
  kill @R;
  intvec saveopt=option(get);
  option(redSB);
  dbprint(ppl,"// -3-2- the final cosmetic std");
  K = engine(K,eng);  //std does the job too
  ideal LD = K;
  export LD;
  if (sParam[1] == ",")
  {
    sParam = sParam[2..size(sParam)];
  }
  //  || ((sParam[1] == " ") && (sParam[2] == ",")))
  execute("ideal Param ="+sParam+";");
  export Param;
  kill sParam;
  option(set,saveopt);
  return(@R2);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = (0,a,b),(x,y),Dp;
  poly F = x^2 - (y-a)*(y-b);
  printlevel = 0;
  def A = annfsParamBM(F);  setring A;
  LD;
  Param;
  setring r;
  poly G = x2-(y-a)^2; // try the exceptional value b=a of parameters
  def B = annfsParamBM(G); setring B;
  LD;
  Param;
}

// *** the following example is nice, but too complicated for the documentation ***
//   ring r = (0,a),(x,y,z),Dp;
//   poly F = x^4+y^4+z^2+a*x*y*z;
//   printlevel = 2; //0
//   def A = annfsParamBM(F);
//   setring A;
//   LD;
//   Param;


proc annfsBMI(ideal F, list #)
"USAGE:  annfsBMI(F [,eng,met,us]);  F an ideal, eng, met, us optional ints
RETURN:  ring
PURPOSE: compute two kinds of Bernstein-Sato ideals, associated to
@* f = F[1]*..*F[P], with the multivariate algorithm by Briancon and Maisonobe.
NOTE:    activate the output ring with the @code{setring} command. In this ring,
@*       - the ideal LD is the annihilator of F[1]^s_1*..*F[P]^s_p,
@*       - the list or ideal BS is a Bernstein-Sato ideal of a polynomial f = F[1]*..*F[P].
@*       If eng <>0, @code{std} is used for Groebner basis computations,
@*       otherwise, and by default @code{slimgb} is used.
@*       If met <0, the B-Sigma ideal (cf. Castro and Ucha,
@*       'On the computation of Bernstein-Sato ideals', 2005) is computed.
@*       If 0 < met < P, then the ideal B_P (cf. the paper) is computed.
@*       Otherwise, and by default, the ideal B (cf. the paper) is computed.
@*       If us<>0, then syzygies-driven method is used.
@*       If the output ideal happens to be principal, the list of factors
@*       with their multiplicities is returned instead of the ideal.
@*       If printlevel=1, progress debug messages will be printed,
@*       if printlevel>=2, all the debug messages will be printed.
EXAMPLE: example annfsBMI; shows examples
"
{
  int eng = 0;
  int met = 0;
  int usesyz = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
    if ( size(#)>1 )
    {
      if ( typeof(#[2]) == "int" )
      {
        met = int(#[2]);
      }
    }
    if ( size(#)>2 )
    {
      if ( typeof(#[3]) == "int" )
      {
        usesyz = int(#[3]);
      }
    }

  }
  // returns a list with a ring and an ideal LD in it
  int ppl = printlevel-voice+2;
  //  printf("plevel :%s, voice: %s",printlevel,voice);
  def save = basering;
  int N = nvars(basering);
  //if F has some generators which are zero, int P = ncols(I);
  int P = size(F);
  if (P < ncols(F))
  {
    F = simplify(F,2);
  }
  P = size(F);
  if (P == 0)
  {
    ERROR("zero ideal in the input");
  }
  int Nnew = 2*N+2*P;
  int i,j;
  string s;
  list RL = ringlist(basering);
  list L, Lord;
  list tmp;
  intvec iv;
  L[1] = RL[1];  //char
  L[4] = RL[4];  //char, minpoly
  // check whether vars have admissible names
  list Name  = RL[2];
  list RName;
  for (j=1; j<=P; j++)
  {
    RName[j] = "t("+string(j)+")";
    RName[j+P] = "s("+string(j)+")";
  }
  for(i=1; i<=N; i++)
  {
    for(j=1; j<=size(RName); j++)
    {
      if (Name[i] == RName[j])
      { ERROR("Variable names should not include t(i),s(i)"); }
    }
  }
  // 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 (lp(P),dp);
  tmp[1] = "lp";  //string
  s      = "iv=";
  for (i=1; i<=2*P; i++)
  {
    s = s+"1,";
  }
  s[size(s)]= ";";
  execute(s);
  tmp[2] = iv;  //intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1] = "dp";  //string
  s      = "iv=";
  for (i=1; i<=Nnew; i++)  //actually i<=2*N
  {
    s = s+"1,";
  }
  s[size(s)]= ";";
  execute(s);
  kill s;
  tmp[2]  = iv;
  Lord[2] = tmp;
  tmp[1]  = "C";
  iv      = 0;
  tmp[2]  = iv;
  Lord[3] = tmp;
  tmp     = 0;
  L[3]    = Lord;
  // we are done with the list
  def @R@ = ring(L);
  setring @R@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=P; i++)
  {
    @D[i,i+P] = t(i);
  }
  for(i=1; i<=N; i++)
  {
    @D[2*P+i,2*P+N+i] = 1;
  }
  // L[5] = matrix(UpOneMatrix(Nnew));
  // L[6] = @D;
  def @R = nc_algebra(1,@D);
  setring @R;
  kill @R@;
  dbprint(ppl,"// -1-1- the ring @R(_t,_s,_x,_Dx) is ready");
  dbprint(ppl-1, @R);
  // create the ideal I
  ideal  F = imap(save,F);
  ideal I = t(1)*F[1]+s(1);
  for (j=2; j<=P; j++)
  {
    I = I, t(j)*F[j]+s(j);
  }
  poly p,q;
  for (i=1; i<=N; i++)
  {
    p=0;
    for (j=1; j<=P; j++)
    {
      q = t(j);
      q = diff(F[j],var(2*P+i))*q;
      p = p + q;
    }
    I = I, var(2*P+N+i) + p;
  }
  // -------- the ideal I is ready ----------
  if (usesyz)
  {
    dbprint(ppl,"// -1-1-1 adding syzygy-driven elements to the ideal");
    // -- form the extended Jacobian matrix; do the comps over K[x]
    setring save;
    module JC = module(transpose(jacob(F))); // NxP
    for (j=1; j<=P; j++)
    {
      JC[j] = JC[j] + F[j]*gen(N+j);
    }
    //    matrix JAC[N+P][P];
    dbprint(ppl,"// -1-1-2 computing syzygies of the extended Jacobian matrix over K[_x]");
    dbprint(ppl-1, matrix(JC));
    matrix SJ = transpose(syz(transpose(JC))); //or of std(syz)?
    dbprint(ppl,"// -1-1-3 finished computing syzygies of the ext Jacobian matrix over K[_x]");
    dbprint(ppl-1, matrix(SJ));
    setring @R;
    // add generators: first N comps d_i, then P comps s_j
    matrix SJ = imap(save, SJ);
    poly pi;
    for (i=1; i<=nrows(SJ); i++)
    {
      pi = 0;
      for (j=1; j<=N; j++)
      {
        pi = pi + SJ[i,j]*var(2*P+N+j);
      }
      for (j=1; j<=P; j++)
      {
        pi = pi + SJ[i,N+j]*s(j);
      }
      dbprint(ppl-1, "// -1-1-4 adding element:" + string(pi));
      I = I, pi;
    }
  }
  dbprint(ppl,"// -1-2- starting the elimination of "+string(t(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 t(i) are eliminated");
  dbprint(ppl-1, K);  //K is without t(i)
  // ----------- the ring @R2 ------------
  // _x, _Dx,s;  elim.ord for _x,_Dx.
  // keep: N, i,j,s, tmp, RL
  setring save;
  Nnew = 2*N+P;
  kill Lord, tmp, iv, RName;
  list Lord, tmp;
  intvec iv;
  L[1] = RL[1];  //char
  L[4] = RL[4];  //char, minpoly
  // check whether vars hava admissible names -> done earlier
  // now, create the names for new var
  for (j=1; j<=P; j++)
  {
    tmp[j] = "s("+string(j)+")";
  }
  // DName is defined earlier
  list NName = Name + DName + tmp;
  L[2] = NName;
  tmp = 0;
  // block ord (dp(N),dp);
  string s = "iv=";
  for (i=1; i<=Nnew-P; i++)
  {
    s = s+"1,";
  }
  s[size(s)]=";";
  execute(s);
  tmp[1] = "dp";  //string
  tmp[2] = iv;    //intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1] = "dp";  //string
  s[size(s)] = ",";
  for (j=1; j<=P; j++)
  {
    s = s+"1,";
  }
  s[size(s)]=";";
  execute(s);
  kill s;
  kill NName;
  tmp[2]  = iv;
  Lord[2] = tmp;
  tmp[1]  = "C";
  iv      = 0;
  tmp[2]  = iv;
  Lord[3] = tmp;
  tmp     = 0;
  L[3]    = Lord;
  // we are done with the list. Now add a Plural part
  def @R2@ = ring(L);
  setring @R2@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R2 = nc_algebra(1,@D);
  setring @R2;
  kill @R2@;
  dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,_s) is ready");
  dbprint(ppl-1, @R2);
//  ideal MM = maxideal(1);
//  MM = 0,s,MM;
//  map R01 = @R, MM;
//  ideal K = R01(K);
  ideal F = imap(save,F);  // maybe ideal F = R01(I); ?
  ideal K = imap(@R,K);    // maybe ideal K = R01(I); ?

  if (met <0)
  {
  //K = K,F;     // to compute Bsigma (see "On the computation of Bernstein-Sato ideals"; Castro, Ucha)
    K = K,F;
    dbprint(ppl,"// -2-1-1 computing the ideal B-Sigma from Castro-Ucha");
  }
  if ( (met>0) && (met<=ncols(F) ) )
  {
  /* compute the ideal B_met */
  //j=2;         // for example
  //K = K,F[j];  // to compute Bj (see "On the computation of Bernstein-Sato ideals"; Castro, Ucha)
    dbprint(ppl,"// -2-1-1 computing the ideal B_i from Castro-Ucha");
    K = K, F[met];
  }
  if ( ( met == 0) || (met > ncols(F) ) )
  {
    // that is met ==0 or met> ncols(F)
    /* compute the ideal B */
    poly f=1;
    for (j=1; j<=P; j++)
    {
      f = f*F[j];
    }
    K = K,f;       // to compute B (Bernstein-Sato ideal)
    dbprint(ppl,"// -2-1-1 computing the ideal B from Castro-Ucha");
  }
  dbprint(ppl,"// -2-2- starting the elimination of _x,_Dx in @R2");
  dbprint(ppl-1, K);
  ideal M = engine(K,eng);
  ideal K2 = nselect(M,1..Nnew-P);
  kill K,M;
  dbprint(ppl,"// -2-3- _x,_Dx are eliminated in @R2");
  dbprint(ppl-1, K2);
  // the ring @R3 and factorize
  ring @R3 = 0,s(1..P),dp;
  dbprint(ppl,"// -3-1- the ring @R3(_s) is ready");
  ideal K3 = imap(@R2,K2);
  if (ncols(K3)==1)
  {
    poly p = K3[1];
    dbprint(ppl,"// -3-2- the Bernstein ideal is principal");
    dbprint(ppl,"// -3-2-1- factorization");
    // Warning: now P is an integer
    list Q = factorize(p);         //with constants and multiplicities
    ideal bs; intvec m;
    for (i=2; i<=size(Q[1]); i++)  //we delete Q[1][1] and Q[2][1]
    {
      bs[i-1] = Q[1][i];
      m[i-1]  = Q[2][i];
    }
    //  "--------- Q-ideal factorizes into ---------";  list(bs,m);
    list BS = bs,m;
  }
  else
  {
    // conjecture for some ideals: the Bernstein ideal is principal
    dbprint(ppl,"// -3-2- the Bernstein ideal is not principal");
    ideal BS = K3;
  }
  // create the ring @R4(_x,_Dx,_s) and put the result into it,
  // _x, _Dx,s;  ord "dp".
  // keep: N, i,j,s, tmp, RL
  setring save;
  Nnew = 2*N+P;
  // list RL = ringlist(save);  //is defined earlier
  kill Lord, tmp, iv;
  L = 0;
  list Lord, tmp;
  intvec iv;
  L[1] = RL[1];  //char
  L[4] = RL[4];  //char, minpoly
  // check whether vars hava admissible names -> done earlier
  // now, create the names for new var
  for (j=1; j<=P; j++)
  {
    tmp[j] = "s("+string(j)+")";
  }
  // DName is defined earlier
  list NName = Name + DName + tmp;
  L[2] = NName;
  tmp = 0;
  // dp ordering;
  string s = "iv=";
  for (i=1; i<=Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)]=";";
  execute(s);
  kill s;
  kill NName;
  tmp[1] = "dp";  //string
  tmp[2] = iv;    //intvec
  Lord[1] = tmp;
  tmp[1]  = "C";
  iv      = 0;
  tmp[2]  = iv;
  Lord[2] = tmp;
  tmp     = 0;
  L[3]    = Lord;
  // we are done with the list. Now add a Plural part
  def @R4@ = ring(L);
  setring @R4@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R4 = nc_algebra(1,@D);
  setring @R4;
  kill @R4@;
  dbprint(ppl,"// -4-1- the ring @R4i(_x,_Dx,_s) is ready");
  dbprint(ppl-1, @R4);
  ideal K4 = imap(@R,K);
  intvec saveopt=option(get);
  option(redSB);
  dbprint(ppl,"// -4-2- the final cosmetic std");
  K4 = engine(K4,eng);  //std does the job too
  // total cleanup
  kill @R;
  kill @R2;
  def BS = imap(@R3,BS);
  export BS;
  kill @R3;
  ideal LD = K4;
  export LD;
  option(set,saveopt);
  return(@R4);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),Dp;
  ideal F = x,y,x+y;
  printlevel = 0;
  // *1* let us compute the B ideal
  def A = annfsBMI(F);    setring A;
  LD; // annihilator
  BS; // Bernstein-Sato ideal
  // *2* now, let us compute B-Sigma ideal
  setring r;
  def Sigma = annfsBMI(F,0,-1); setring Sigma;
  print(matrix(lead(LD))); // compact form of leading
  //  monomials from the annihilator
  BS; // Bernstein-Sato B-Sigma ideal: it is principal,
  // so factors and multiplicities are returned
  // *3* and now, let us compute B-i ideal
  setring r;
  def Bi = annfsBMI(F,0,3); // that is F[3]=x+y is taken
  setring Bi;
  print(matrix(lead(LD)));   // compact form of leading
  //  monomials from the annihilator
  BS; // the B_3 ideal: it is principal, so factors
  // and multiplicities are returned
}

proc annfsOT(poly F, list #)
"USAGE:  annfsOT(f [,eng]);  f a poly, eng an optional int
RETURN:  ring
PURPOSE: compute the D-module structure of basering[1/f]*f^s,
@* according to the algorithm by Oaku and Takayama
NOTE:    activate the output ring with the @code{setring} command. In this ring,
@*       - the ideal LD (which is a Groebner basis) is the needed D-module structure,
@*         which is obtained by substituting the minimal integer root of a Bernstein
@*         polynomial into the s-parametric ideal;
@*       - the list BS contains roots with multiplicities of a Bernstein polynomial of f.
@*       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 annfsOT; shows examples
"
{
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  // returns a list with a ring and an ideal LD in it
  int ppl = printlevel-voice+2;
  //  printf("plevel :%s, voice: %s",printlevel,voice);
  def save = basering;
  int N = nvars(basering);
  int Nnew = 2*(N+2);
  int i,j;
  string s;
  list RL = ringlist(basering);
  list L, Lord;
  list tmp;
  intvec iv;
  L[1] = RL[1]; // char
  L[4] = RL[4]; // char, minpoly
  // check whether vars have admissible names
  list Name  = RL[2];
  list RName;
  RName[1] = "u";
  RName[2] = "v";
  RName[3] = "t";
  RName[4] = "Dt";
  for(i=1;i<=N;i++)
  {
    for(j=1; j<=size(RName);j++)
    {
      if (Name[i] == RName[j])
      {
        ERROR("Variable names should not include u,v,t,Dt");
      }
    }
  }
  // now, create the names for new vars
  tmp[1]     = "u";
  tmp[2]     = "v";
  list UName = tmp;
  list DName;
  for(i=1;i<=N;i++)
  {
    DName[i] = "D"+Name[i]; // concat
  }
  tmp    =  0;
  tmp[1] = "t";
  tmp[2] = "Dt";
  list NName = UName +  tmp + Name + DName;
  L[2]   = NName;
  tmp    = 0;
  // Name, Dname will be used further
  kill UName;
  kill NName;
  // block ord (a(1,1),dp);
  tmp[1]  = "a"; // string
  iv      = 1,1;
  tmp[2]  = iv; //intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1]  = "dp"; // string
  s       = "iv=";
  for(i=1;i<=Nnew;i++)
  {
    s = s+"1,";
  }
  s[size(s)]= ";";
  execute(s);
  tmp[2]    = iv;
  Lord[2]   = tmp;
  tmp[1]    = "C";
  iv        = 0;
  tmp[2]    = iv;
  Lord[3]   = tmp;
  tmp       = 0;
  L[3]      = Lord;
  // we are done with the list
  def @R@ = ring(L);
  setring @R@;
  matrix @D[Nnew][Nnew];
  @D[3,4]=1;
  for(i=1; i<=N; i++)
  {
    @D[4+i,N+4+i]=1;
  }
  //  @D[N+3,2*(N+2)]=1; old t,Dt stuff
  //  L[5] = matrix(UpOneMatrix(Nnew));
  //  L[6] = @D;
  def @R = nc_algebra(1,@D);
  setring @R;
  kill @R@;
  dbprint(ppl,"// -1-1- the ring @R(u,v,t,Dt,_x,_Dx) is ready");
  dbprint(ppl-1, @R);
  // create the ideal I
  poly  F = imap(save,F);
  ideal I = u*F-t,u*v-1;
  poly p;
  for(i=1; i<=N; i++)
  {
    p = u*Dt; // u*Dt
    p = diff(F,var(4+i))*p;
    I = I, var(N+4+i) + p;
  }
  // -------- the ideal I is ready ----------
  dbprint(ppl,"// -1-2- starting the elimination of u,v in @R");
  dbprint(ppl-1, I);
  ideal J = engine(I,eng);
  ideal K = nselect(J,1..2);
  dbprint(ppl,"// -1-3- u,v are eliminated");
  dbprint(ppl-1, K);  // K is without u,v
  setring save;
  // ------------ new ring @R2 ------------------
  // without u,v and with the elim.ord for t,Dt
  // tensored with the K[s]
  // keep: N, i,j,s, tmp, RL
  Nnew = 2*N+2+1;
  //  list RL = ringlist(save); // is defined earlier
  L = 0;  //  kill L;
  kill Lord, tmp, iv, RName;
  list Lord, tmp;
  intvec iv;
  L[1] = RL[1]; L[4] = RL[4]; // char, minpoly
  // check whether vars have admissible names -> done earlier
  //  list Name  = RL[2];
  list RName;
  RName[1] = "t";
  RName[2] = "Dt";
  // now, create the names for new var (here, s only)
  tmp[1]     = "s";
  // DName is defined earlier
  list NName = RName + Name + DName + tmp;
  L[2]   = NName;
  tmp    = 0;
  // block ord (a(1,1),dp);
  tmp[1]  = "a";  iv = 1,1; tmp[2]  = iv; //intvec
  Lord[1] = tmp;
  // continue with a(1,1,1,1)...
  tmp[1]  = "dp"; s  = "iv=";
  for(i=1; i<= Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)]= ";";  execute(s);
  kill NName;
  tmp[2]    = iv;
  Lord[2]   = tmp;
  // extra block for s
  // tmp[1] = "dp"; iv = 1;
  //  s[size(s)]= ",";  s = s + "1,1,1;";  execute(s);  tmp[2]    = iv;
  //  Lord[3]   = tmp;
  kill s;
  tmp[1]    = "C";   iv  = 0; tmp[2] = iv;
  Lord[3]   = tmp;   tmp = 0;
  L[3]      = Lord;
  // we are done with the list. Now, add a Plural part
  def @R2@ = ring(L);
  setring @R2@;
  matrix @D[Nnew][Nnew];
  @D[1,2] = 1;
  for(i=1; i<=N; i++)
  {
    @D[2+i,2+N+i] = 1;
  }
  def @R2 = nc_algebra(1,@D);
  setring @R2;
  kill @R2@;
  dbprint(ppl,"// -2-1- the ring @R2(t,Dt,_x,_Dx,s) is ready");
  dbprint(ppl-1, @R2);
  ideal MM = maxideal(1);
  MM = 0,0,MM;
  map R01 = @R, MM;
  ideal K = R01(K);
  //  ideal K = imap(@R,K); // names of vars are important!
  poly G = t*Dt+s+1; // s is a variable here
  K = NF(K,std(G)),G;
  // -------- the ideal K_(@R2) is ready ----------
  dbprint(ppl,"// -2-2- starting the elimination of t,Dt in @R2");
  dbprint(ppl-1, K);
  ideal M  = engine(K,eng);
  ideal K2 = nselect(M,1..2);
  dbprint(ppl,"// -2-3- t,Dt are eliminated");
  dbprint(ppl-1, K2);
  //  dbprint(ppl-1+1," -2-4- std of K2");
  //  option(redSB);  option(redTail);  K2 = std(K2);
  //  K2; // without t,Dt, and with s
  // -------- the ring @R3 ----------
  // _x, _Dx, s; elim.ord for _x,_Dx.
  // keep: N, i,j,s, tmp, RL
  setring save;
  Nnew = 2*N+1;
  //  list RL = ringlist(save); // is defined earlier
  //  kill L;
  kill Lord, tmp, iv, RName;
  list Lord, tmp;
  intvec iv;
  L[1] = RL[1]; L[4] = RL[4]; // char, minpoly
  // check whether vars have admissible names -> done earlier
  //  list Name  = RL[2];
  // now, create the names for new var (here, s only)
  tmp[1]     = "s";
  // DName is defined earlier
  list NName = Name + DName + tmp;
  L[2]   = NName;
  tmp    = 0;
  // block ord (a(1,1...),dp);
  string  s = "iv=";
  for(i=1; i<=Nnew-1; i++)
  {
    s = s+"1,";
  }
  s[size(s)]= ";";
  execute(s);
  tmp[1]  = "a"; // string
  tmp[2]  = iv; //intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1]  = "dp"; // string
  s[size(s)]=","; s= s+"1;";
  execute(s);
  kill s;
  kill NName;
  tmp[2]    = iv;
  Lord[2]   = tmp;
  tmp[1]    = "C";  iv  = 0;  tmp[2] = iv;
  Lord[3]   = tmp;  tmp = 0;
  L[3]      = Lord;
  // we are done with the list. Now add a Plural part
  def @R3@ = ring(L);
  setring @R3@;
  matrix @D[Nnew][Nnew];
  for(i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R3 = nc_algebra(1,@D);
  setring @R3;
  kill @R3@;
  dbprint(ppl,"// -3-1- the ring @R3(_x,_Dx,s) is ready");
  dbprint(ppl-1, @R3);
  ideal MM = maxideal(1);
  MM = 0,0,MM;
  map R12 = @R2, MM;
  ideal K = R12(K2);
  poly  F = imap(save,F);
  K = K,F;
  dbprint(ppl,"// -3-2- starting the elimination of _x,_Dx in @R3");
  dbprint(ppl-1, K);
  ideal M = engine(K,eng);
  ideal K3 = nselect(M,1..Nnew-1);
  dbprint(ppl,"// -3-3-  _x,_Dx are eliminated in @R3");
  dbprint(ppl-1, K3);
  // the ring @R4  and the search for minimal negative int s
  ring @R4 = 0,(s),dp;
  dbprint(ppl,"// -4-1- the ring @R4 is ready");
  ideal K4 = imap(@R3,K3);
  poly p = K4[1];
  dbprint(ppl,"// -4-2- factorization");
////  ideal P = factorize(p,1);  // without constants and multiplicities
  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];
  }
  //  "------ b-function factorizes into ----------";  P;
////  int sP  = minIntRoot(P, 1);
  int sP = minIntRoot(bs,1);
  dbprint(ppl,"// -4-3- minimal integer root found");
  dbprint(ppl-1, sP);
  // convert factors to a list of their roots
  // assume all factors are linear
////  ideal BS = normalize(P);
////  BS = subst(BS,s,0);
////  BS = -BS;
  bs = normalize(bs);
  bs = subst(bs,s,0);
  bs = -bs;
  list BS = bs,m;
  // TODO: sort BS!
  // ------ substitute s found in the ideal ------
  // ------- going back to @R2 and substitute --------
  setring @R2;
  ideal K3 = subst(K2,s,sP);
  // create the ordinary Weyl algebra and put the result into it,
  // thus creating the ring @R5
  // keep: N, i,j,s, tmp, RL
  setring save;
  Nnew = 2*N;
  //  list RL = ringlist(save); // is defined earlier
  kill Lord, tmp, iv;
  L = 0;
  list Lord, tmp;
  intvec iv;
  L[1] = RL[1];   L[4] = RL[4]; // char, minpoly
  // check whether vars have admissible names -> done earlier
  //  list Name  = RL[2];
  // DName is defined earlier
  list NName = Name + DName;
  L[2]   = NName;
  // dp ordering;
  string   s       = "iv=";
  for(i=1;i<=Nnew;i++)
  {
    s = s+"1,";
  }
  s[size(s)]= ";";
  execute(s);
  tmp     = 0;
  tmp[1]  = "dp"; // string
  tmp[2]  = iv; //intvec
  Lord[1] = tmp;
  kill s;
  tmp[1]    = "C";
  iv        = 0;
  tmp[2]    = iv;
  Lord[2]   = tmp;
  tmp       = 0;
  L[3]      = Lord;
  // we are done with the list
  // Add: Plural part
  def @R5@ = ring(L);
  setring @R5@;
  matrix @D[Nnew][Nnew];
  for(i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R5 = nc_algebra(1,@D);
  setring @R5;
  kill @R5@;
  dbprint(ppl,"// -5-1- the ring @R5 is ready");
  dbprint(ppl-1, @R5);
  ideal K5 = imap(@R2,K3);
  intvec saveopt=option(get);
  option(redSB);
  dbprint(ppl,"// -5-2- the final cosmetic std");
  K5 = engine(K5,eng); // std does the job too
  // total cleanup
  kill @R;
  kill @R2;
  kill @R3;
////  ideal BS = imap(@R4,BS);
  list BS = imap(@R4,BS);
  export BS;
  ideal LD = K5;
  kill @R4;
  export LD;
  option(set,saveopt);
  return(@R5);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),Dp;
  poly F = x^2+y^3+z^5;
  printlevel = 0;
  def A  = annfsOT(F);
  setring A;
  LD;
  BS;
}


proc SannfsOT(poly F, list #)
"USAGE:  SannfsOT(f [,eng]);  f a poly, eng an optional int
RETURN:  ring
PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the
@* 1st step of the algorithm by Oaku and Takayama in the ring D[s]
NOTE:    activate the output ring with the @code{setring} command.
@*  In the output ring D[s], the ideal LD (which is NOT a Groebner basis)
@*  is the needed D-module structure.
@*  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 SannfsOT; shows examples
"
{
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  // returns a list with a ring and an ideal LD in it
  int ppl = printlevel-voice+2;
  //  printf("plevel :%s, voice: %s",printlevel,voice);
  def save = basering;
  int N = nvars(basering);
  int Nnew = 2*(N+2);
  int i,j;
  string s;
  list RL = ringlist(basering);
  list L, Lord;
  list tmp;
  intvec iv;
  L[1] = RL[1]; // char
  L[4] = RL[4]; // char, minpoly
  // check whether vars have admissible names
  list Name  = RL[2];
  list RName;
  RName[1] = "u";
  RName[2] = "v";
  RName[3] = "t";
  RName[4] = "Dt";
  for(i=1;i<=N;i++)
  {
    for(j=1; j<=size(RName);j++)
    {
      if (Name[i] == RName[j])
      {
        ERROR("Variable names should not include u,v,t,Dt");
      }
    }
  }
  // now, create the names for new vars
  tmp[1]     = "u";
  tmp[2]     = "v";
  list UName = tmp;
  list DName;
  for(i=1;i<=N;i++)
  {
    DName[i] = "D"+Name[i]; // concat
  }
  tmp    =  0;
  tmp[1] = "t";
  tmp[2] = "Dt";
  list NName = UName +  tmp + Name + DName;
  L[2]   = NName;
  tmp    = 0;
  // Name, Dname will be used further
  kill UName;
  kill NName;
  // block ord (a(1,1),dp);
  tmp[1]  = "a"; // string
  iv      = 1,1;
  tmp[2]  = iv; //intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1]  = "dp"; // string
  s       = "iv=";
  for(i=1;i<=Nnew;i++)
  {
    s = s+"1,";
  }
  s[size(s)]= ";";
  execute(s);
  tmp[2]    = iv;
  Lord[2]   = tmp;
  tmp[1]    = "C";
  iv        = 0;
  tmp[2]    = iv;
  Lord[3]   = tmp;
  tmp       = 0;
  L[3]      = Lord;
  // we are done with the list
  def @R@ = ring(L);
  setring @R@;
  matrix @D[Nnew][Nnew];
  @D[3,4]=1;
  for(i=1; i<=N; i++)
  {
    @D[4+i,N+4+i]=1;
  }
  //  @D[N+3,2*(N+2)]=1; old t,Dt stuff
  //  L[5] = matrix(UpOneMatrix(Nnew));
  //  L[6] = @D;
  def @R =  nc_algebra(1,@D);
  setring @R;
  kill @R@;
  dbprint(ppl,"// -1-1- the ring @R(u,v,t,Dt,_x,_Dx) is ready");
  dbprint(ppl-1, @R);
  // create the ideal I
  poly  F = imap(save,F);
  ideal I = u*F-t,u*v-1;
  poly p;
  for(i=1; i<=N; i++)
  {
    p = u*Dt; // u*Dt
    p = diff(F,var(4+i))*p;
    I = I, var(N+4+i) + p;
  }
  // -------- the ideal I is ready ----------
  dbprint(ppl,"// -1-2- starting the elimination of u,v in @R");
  dbprint(ppl-1, I);
  ideal J = engine(I,eng);
  ideal K = nselect(J,1..2);
  dbprint(ppl,"// -1-3- u,v are eliminated");
  dbprint(ppl-1, K);  // K is without u,v
  setring save;
  // ------------ new ring @R2 ------------------
  // without u,v and with the elim.ord for t,Dt
  // tensored with the K[s]
  // keep: N, i,j,s, tmp, RL
  Nnew = 2*N+2+1;
  //  list RL = ringlist(save); // is defined earlier
  L = 0;  //  kill L;
  kill Lord, tmp, iv, RName;
  list Lord, tmp;
  intvec iv;
  L[1] = RL[1]; L[4] = RL[4]; // char, minpoly
  // check whether vars have admissible names -> done earlier
  //  list Name  = RL[2];
  list RName;
  RName[1] = "t";
  RName[2] = "Dt";
  // now, create the names for new var (here, s only)
  tmp[1]     = "s";
  // DName is defined earlier
  list NName = RName + Name + DName + tmp;
  L[2]   = NName;
  tmp    = 0;
  // block ord (a(1,1),dp);
  tmp[1]  = "a";  iv = 1,1; tmp[2]  = iv; //intvec
  Lord[1] = tmp;
  // continue with a(1,1,1,1)...
  tmp[1]  = "dp"; s  = "iv=";
  for(i=1; i<= Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)]= ";";  execute(s);
  kill NName;
  tmp[2]    = iv;
  Lord[2]   = tmp;
  // extra block for s
  // tmp[1] = "dp"; iv = 1;
  //  s[size(s)]= ",";  s = s + "1,1,1;";  execute(s);  tmp[2]    = iv;
  //  Lord[3]   = tmp;
  kill s;
  tmp[1]    = "C";   iv  = 0; tmp[2] = iv;
  Lord[3]   = tmp;   tmp = 0;
  L[3]      = Lord;
  // we are done with the list. Now, add a Plural part
  def @R2@ = ring(L);
  setring @R2@;
  matrix @D[Nnew][Nnew];
  @D[1,2] = 1;
  for(i=1; i<=N; i++)
  {
    @D[2+i,2+N+i] = 1;
  }
  def @R2 = nc_algebra(1,@D);
  setring @R2;
  kill @R2@;
  dbprint(ppl,"// -2-1- the ring @R2(t,Dt,_x,_Dx,s) is ready");
  dbprint(ppl-1, @R2);
  ideal MM = maxideal(1);
  MM = 0,0,MM;
  map R01 = @R, MM;
  ideal K = R01(K);
  //  ideal K = imap(@R,K); // names of vars are important!
  poly G = t*Dt+s+1; // s is a variable here
  K = NF(K,std(G)),G;
  // -------- the ideal K_(@R2) is ready ----------
  dbprint(ppl,"// -2-2- starting the elimination of t,Dt in @R2");
  dbprint(ppl-1, K);
  ideal M  = engine(K,eng);
  ideal K2 = nselect(M,1..2);
  dbprint(ppl,"// -2-3- t,Dt are eliminated");
  dbprint(ppl-1, K2);
  K2 = engine(K2,eng);
  setring save;
  // ----------- the ring @R3 ------------
  // _x, _Dx,s;  elim.ord for _x,_Dx.
  // keep: N, i,j,s, tmp, RL
  Nnew = 2*N+1;
  kill Lord, tmp, iv, RName;
  list Lord, tmp;
  intvec iv;
  L[1] = RL[1];
  L[4] = RL[4];  // char, minpoly
  // check whether vars hava admissible names -> done earlier
  // now, create the names for new var
  tmp[1] = "s";
  // DName is defined earlier
  list NName = Name + DName + tmp;
  L[2] = NName;
  tmp = 0;
  // block ord (dp(N),dp);
  string s = "iv=";
  for (i=1; i<=Nnew-1; i++)
  {
    s = s+"1,";
  }
  s[size(s)]=";";
  execute(s);
  tmp[1] = "dp";  // string
  tmp[2] = iv;   // intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1] = "dp";  // string
  s[size(s)] = ",";
  s = s+"1;";
  execute(s);
  kill s;
  kill NName;
  tmp[2]      = iv;
  Lord[2]     = tmp;
  tmp[1]      = "C";  iv  = 0;  tmp[2]=iv;
  Lord[3]     = tmp;  tmp = 0;
  L[3]        = Lord;
  // we are done with the list. Now add a Plural part
  def @R3@ = ring(L);
  setring @R3@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R3 = nc_algebra(1,@D);
  setring @R3;
  kill @R3@;
  dbprint(ppl,"// -3-1- the ring @R3(_x,_Dx,s) is ready");
  dbprint(ppl-1, @R3);
  ideal MM = maxideal(1);
  MM = 0,s,MM;
  map R01 = @R2, MM;
  ideal K2 = R01(K2);
  // total cleanup
  ideal LD = K2;
  // make leadcoeffs positive
  for (i=1; i<= ncols(LD); i++)
  {
    if (leadcoef(LD[i]) <0 )
    {
      LD[i] = -LD[i];
    }
  }
  export LD;
  kill @R;
  kill @R2;
  return(@R3);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),Dp;
  poly F = x^3+y^3+z^3;
  printlevel = 0;
  def A  = SannfsOT(F);
  setring A;
  LD;
}

proc SannfsBM(poly F, list #)
"USAGE:  SannfsBM(f [,eng]);  f a poly, eng an optional int
RETURN:  ring
PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the
@* 1st step of the algorithm by Briancon and Maisonobe in the ring D[s].
NOTE:  activate the output ring with the @code{setring} command.
@*   In the output ring D[s], the ideal LD (which is NOT a Groebner basis) is
@*   the needed D-module structure.
@*   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 SannfsBM; shows examples
"
{
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  // returns a list with a ring and an ideal LD in it
  int ppl = printlevel-voice+2;
  //  printf("plevel :%s, voice: %s",printlevel,voice);
  def save = basering;
  int N = nvars(basering);
  int Nnew = 2*N+2;
  int i,j;
  string s;
  list RL = ringlist(basering);
  list L, Lord;
  list tmp;
  intvec iv;
  L[1] = RL[1]; // char
  L[4] = RL[4]; // char, minpoly
  // check whether vars have admissible names
  list Name  = RL[2];
  list RName;
  RName[1] = "t";
  RName[2] = "s";
  for(i=1;i<=N;i++)
  {
    for(j=1; j<=size(RName);j++)
    {
      if (Name[i] == RName[j])
      {
        ERROR("Variable names should not include t,s");
      }
    }
  }
  // now, create the names for new vars
  list DName;
  for(i=1;i<=N;i++)
  {
    DName[i] = "D"+Name[i]; // concat
  }
  tmp[1] = "t";
  tmp[2] = "s";
  list NName = tmp + Name + DName;
  L[2]   = NName;
  // Name, Dname will be used further
  kill NName;
  // block ord (lp(2),dp);
  tmp[1]  = "lp"; // string
  iv      = 1,1;
  tmp[2]  = iv; //intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1]  = "dp"; // string
  s       = "iv=";
  for(i=1;i<=Nnew;i++)
  {
    s = s+"1,";
  }
  s[size(s)]= ";";
  execute(s);
  kill s;
  tmp[2]    = iv;
  Lord[2]   = tmp;
  tmp[1]    = "C";
  iv        = 0;
  tmp[2]    = iv;
  Lord[3]   = tmp;
  tmp       = 0;
  L[3]      = Lord;
  // we are done with the list
  def @R@ = ring(L);
  setring @R@;
  matrix @D[Nnew][Nnew];
  @D[1,2]=t;
  for(i=1; i<=N; i++)
  {
    @D[2+i,N+2+i]=1;
  }
  //  L[5] = matrix(UpOneMatrix(Nnew));
  //  L[6] = @D;
  def @R = nc_algebra(1,@D);
  setring @R;
  kill @R@;
  dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready");
  dbprint(ppl-1, @R);
  // create the ideal I
  poly  F = imap(save,F);
  ideal I = t*F+s;
  poly p;
  for(i=1; i<=N; i++)
  {
    p = t; // t
    p = diff(F,var(2+i))*p;
    I = I, var(N+2+i) + p;
  }
  // -------- the ideal I is ready ----------
  dbprint(ppl,"// -1-2- starting the elimination of t in @R");
  dbprint(ppl-1, I);
  ideal J = engine(I,eng);
  ideal K = nselect(J,1);
  dbprint(ppl,"// -1-3- t is eliminated");
  dbprint(ppl-1, K);  // K is without t
  K = engine(K,eng);  // std does the job too
  // now, we must change the ordering
  // and create a ring without t, Dt
  //  setring S;
  // ----------- the ring @R3 ------------
  // _x, _Dx,s;  elim.ord for _x,_Dx.
  // keep: N, i,j,s, tmp, RL
  Nnew = 2*N+1;
  kill Lord, tmp, iv, RName;
  list Lord, tmp;
  intvec iv;
  list L=imap(save,L);
  list RL=imap(save,RL);
  L[1] = RL[1];
  L[4] = RL[4];  // char, minpoly
  // check whether vars hava admissible names -> done earlier
  // now, create the names for new var
  tmp[1] = "s";
  // DName is defined earlier
  list NName = Name + DName + tmp;
  L[2] = NName;
  tmp = 0;
  // block ord (dp(N),dp);
  string s = "iv=";
  for (i=1; i<=Nnew-1; i++)
  {
    s = s+"1,";
  }
  s[size(s)]=";";
  execute(s);
  tmp[1] = "dp";  // string
  tmp[2] = iv;   // intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1] = "dp";  // string
  s[size(s)] = ",";
  s = s+"1;";
  execute(s);
  kill s;
  kill NName;
  tmp[2]      = iv;
  Lord[2]     = tmp;
  tmp[1]      = "C";  iv  = 0;  tmp[2]=iv;
  Lord[3]     = tmp;  tmp = 0;
  L[3]        = Lord;
  // we are done with the list. Now add a Plural part
  def @R2@ = ring(L);
  setring @R2@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R2 = nc_algebra(1,@D);
  setring @R2;
  kill @R2@;
  dbprint(ppl,"//  -2-1- the ring @R2(_x,_Dx,s) is ready");
  dbprint(ppl-1, @R2);
  ideal MM = maxideal(1);
  MM = 0,s,MM;
  map R01 = @R, MM;
  ideal K = R01(K);
  // total cleanup
  ideal LD = K;
  // make leadcoeffs positive
  for (i=1; i<= ncols(LD); i++)
  {
    if (leadcoef(LD[i]) <0 )
    {
      LD[i] = -LD[i];
    }
  }
  export LD;
  kill @R;
  return(@R2);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),Dp;
  poly F = x^3+y^3+z^3;
  printlevel = 0;
  def A = SannfsBM(F);
  setring A;
  LD;
}

static proc safeVarName (string s, list #)
{
  string S;
  int cv = 1;
  if (size(#)>1)
  {
    if (#[1]=="v") { cv = 0; S = varstr(basering); }
    if (#[1]=="c") { cv = 0; S = charstr(basering); }
  }
  if (cv) { S = charstr(basering) + "," + varstr(basering); }
  S = "," + S + ",";
  s = "," + s + ",";
  while (find(S,s) <> 0)
  {
    s[1] = "@";
    s = "," + s;
  }
  s = s[2..size(s)-1];
  return(s)
}

proc SannfsBFCT(poly F, list #)
"USAGE:  SannfsBFCT(f [,a,b,c]);  f a poly, a,b,c optional ints
RETURN:  ring
PURPOSE: compute a Groebner basis either of Ann(f^s)+<f> or of
@*       Ann(f^s)+<f,f_1,...,f_n> in D[s]
NOTE:    Activate the output ring with the @code{setring} command.
@*  This procedure, unlike SannfsBM, returns the ring D[s] with an anti-
@*  elimination ordering for s.
@*  The output ring contains an ideal @code{LD}, being a Groebner basis
@*  either of  Ann(f^s)+<f>, if a=0 (and by default), or of
@*  Ann(f^s)+<f,f_1,...,f_n>, otherwise.
@*  Here, f_i stands for the i-th  partial derivative of f.
@*  If b<>0, @code{std} is used for Groebner basis computations,
@*  otherwise, and by default @code{slimgb} is used.
@*  If c<>0, @code{std} is used for Groebner basis computations of
@*  ideals <I+J> when I is already a Groebner basis of <I>.
@*  Otherwise, and by default the engine determined by the switch b is
@*  used. Note that in the case c<>0, the choice for b will be
@*  overwritten only for the types of ideals mentioned above.
@*  This means that if b<>0, specifying c has no effect.
DISPLAY: If printlevel=1, progress debug messages will be printed,
@*       if printlevel>=2, all the debug messages will be printed.
EXAMPLE: example SannfsBFCT; shows examples
"
{
  int addPD,eng,stdsum;
  if (size(#)>0)
  {
    if (typeof(#[1])=="int" || typeof(#[1])=="number")
    {
      addPD = int(#[1]);
    }
    if (size(#)>1)
    {
      if (typeof(#[2])=="int" || typeof(#[2])=="number")
      {
        eng = int(#[2]);
      }
      if (size(#)>2)
      {
        if (typeof(#[3])=="int" || typeof(#[3])=="number")
        {
          stdsum = int(#[3]);
        }
      }
    }
  }
  int ppl = printlevel-voice+2;
  def save = basering;
  int N = nvars(save);
  intvec optSave = option(get);
  int i,j;
  list RL = ringlist(save);
  // ----- step 1: compute syzigies
  intvec iv;
  list L,Lord;
  iv = 1:N; Lord[1] = list("dp",iv);
  iv = 0;   Lord[2] = list("C",iv);
  L = RL;
  L[3] = Lord;
  def @RM = ring(L);
  kill L,Lord;
  setring @RM;
  option(redSB);
  option(redTail);
  def RM = makeModElimRing(@RM);
  setring RM;
  poly F = imap(save,F);
  ideal J = jacob(F);
  J = F,J;
  dbprint(ppl,"// -1-1- Starting the computation of syz(F,_Dx(F))");
  dbprint(ppl-1, J);
  module M = syz(J);
  dbprint(ppl,"// -1-2- The module syz(F,_Dx(F)) has been computed");
  dbprint(ppl-1, M);
  dbprint(ppl,"// -1-3- Starting GB computation of syz(F,_Dx(F))");
  M = engine(M,eng);
  dbprint(ppl,"// -1-4- GB computation finished");
  dbprint(ppl-1, M);
  // ----- step 2: compute part of Ann(F^s)
  setring save;
  option(set,optSave);
  module M = imap(RM,M);
  kill optSave,RM;
  // ----- create D[s]
  int Nnew = 2*N+1;
  list L, Lord;
  // ----- keep char, minpoly
  L[1] = RL[1];
  L[4] = RL[4];
  // ----- create names for new vars
  list Name  = RL[2];
  string newVar@s = safeVarName("s");
  if (newVar@s[1] == "@")
  {
    print("Name s already assigned to parameter/ringvar.");
    print("Using " + newVar@s + " instead.")
  }
  list DName;
  for (i=1; i<=N; i++)
  {
    DName[i] = safeVarName("D" + Name[i]);
  }
  L[2] = list(newVar@s) + Name + DName;
  // ----- create ordering
  // --- anti-elimination ordering for s
  iv = 1;       Lord[1] = list("dp",iv);
  iv = 1:(2*N); Lord[2] = list("dp",iv);
  iv = 0;       Lord[3] = list("C",iv);
  L[3] = Lord;
  // ----- create commutative ring
  def @Ds = ring(L);
  kill L,Lord;
  setring @Ds;
  // ----- create nc relations
   matrix Drel[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    Drel[i+1,N+1+i] = 1;
  }
  def Ds = nc_algebra(1,Drel);
  setring Ds;
  kill @Ds;
  dbprint(ppl,"// -2-1- The ring D[s] is ready");
  dbprint(ppl-1, Ds);
  matrix M = imap(save,M);
  vector v = var(1)*gen(1);
  for (i=1; i<=N; i++)
  {
    v = v + var(i+1+N)*gen(i+1); //[s,_Dx]
  }
  ideal J = transpose(M)*v;
  kill M,v;
  dbprint(ppl,"// -2-2- Compute part of Ann(F^s)");
  dbprint(ppl-1, J);
  J = engine(J,eng);
  dbprint(ppl,"// -2-3- GB computation finished");
  dbprint(ppl-1, J);
  // ----- step 3: the full annihilator
  // ----- create D<t,s>
  setring save;
  Nnew = 2*N+2;
  list L, Lord;
  // ----- keep char, minpoly
  L[1] = RL[1];
  L[4] = RL[4];
  // ----- create vars
  string newVar@t = safeVarName("t");
  L[2] = list(newVar@t,newVar@s) + DName + Name;
  // ----- create ordering for elimination of t
  // block ord (lp(2),dp);
  iv = 1,1;    Lord[1] = list("lp",iv);
  iv = 1:Nnew; Lord[2] = list("dp",iv);
  iv = 0;      Lord[3] = list("C",iv);
  L[3] = Lord;
  def @Dts = ring(L);
  kill RL,L,Lord,Name,DName,newVar@s,newVar@t;
  setring @Dts;
  // ----- create nc relations
  matrix Drel[Nnew][Nnew];
  Drel[1,2] = var(1);
  for(i=1; i<=N; i++)
  {
    Drel[2+i,N+2+i]=-1;
  }
  def Dts = nc_algebra(1,Drel);
  setring Dts;
  kill @Dts;
  dbprint(ppl,"// -3-1- The ring D<t,s> is ready");
  dbprint(ppl-1, Dts);
  // ----- create the ideal I following BM
  poly  F = imap(save,F);
  ideal I = var(1)*F + var(2); // = t*F + s
  poly p;
  for(i=1; i<=N; i++)
  {
    p = var(1)*diff(F,var(N+2+i)) + var(2+i); // = t*F_i + D_i
    I[i+1] = p;
  }
  // ----- add already computed part to it
  ideal MM = var(2);      // s
  for (i=1; i<=N; i++)
  {
    MM[1+i] = var(2+N+i); // _x
    MM[1+N+i] = var(2+i); // _Dx
  }
  map Ds2Dts = Ds,MM;
  ideal J = Ds2Dts(J);
  attrib(J,"isSB",1);
  kill MM,Ds2Dts;
  // ----- start the elimination
  dbprint(ppl,"// -3-2- Starting the elimination of t in D<t,s>");
  dbprint(ppl-1, I);
  if (stdsum || eng <> 0)
  {
    I = std(J,I);
  }
  else
  {
    I = J,I;
    I = engine(I,eng);
  }
  kill J;
  I = nselect(I,1);
  dbprint(ppl,"// -3-3- t is eliminated");
  dbprint(ppl-1, I);  // I is without t
  // ----- step 4: add F
  // ----- back to D[s]
  setring Ds;
  ideal MM = 0,var(1);      // 0,s
  for (i=1; i<=N; i++)
  {
    MM[2+i]   = var(1+N+i); // _Dx
    MM[2+N+i] = var(1+i);   // _x
  }
  map Dts2Ds = Dts, MM;
  ideal LD = Dts2Ds(I);
  kill J,Dts,Dts2Ds,MM;
  dbprint(ppl,"// -4-1- Starting cosmetic Groebner computation");
  LD = engine(LD,eng);
  dbprint(ppl,"// -4-2- Finished cosmetic Groebner computation");
  dbprint(ppl-1, LD);
  // ----- use reduction trick as Macaulay2 does: compute b(s)/(s+1) by adding all partial derivations also
  ideal J;
  if (addPD)
  {
    setring @RM;
    poly F = imap(save,F);
    ideal J = jacob(F);
    J = F,J;
    dbprint(ppl,"// -4-2-1- Start GB computation <f, f_i>");
    J = engine(J,eng);
    dbprint(ppl,"// -4-2-2- Finished GB computation <f, f_i>");
    dbprint(ppl-1, J);
    setring Ds;
    J = imap(@RM,J);
    attrib(J,"isSB",1);
    dbprint(ppl,"// -4-3- Start GB computations for Ann f^s + <f, f_i>");
  }
  else
  {
    J = imap(save,F);
    dbprint(ppl,"// -4-3- Start GB computations for Ann f^s + <f>");
  }
  kill @RM;
  // ----- the really hard part
  if (stdsum || eng <> 0)
  {
    LD = std(LD,J);
  }
  else
  {
    LD = LD,J;
    LD = engine(LD,eng);
  }
  if (addPD) { dbprint(ppl,"// -4-4- Finished GB computations for Ann f^s + <f, f_i>"); }
  else       { dbprint(ppl,"// -4-4- Finished GB computations for Ann f^s + <f>"); }
  dbprint(ppl-1, LD);
  export LD;
  return(Ds);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z,w),Dp;
  poly F = x^3+y^3+z^3*w;
  // compute Ann(F^s)+<F> using slimgb only
  def A = SannfsBFCT(F);
  setring A; A;
  LD;
  // the Bernstein-Sato poly of F:
  vec2poly(pIntersect(s,LD));
  // a fancier example:
  def R = reiffen(4,5); setring R;
  RC; // the Reiffen curve in 4,5
  // compute Ann(RC^s)+<RC,diff(RC,x),diff(RC,y)>
  // using std for GB computations of ideals <I+J>
  // where I is already a GB of <I>
  // and slimgb for other ideals
  def B = SannfsBFCT(RC,1,0,1);
  setring B;
  // the Bernstein-Sato poly of RC:
  (s-1)*vec2poly(pIntersect(s,LD));
}


proc SannfsBFCTstd(poly F, list #)
"USAGE:  SannfsBFCTstd(f [,eng]);  f a poly, eng an optional int
RETURN:  ring
PURPOSE: compute Ann f^s and Groebner basis of Ann f^s+f in D[s]
NOTE:    activate the output ring with the @code{setring} command.
@*    This procedure, unlike SannfsBM, returns a ring with the degrevlex
@*    ordering in all variables.
@*    In the ring D[s], the ideal LD (which IS a Groebner basis) is the needed ideal.
@*    In this procedure @code{std} is used for Groebner basis computations.
DISPLAY: If printlevel=1, progress debug messages will be printed,
@*       if printlevel>=2, all the debug messages will be printed.
EXAMPLE: example SannfsBFCTstd; shows examples
"
{
  // DEBUG INFO: ordering on the output ring = dp,
  // use std(K,F); for reusing the std property of K

  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  // returns a list with a ring and an ideal LD in it
  int ppl = printlevel-voice+2;
  //  printf("plevel :%s, voice: %s",printlevel,voice);
  def save = basering;
  int N = nvars(basering);
  int Nnew = 2*N+2;
  int i,j;
  string s;
  list RL = ringlist(basering);
  list L, Lord;
  list tmp;
  intvec iv;
  L[1] = RL[1]; // char
  L[4] = RL[4]; // char, minpoly
  // check whether vars have admissible names
  list Name  = RL[2];
  list RName;
  RName[1] = "@t";
  RName[2] = "@s";
  for(i=1;i<=N;i++)
  {
    for(j=1; j<=size(RName);j++)
    {
      if (Name[i] == RName[j])
      {
        ERROR("Variable names should not include @t,@s");
      }
    }
  }
  // now, create the names for new vars
  list DName;
  for(i=1;i<=N;i++)
  {
    DName[i] = "D"+Name[i]; // concat
  }
  tmp[1] = "t";
  tmp[2] = "s";
  list NName = tmp + DName + Name ;
  L[2]   = NName;
  // Name, Dname will be used further
  kill NName;
  // block ord (lp(2),dp);
  tmp[1]  = "lp"; // string
  iv      = 1,1;
  tmp[2]  = iv; //intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1]  = "dp"; // string
  s       = "iv=";
  for(i=1;i<=Nnew;i++)
  {
    s = s+"1,";
  }
  s[size(s)]= ";";
  execute(s);
  kill s;
  tmp[2]    = iv;
  Lord[2]   = tmp;
  tmp[1]    = "C";
  iv        = 0;
  tmp[2]    = iv;
  Lord[3]   = tmp;
  tmp       = 0;
  L[3]      = Lord;
  // we are done with the list
  def @R@ = ring(L);
  setring @R@;
  matrix @D[Nnew][Nnew];
  @D[1,2]=t;
  for(i=1; i<=N; i++)
  {
    @D[2+i,N+2+i]=-1;
  }
  //  L[5] = matrix(UpOneMatrix(Nnew));
  //  L[6] = @D;
  def @R = nc_algebra(1,@D);
  setring @R;
  kill @R@;
  dbprint(ppl,"// -1-1- the ring @R(t,s,_Dx,_x) is ready");
  dbprint(ppl-1, @R);
  // create the ideal I
  poly  F = imap(save,F);
  ideal I = t*F+s;
  poly p;
  for(i=1; i<=N; i++)
  {
    p = t; // t
    p = diff(F,var(N+2+i))*p;
    I = I, var(2+i) + p;
  }
  // -------- the ideal I is ready ----------
  dbprint(ppl,"// -1-2- starting the elimination of t in @R");
  dbprint(ppl-1, I);
  ideal J = engine(I,eng);
  ideal K = nselect(J,1);
  dbprint(ppl,"// -1-3- t is eliminated");
  dbprint(ppl-1, K);  // K is without t
  K = engine(K,eng);  // std does the job too
  // now, we must change the ordering
  // and create a ring without t
  //  setring S;
  // ----------- the ring @R3 ------------
  // _Dx,_x,s;  +fast ord !
  // keep: N, i,j,s, tmp, RL
  Nnew = 2*N+1;
  kill Lord, tmp, iv, RName;
  list Lord, tmp;
  intvec iv;
  list L=imap(save,L);
  list RL=imap(save,RL);
  L[1] = RL[1];
  L[4] = RL[4];  // char, minpoly
  // check whether vars hava admissible names -> done earlier
  // now, create the names for new var
  tmp[1] = "s";
  // DName is defined earlier
  list NName = DName + Name + tmp;
  L[2] = NName;
  tmp = 0;
  // just dp
  // block ord (dp(N),dp);
  string s = "iv=";
  for (i=1; i<=Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)]=";";
  execute(s);
  tmp[1] = "dp";  // string
  tmp[2] = iv;   // intvec
  Lord[1] = tmp;
  kill s;
  kill NName;
  tmp[1]      = "C";
  Lord[2]     = tmp;  tmp = 0;
  L[3]        = Lord;
  // we are done with the list. Now add a Plural part
  def @R2@ = ring(L);
  setring @R2@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=-1;
  }
  def @R2 = nc_algebra(1,@D);
  setring @R2;
  kill @R2@;
  dbprint(ppl,"//  -2-1- the ring @R2(_Dx,_x,s) is ready");
  dbprint(ppl-1, @R2);
  ideal MM = maxideal(1);
  MM = 0,s,MM;
  map R01 = @R, MM;
  ideal K = R01(K);
  // total cleanup
  poly F = imap(save, F);
  //  ideal LD = K,F;
  dbprint(ppl,"//  -2-2- start GB computations for Ann f^s + f");
  //  dbprint(ppl-1, LD);
  ideal LD = std(K,F);
  //  LD = engine(LD,eng);
  dbprint(ppl,"//  -2-3- finished GB computations for Ann f^s + f");
  dbprint(ppl-1, LD);
  // make leadcoeffs positive
  for (i=1; i<= ncols(LD); i++)
  {
    if (leadcoef(LD[i]) <0 )
    {
      LD[i] = -LD[i];
    }
  }
  export LD;
  kill @R;
  return(@R2);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z,w),Dp;
  poly F = x^3+y^3+z^3*w;
  printlevel = 0;
  def A = SannfsBFCT(F); setring A;
  intvec v = 1,2,3,4,5,6,7,8;
  // are there polynomials, depending on s only?
  nselect(LD,v);
  // a fancier example:
  def R = reiffen(4,5); setring R;
  v = 1,2,3,4;
  RC; // the Reiffen curve in 4,5
  def B = SannfsBFCT(RC);
  setring B;
  // Are there polynomials, depending on s only?
  nselect(LD,v);
  // It is not the case. Are there leading monomials in s only?
  nselect(lead(LD),v);
}

// use a finer ordering

proc SannfsLOT(poly F, list #)
"USAGE:  SannfsLOT(f [,eng]);  f a poly, eng an optional int
RETURN:  ring
PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the
@* Levandovskyy's modification of the algorithm by Oaku and Takayama in D[s]
NOTE:    activate the output ring with the @code{setring} command.
@*    In the ring D[s], the ideal LD (which is NOT a Groebner basis) is
@*    the needed D-module structure.
@*    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 SannfsLOT; shows examples
"
{
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  // returns a list with a ring and an ideal LD in it
  int ppl = printlevel-voice+2;
  //  printf("plevel :%s, voice: %s",printlevel,voice);
  def save = basering;
  int N = nvars(basering);
  //  int Nnew = 2*(N+2);
  int Nnew = 2*(N+1)+1; //removed u,v; added s
  int i,j;
  string s;
  list RL = ringlist(basering);
  list L, Lord;
  list tmp;
  intvec iv;
  L[1] = RL[1]; // char
  L[4] = RL[4]; // char, minpoly
  // check whether vars have admissible names
  list Name  = RL[2];
  list RName;
//   RName[1] = "u";
//   RName[2] = "v";
  RName[1] = "t";
  RName[2] = "Dt";
  for(i=1;i<=N;i++)
  {
    for(j=1; j<=size(RName);j++)
    {
      if (Name[i] == RName[j])
      {
        ERROR("Variable names should not include t,Dt");
      }
    }
  }
  // now, create the names for new vars
//   tmp[1]     = "u";
//   tmp[2]     = "v";
//   list UName = tmp;
  list DName;
  for(i=1;i<=N;i++)
  {
    DName[i] = "D"+Name[i]; // concat
  }
  tmp    =  0;
  tmp[1] = "t";
  tmp[2] = "Dt";
  list SName ; SName[1] = "s";
  //  list NName = tmp + Name + DName + SName;
  list NName = tmp + SName + Name + DName;
  L[2]   = NName;
  tmp    = 0;
  // Name, Dname will be used further
  //  kill UName;
  kill NName;
  // block ord (a(1,1),dp);
  tmp[1]  = "a"; // string
  iv      = 1,1;
  tmp[2]  = iv; //intvec
  Lord[1] = tmp;
  // continue with a(0,0,1)
  tmp[1]  = "a"; // string
  iv      = 0,0,1;
  tmp[2]  = iv; //intvec
  Lord[2] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1]  = "dp"; // string
  s       = "iv=";
  for(i=1;i<=Nnew;i++)
  {
    s = s+"1,";
  }
  s[size(s)]= ";";
  execute(s);
  tmp[2]    = iv;
  Lord[3]   = tmp;
  tmp[1]    = "C";
  iv        = 0;
  tmp[2]    = iv;
  Lord[4]   = tmp;
  tmp       = 0;
  L[3]      = Lord;
  // we are done with the list
  def @R@ = ring(L);
  setring @R@;
  matrix @D[Nnew][Nnew];
  @D[1,2]=1;
  for(i=1; i<=N; i++)
  {
    @D[3+i,N+3+i]=1;
  }
  // ADD [s,t]=-t, [s,Dt]=Dt
  @D[1,3] = -var(1);
  @D[2,3] = var(2);
  //  @D[N+3,2*(N+2)]=1; old t,Dt stuff
  //  L[5] = matrix(UpOneMatrix(Nnew));
  //  L[6] = @D;
  def @R = nc_algebra(1,@D);
  setring @R;
  kill @R@;
  dbprint(ppl,"// -1-1- the ring @R(t,Dt,s,_x,_Dx) is ready");
  dbprint(ppl-1, @R);
  // create the ideal I
  poly  F = imap(save,F);
  //  ideal I = u*F-t,u*v-1;
  ideal I = F-t;
  poly p;
  for(i=1; i<=N; i++)
  {
    //    p = u*Dt; // u*Dt
    p = Dt;
    p = diff(F,var(3+i))*p;
    I = I, var(N+3+i) + p;
  }
  //  I = I, var(1)*var(2) + var(Nnew) +1; // reduce it with t-f!!!
  // t*Dt + s +1 reduced with t-f gives f*Dt + s
  I = I, F*var(2) + var(3);
  // -------- the ideal I is ready ----------
  dbprint(ppl,"// -1-2- starting the elimination of t,Dt in @R");
  dbprint(ppl-1, I);
  ideal J = engine(I,eng);
  ideal K = nselect(J,1..2);
  dbprint(ppl,"// -1-3- t,Dt are eliminated");
  dbprint(ppl-1, K);  // K is without t, Dt
  K = engine(K,eng);  // std does the job too
  // now, we must change the ordering
  // and create a ring without t, Dt
  setring save;
  // ----------- the ring @R3 ------------
  // _x, _Dx,s;  elim.ord for _x,_Dx.
  // keep: N, i,j,s, tmp, RL
  Nnew = 2*N+1;
  kill Lord, tmp, iv, RName;
  list Lord, tmp;
  intvec iv;
  L[1] = RL[1];
  L[4] = RL[4];  // char, minpoly
  // check whether vars hava admissible names -> done earlier
  // now, create the names for new var
  tmp[1] = "s";
  // DName is defined earlier
  list NName = Name + DName + tmp;
  L[2] = NName;
  tmp = 0;
  // block ord (dp(N),dp);
  // string s is already defined
  s = "iv=";
  for (i=1; i<=Nnew-1; i++)
  {
    s = s+"1,";
  }
  s[size(s)]=";";
  execute(s);
  tmp[1] = "dp";  // string
  tmp[2] = iv;   // intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1] = "dp";  // string
  s[size(s)] = ",";
  s = s+"1;";
  execute(s);
  kill s;
  kill NName;
  tmp[2]      = iv;
  Lord[2]     = tmp;
  tmp[1]      = "C";  iv  = 0;  tmp[2]=iv;
  Lord[3]     = tmp;  tmp = 0;
  L[3]        = Lord;
  // we are done with the list. Now add a Plural part
  def @R2@ = ring(L);
  setring @R2@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R2 = nc_algebra(1,@D);
  setring @R2;
  kill @R2@;
  dbprint(ppl,"//  -2-1- the ring @R2(_x,_Dx,s) is ready");
  dbprint(ppl-1, @R2);
  ideal MM = maxideal(1);
  //  MM = 0,s,MM;
  MM = 0,0,s,MM[1..size(MM)-1];
  map R01 = @R, MM;
  ideal K = R01(K);
  // total cleanup
  ideal LD = K;
  // make leadcoeffs positive
  for (i=1; i<= ncols(LD); i++)
  {
    if (leadcoef(LD[i]) <0 )
    {
      LD[i] = -LD[i];
    }
  }
  export LD;
  kill @R;
  return(@R2);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),Dp;
  poly F = x^3+y^3+z^3;
  printlevel = 0;
  def A  = SannfsLOT(F);
  setring A;
  LD;
}

/*
proc SannfsLOTold(poly F, list #)
"USAGE:  SannfsLOT(f [,eng]);  f a poly, eng an optional int
RETURN:  ring
PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the Levandovskyy's modification of the algorithm by Oaku and Takayama in the ring D[s], where D is the Weyl algebra
NOTE:    activate the output ring with the @code{setring} command.
@*       In the ring D[s], the ideal LD (which is NOT a Groebner basis) is the needed D-module structure.
@*       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 SannfsLOT; shows examples
"
{
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  // returns a list with a ring and an ideal LD in it
  int ppl = printlevel-voice+2;
  //  printf("plevel :%s, voice: %s",printlevel,voice);
  def save = basering;
  int N = nvars(basering);
  //  int Nnew = 2*(N+2);
  int Nnew = 2*(N+1)+1; //removed u,v; added s
  int i,j;
  string s;
  list RL = ringlist(basering);
  list L, Lord;
  list tmp;
  intvec iv;
  L[1] = RL[1]; // char
  L[4] = RL[4]; // char, minpoly
  // check whether vars have admissible names
  list Name  = RL[2];
  list RName;
//   RName[1] = "u";
//   RName[2] = "v";
  RName[1] = "t";
  RName[2] = "Dt";
  for(i=1;i<=N;i++)
  {
    for(j=1; j<=size(RName);j++)
    {
      if (Name[i] == RName[j])
      {
        ERROR("Variable names should not include t,Dt");
      }
    }
  }
  // now, create the names for new vars
//   tmp[1]     = "u";
//   tmp[2]     = "v";
//   list UName = tmp;
  list DName;
  for(i=1;i<=N;i++)
  {
    DName[i] = "D"+Name[i]; // concat
  }
  tmp    =  0;
  tmp[1] = "t";
  tmp[2] = "Dt";
  list SName ; SName[1] = "s";
  //  list NName = UName +  tmp + Name + DName;
  list NName = tmp + Name + DName + SName;
  L[2]   = NName;
  tmp    = 0;
  // Name, Dname will be used further
  //  kill UName;
  kill NName;
  // block ord (a(1,1),dp);
  tmp[1]  = "a"; // string
  iv      = 1,1;
  tmp[2]  = iv; //intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1]  = "dp"; // string
  s       = "iv=";
  for(i=1;i<=Nnew;i++)
  {
    s = s+"1,";
  }
  s[size(s)]= ";";
  execute(s);
  tmp[2]    = iv;
  Lord[2]   = tmp;
  tmp[1]    = "C";
  iv        = 0;
  tmp[2]    = iv;
  Lord[3]   = tmp;
  tmp       = 0;
  L[3]      = Lord;
  // we are done with the list
  def @R@ = ring(L);
  setring @R@;
  matrix @D[Nnew][Nnew];
  @D[1,2]=1;
  for(i=1; i<=N; i++)
  {
    @D[2+i,N+2+i]=1;
  }
  // ADD [s,t]=-t, [s,Dt]=Dt
  @D[1,Nnew] = -var(1);
  @D[2,Nnew] = var(2);
  //  @D[N+3,2*(N+2)]=1; old t,Dt stuff
  //  L[5] = matrix(UpOneMatrix(Nnew));
  //  L[6] = @D;
  def @R = nc_algebra(1,@D);
  setring @R;
  kill @R@;
  dbprint(ppl,"// -1-1- the ring @R(t,Dt,_x,_Dx,s) is ready");
  dbprint(ppl-1, @R);
  // create the ideal I
  poly  F = imap(save,F);
  //  ideal I = u*F-t,u*v-1;
  ideal I = F-t;
  poly p;
  for(i=1; i<=N; i++)
  {
    //    p = u*Dt; // u*Dt
    p = Dt;
    p = diff(F,var(2+i))*p;
    I = I, var(N+2+i) + p;
  }
  //  I = I, var(1)*var(2) + var(Nnew) +1; // reduce it with t-f!!!
  // t*Dt + s +1 reduced with t-f gives f*Dt + s
  I = I, F*var(2) + var(Nnew);
  // -------- the ideal I is ready ----------
  dbprint(ppl,"// -1-2- starting the elimination of t,Dt in @R");
  dbprint(ppl-1, I);
  ideal J = engine(I,eng);
  ideal K = nselect(J,1..2);
  dbprint(ppl,"// -1-3- t,Dt are eliminated");
  dbprint(ppl-1, K);  // K is without t, Dt
  K = engine(K,eng);  // std does the job too
  // now, we must change the ordering
  // and create a ring without t, Dt
  setring save;
  // ----------- the ring @R3 ------------
  // _x, _Dx,s;  elim.ord for _x,_Dx.
  // keep: N, i,j,s, tmp, RL
  Nnew = 2*N+1;
  kill Lord, tmp, iv, RName;
  list Lord, tmp;
  intvec iv;
  L[1] = RL[1];
  L[4] = RL[4];  // char, minpoly
  // check whether vars hava admissible names -> done earlier
  // now, create the names for new var
  tmp[1] = "s";
  // DName is defined earlier
  list NName = Name + DName + tmp;
  L[2] = NName;
  tmp = 0;
  // block ord (dp(N),dp);
  // string s is already defined
  s = "iv=";
  for (i=1; i<=Nnew-1; i++)
  {
    s = s+"1,";
  }
  s[size(s)]=";";
  execute(s);
  tmp[1] = "dp";  // string
  tmp[2] = iv;   // intvec
  Lord[1] = tmp;
  // continue with dp 1,1,1,1...
  tmp[1] = "dp";  // string
  s[size(s)] = ",";
  s = s+"1;";
  execute(s);
  kill s;
  kill NName;
  tmp[2]      = iv;
  Lord[2]     = tmp;
  tmp[1]      = "C";  iv  = 0;  tmp[2]=iv;
  Lord[3]     = tmp;  tmp = 0;
  L[3]        = Lord;
  // we are done with the list. Now add a Plural part
  def @R2@ = ring(L);
  setring @R2@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R2 = nc_algebra(1,@D);
  setring @R2;
  kill @R2@;
  dbprint(ppl,"//  -2-1- the ring @R2(_x,_Dx,s) is ready");
  dbprint(ppl-1, @R2);
  ideal MM = maxideal(1);
  MM = 0,s,MM;
  map R01 = @R, MM;
  ideal K = R01(K);
  // total cleanup
  ideal LD = K;
  // make leadcoeffs positive
  for (i=1; i<= ncols(LD); i++)
  {
    if (leadcoef(LD[i]) <0 )
    {
      LD[i] = -LD[i];
    }
  }
  export LD;
  kill @R;
  return(@R2);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),Dp;
  poly F = x^3+y^3+z^3;
  printlevel = 0;
  def A  = SannfsLOTold(F);
  setring A;
  LD;
}

*/

proc annfsLOT(poly F, list #)
"USAGE:  annfsLOT(F [,eng]);  F a poly, eng an optional int
RETURN:  ring
PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to
@* the Levandovskyy's modification of the algorithm by Oaku and Takayama
NOTE:    activate the output ring with the @code{setring} command. In this ring,
@*  - the ideal LD (which is a Groebner basis) is the needed D-module structure,
@*    which is obtained by substituting the minimal integer root of a Bernstein
@*    polynomial into the s-parametric ideal;
@*  - the list BS contains the roots with multiplicities of BS polynomial of f.
@*    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 annfsLOT; shows examples
"
{
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  printlevel=printlevel+1;
  def save = basering;
  def @A = SannfsLOT(F,eng);
  setring @A;
  poly F = imap(save,F);
  def B  = annfs0(LD,F,eng);
  return(B);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),Dp;
  poly F = z*x^2+y^3;
  printlevel = 0;
  def A  = annfsLOT(F);
  setring A;
  LD;
  BS;
}

proc annfs0(ideal I, poly F, list #)
"USAGE:  annfs0(I, F [,eng]);  I an ideal, F a poly, eng an optional int
RETURN:  ring
PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, based
@* on the output of Sannfs-like procedure
NOTE:    activate the output ring with the @code{setring} command. In this ring,
@* - the ideal LD (which is a Groebner basis) is the annihilator of f^s,
@* - the list BS contains the roots with multiplicities of BS polynomial of f.
@*  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 annfs0; shows examples
"
{
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  def @R2 = basering;
  // we're in D_n[s], where the elim ord for s is set
  ideal J = NF(I,std(F));
  // make leadcoeffs positive
  int i;
  for (i=1; i<= ncols(J); i++)
  {
    if (leadcoef(J[i]) <0 )
    {
      J[i] = -J[i];
    }
  }
  J = J,F;
  ideal M = engine(J,eng);
  int Nnew = nvars(@R2);
  ideal K2 = nselect(M,1..Nnew-1);
  int ppl = printlevel-voice+2;
  dbprint(ppl,"// -1-1- _x,_Dx are eliminated in basering");
  dbprint(ppl-1, K2);
  // the ring @R3 and the search for minimal negative int s
  ring @R3 = 0,s,dp;
  dbprint(ppl,"// -2-1- the ring @R3 i.e. K[s] is ready");
  ideal K3 = imap(@R2,K2);
  poly p = K3[1];
  dbprint(ppl,"// -2-2- factorization");
  //  ideal P = factorize(p,1);  //without constants and multiplicities
  //  "--------- b-function factorizes into ---------";   P;
  // convert factors to the list of their roots with mults
  // assume all factors are linear
  //  ideal BS = normalize(P);
  //  BS = subst(BS,s,0);
  //  BS = -BS;
  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];
  }
  int sP = minIntRoot(bs,1);
  bs =  normalize(bs);
  bs = -subst(bs,s,0);
  dbprint(ppl,"// -2-3- minimal integer root found");
  dbprint(ppl-1, sP);
 //TODO: sort BS!
  // --------- substitute s found in the ideal ---------
  // --------- going back to @R and substitute ---------
  setring @R2;
  K2 = subst(I,s,sP);
  // create the ordinary Weyl algebra and put the result into it,
  // thus creating the ring @R5
  // keep: N, i,j,s, tmp, RL
  Nnew = Nnew - 1; // former 2*N;
  // list RL = ringlist(save);  // is defined earlier
  //  kill Lord, tmp, iv;
  list L = 0;
  list Lord, tmp;
  intvec iv;
  list RL = ringlist(basering);
  L[1] = RL[1];
  L[4] = RL[4];  //char, minpoly
  // check whether vars have admissible names -> done earlier
  // list Name = RL[2]M
  // DName is defined earlier
  list NName; // = RL[2]; // skip the last var 's'
  for (i=1; i<=Nnew; i++)
  {
    NName[i] =  RL[2][i];
  }
  L[2] = NName;
  // dp ordering;
  string s = "iv=";
  for (i=1; i<=Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)] = ";";
  execute(s);
  tmp     = 0;
  tmp[1]  = "dp";  // string
  tmp[2]  = iv;  // intvec
  Lord[1] = tmp;
  kill s;
  tmp[1]  = "C";
  iv = 0;
  tmp[2]  = iv;
  Lord[2] = tmp;
  tmp     = 0;
  L[3]    = Lord;
  // we are done with the list
  // Add: Plural part
  def @R4@ = ring(L);
  setring @R4@;
  int N = Nnew div 2;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R4 = nc_algebra(1,@D);
  setring @R4;
  kill @R4@;
  dbprint(ppl,"// -3-1- the ring @R4 is ready");
  dbprint(ppl-1, @R4);
  ideal K4 = imap(@R2,K2);
  intvec saveopt=option(get);
  option(redSB);
  dbprint(ppl,"// -3-2- the final cosmetic std");
  K4 = engine(K4,eng);  // std does the job too
  // total cleanup
  ideal bs = imap(@R3,bs);
  kill @R3;
  list BS = bs,m;
  export BS;
  ideal LD = K4;
  export LD;
  option(set,saveopt);
  return(@R4);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),Dp;
  poly F = x^3+y^3+z^3;
  printlevel = 0;
  def A = SannfsBM(F);   setring A;
  // alternatively, one can use SannfsOT or SannfsLOT
  LD;
  poly F = imap(r,F);
  def B  = annfs0(LD,F);  setring B;
  LD;
  BS;
}

// proc annfsgms(poly F, list #)
// "USAGE:  annfsgms(f [,eng]);  f a poly, eng an optional int
// ASSUME:  f has an isolated critical point at 0
// RETURN:  ring
// PURPOSE: compute the D-module structure of basering[1/f]*f^s
// NOTE:    activate the output ring with the @code{setring} command. In this ring,
// @*       - the ideal LD is the needed D-mod structure,
// @*       - the ideal BS is the list of roots of a Bernstein polynomial of f.
// @*       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 annfsgms; shows examples
// "
// {
//   LIB "gmssing.lib";
//   int eng = 0;
//   if ( size(#)>0 )
//   {
//     if ( typeof(#[1]) == "int" )
//     {
//       eng = int(#[1]);
//     }
//   }
//   int ppl = printlevel-voice+2;
//   // returns a ring with the ideal LD in it
//   def save = basering;
//   // compute the Bernstein polynomial from gmssing.lib
//   list RL = ringlist(basering);
//   // in the descr. of the ordering, replace "p" by "s"
//   list NL = convloc(RL);
//   // create a ring with the ordering, converted to local
//   def @LR = ring(NL);
//   setring @LR;
//   poly F  = imap(save, F);
//   ideal B = bernstein(F)[1];
//   // since B may not contain (s+1) [following gmssing.lib]
//   // add it!
//   B = B,-1;
//   B = simplify(B,2+4); // erase zero and repeated entries
//   // find the minimal integer value
//   int   S = minIntRoot(B,0);
//   dbprint(ppl,"// -0- minimal integer root found");
//   dbprint(ppl-1,S);
//   setring save;
//   int N = nvars(basering);
//   int Nnew = 2*(N+2);
//   int i,j;
//   string s;
//   //  list RL = ringlist(basering);
//   list L, Lord;
//   list tmp;
//   intvec iv;
//   L[1] = RL[1]; // char
//   L[4] = RL[4]; // char, minpoly
//   // check whether vars have admissible names
//   list Name  = RL[2];
//   list RName;
//   RName[1] = "u";
//   RName[2] = "v";
//   RName[3] = "t";
//   RName[4] = "Dt";
//   for(i=1;i<=N;i++)
//   {
//     for(j=1; j<=size(RName);j++)
//     {
//       if (Name[i] == RName[j])
//       {
//         ERROR("Variable names should not include u,v,t,Dt");
//       }
//     }
//   }
//   // now, create the names for new vars
//   //  tmp[1]     = "u";  tmp[2]     = "v";  tmp[3]     = "t";  tmp[4]     = "Dt";
//   list UName = RName;
//   list DName;
//   for(i=1;i<=N;i++)
//   {
//     DName[i] = "D"+Name[i]; // concat
//   }
//   list NName = UName + Name + DName;
//   L[2]       = NName;
//   tmp        = 0;
//   // Name, Dname will be used further
//   kill UName;
//   kill NName;
//   // block ord (a(1,1),dp);
//   tmp[1]  = "a"; // string
//   iv      = 1,1;
//   tmp[2]  = iv; //intvec
//   Lord[1] = tmp;
//   // continue with dp 1,1,1,1...
//   tmp[1]  = "dp"; // string
//   s       = "iv=";
//   for(i=1; i<=Nnew; i++) // need really all vars!
//   {
//     s = s+"1,";
//   }
//   s[size(s)]= ";";
//   execute(s);
//   tmp[2]    = iv;
//   Lord[2]   = tmp;
//   tmp[1]    = "C";
//   iv        = 0;
//   tmp[2]    = iv;
//   Lord[3]   = tmp;
//   tmp       = 0;
//   L[3]      = Lord;
//   // we are done with the list
//   def @R = ring(L);
//   setring @R;
//   matrix @D[Nnew][Nnew];
//   @D[3,4] = 1; // t,Dt
//   for(i=1; i<=N; i++)
//   {
//     @D[4+i,4+N+i]=1;
//   }
//   //  L[5] = matrix(UpOneMatrix(Nnew));
//   //  L[6] = @D;
//   nc_algebra(1,@D);
//   dbprint(ppl,"// -1-1- the ring @R is ready");
//   dbprint(ppl-1,@R);
//   // create the ideal
//   poly F  = imap(save,F);
//   ideal I = u*F-t,u*v-1;
//   poly p;
//   for(i=1; i<=N; i++)
//   {
//     p = u*Dt; // u*Dt
//     p = diff(F,var(4+i))*p;
//     I = I, var(N+4+i) + p; // Dx, Dy
//   }
//   // add the relations between t,Dt and s
//   //  I = I, t*Dt+1+S;
//   // -------- the ideal I is ready ----------
//   dbprint(ppl,"// -1-2- starting the elimination of u,v in @R");
//   ideal J = engine(I,eng);
//   ideal K = nselect(J,1..2);
//   dbprint(ppl,"// -1-3- u,v are eliminated in @R");
//   dbprint(ppl-1,K); // without u,v: not yet our answer
//   //----- create a ring with elim.ord for t,Dt -------
//   setring save;
//   // ------------ new ring @R2 ------------------
//   // without u,v and with the elim.ord for t,Dt
//   // keep: N, i,j,s, tmp, RL
//   Nnew = 2*N+2;
//   //  list RL = ringlist(save); // is defined earlier
//   kill Lord,tmp,iv, RName;
//   L = 0;
//   list Lord, tmp;
//   intvec iv;
//   L[1] = RL[1]; // char
//   L[4] = RL[4]; // char, minpoly
//   // check whether vars have admissible names -> done earlier
//   //  list Name  = RL[2];
//   list RName;
//   RName[1] = "t";
//   RName[2] = "Dt";
//   // DName is defined earlier
//   list NName = RName + Name + DName;
//   L[2]   = NName;
//   tmp    = 0;
//   // block ord (a(1,1),dp);
//   tmp[1]  = "a"; // string
//   iv      = 1,1;
//   tmp[2]  = iv; //intvec
//   Lord[1] = tmp;
//   // continue with dp 1,1,1,1...
//   tmp[1]  = "dp"; // string
//   s       = "iv=";
//   for(i=1;i<=Nnew;i++)
//   {
//     s = s+"1,";
//   }
//   s[size(s)]= ";";
//   execute(s);
//   kill s;
//   kill NName;
//   tmp[2]    = iv;
//   Lord[2]   = tmp;
//   tmp[1]    = "C";
//   iv        = 0;
//   tmp[2]    = iv;
//   Lord[3]   = tmp;
//   tmp       = 0;
//   L[3]      = Lord;
//   // we are done with the list
//   // Add: Plural part
//   def @R2 = ring(L);
//   setring @R2;
//   matrix @D[Nnew][Nnew];
//   @D[1,2]=1;
//   for(i=1; i<=N; i++)
//   {
//     @D[2+i,2+N+i]=1;
//   }
//   nc_algebra(1,@D);
//   dbprint(ppl,"// -2-1- the ring @R2 is ready");
//   dbprint(ppl-1,@R2);
//   ideal MM = maxideal(1);
//   MM = 0,0,MM;
//   map R01 = @R, MM;
//   ideal K2 = R01(K);
//   // add the relations between t,Dt and s
//   //  K2       = K2, t*Dt+1+S;
//   poly G = t*Dt+S+1;
//   K2 = NF(K2,std(G)),G;
//   dbprint(ppl,"// -2-2- starting elimination for t,Dt in @R2");
//   ideal J  = engine(K2,eng);
//   ideal K  = nselect(J,1..2);
//   dbprint(ppl,"// -2-3- t,Dt are eliminated");
//   dbprint(ppl-1,K);
//   //  "------- produce a final result --------";
//   // ----- create the ordinary Weyl algebra and put
//   // ----- the result into it
//   // ------ create the ring @R5
//   // keep: N, i,j,s, tmp, RL
//   setring save;
//   Nnew = 2*N;
//   //  list RL = ringlist(save); // is defined earlier
//   kill Lord, tmp, iv;
//   //  kill L;
//   L = 0;
//   list Lord, tmp;
//   intvec iv;
//   L[1] = RL[1]; // char
//   L[4] = RL[4]; // char, minpoly
//   // check whether vars have admissible names -> done earlier
//   //  list Name  = RL[2];
//   // DName is defined earlier
//   list NName = Name + DName;
//   L[2]   = NName;
//   // dp ordering;
//   string   s       = "iv=";
//   for(i=1;i<=2*N;i++)
//   {
//     s = s+"1,";
//   }
//   s[size(s)]= ";";
//   execute(s);
//   tmp     = 0;
//   tmp[1]  = "dp"; // string
//   tmp[2]  = iv; //intvec
//   Lord[1] = tmp;
//   kill s;
//   tmp[1]    = "C";
//   iv        = 0;
//   tmp[2]    = iv;
//   Lord[2]   = tmp;
//   tmp       = 0;
//   L[3]      = Lord;
//   // we are done with the list
//   // Add: Plural part
//   def @R5 = ring(L);
//   setring @R5;
//   matrix @D[Nnew][Nnew];
//   for(i=1; i<=N; i++)
//   {
//     @D[i,N+i]=1;
//   }
//   nc_algebra(1,@D);
//   dbprint(ppl,"// -3-1- the ring @R5 is ready");
//   dbprint(ppl-1,@R5);
//   ideal K5 = imap(@R2,K);
//   option(redSB);
//   dbprint(ppl,"// -3-2- the final cosmetic std");
//   K5 = engine(K5,eng); // std does the job too
//   // total cleanup
//   kill @R;
//   kill @R2;
//   ideal LD = K5;
//   ideal BS = imap(@LR,B);
//   kill @LR;
//   export BS;
//   export LD;
//   return(@R5);
// }
// example
// {
//   "EXAMPLE:"; echo = 2;
//   ring r = 0,(x,y,z),Dp;
//   poly F = x^2+y^3+z^5;
//   def A = annfsgms(F);
//   setring A;
//   LD;
//   print(matrix(BS));
// }


proc convloc(list @NL)
"USAGE:  convloc(L); L a list
RETURN:  list
PURPOSE: convert a ringlist L into another ringlist,
@* where all the 'p' orderings are replaced with the 's' orderings, e.g. @code{dp} by @code{ds}.
ASSUME: L is a result of a ringlist command
EXAMPLE: example convloc; shows examples
"
{
  list NL = @NL;
  // given a ringlist, returns a new ringlist,
  // where all the p-orderings are replaced with s-ord's
  int i,j,k,found;
  int nblocks = size(NL[3]);
  for(i=1; i<=nblocks; i++)
  {
    for(j=1; j<=size(NL[3][i]); j++)
    {
      if (typeof(NL[3][i][j]) == "string" )
      {
        found = 0;
        for (k=1; k<=size(NL[3][i][j]); k++)
        {
          if (NL[3][i][j][k]=="p")
          {
            NL[3][i][j][k]="s";
            found = 1;
            //    printf("replaced at %s,%s,%s",i,j,k);
          }
        }
      }
    }
  }
  return(NL);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),(Dp(2),dp(1));
  list L = ringlist(r);
  list N = convloc(L);
  def rs = ring(N);
  setring rs;
  rs;
}

proc annfspecial(ideal I, poly F, int mir, number n)
"USAGE:  annfspecial(I,F,mir,n);  I an ideal, F a poly, int mir, number n
RETURN:  ideal
PURPOSE: compute the annihilator ideal of F^n in the Weyl Algebra
@*           for the given rational number n
ASSUME:  The basering is D[s] and contains 's' explicitly as a variable,
@*   the ideal I is the Ann F^s in D[s] (obtained with e.g. SannfsBM),
@*   the integer 'mir' is the minimal integer root of the BS polynomial of F,
@*   and the number n is rational.
NOTE: We compute the real annihilator for any rational value of n (both
@*       generic and exceptional). The implementation goes along the lines of
@*       the Algorithm 5.3.15 from Saito-Sturmfels-Takayama.
DISPLAY: If printlevel=1, progress debug messages will be printed,
@*       if printlevel>=2, all the debug messages will be printed.
EXAMPLE: example annfspecial; shows examples
"
{

  if (!isRational(n))
  {
    "ERROR: n must be a rational number!";
  }
  int ppl = printlevel-voice+2;
  //  int mir =  minIntRoot(L[1],0);
  int ns   = varNum("s");
  if (!ns)
  {
    ERROR("Variable s expected in the ideal AnnFs");
  }
  int d;
  ideal P = subst(I,var(ns),n);
  if (denominator(n) == 1)
  {
    // n is fraction-free
    d = int(numerator(n));
    if ( (!d) && (n!=0))
    {
      ERROR("no parametric special values are allowed");
    }
    d = d - mir;
    if (d>0)
    {
      dbprint(ppl,"// -1-1- starting syzygy computations");
      matrix J[1][1] = F^d;
      dbprint(ppl-1,"// -1-1-1- of the polynomial ideal");
      dbprint(ppl-1,J);
      matrix K[1][size(I)] = subst(I,var(ns),mir);
      dbprint(ppl-1,"// -1-1-2- modulo ideal:");
      dbprint(ppl-1, K);
      module M = modulo(J,K);
      dbprint(ppl-1,"// -1-1-3- getting the result:");
      dbprint(ppl-1, M);
      P  = P, ideal(M);
      dbprint(ppl,"// -1-2- finished syzygy computations");
    }
    else
    {
      dbprint(ppl,"// -1-1- d<=0, no syzygy computations needed");
      dbprint(ppl-1,"// -1-2- for d =");
      dbprint(ppl-1, d);
    }
  }
  // also the else case: d<=0 or n is not an integer
  dbprint(ppl,"// -2-1- starting final Groebner basis");
  P = groebner(P);
  dbprint(ppl,"// -2-2- finished final Groebner basis");
  return(P);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  poly F = x3-y2;
  def  B = annfs(F);  setring B;
  minIntRoot(BS[1],0);
  // So, the minimal integer root is -1
  setring r;
  def  A = SannfsBM(F);
  setring A;
  poly F = x3-y2;
  annfspecial(LD,F,-1,3/4); // generic root
  annfspecial(LD,F,-1,-2);  // integer but still generic root
  annfspecial(LD,F,-1,1);   // exceptional root
}

/*
//static proc new_ex_annfspecial()
{
  //another example for annfspecial: x3+y3+z3
  ring r = 0,(x,y,z),dp;
  poly F =  x3+y3+z3;
  def  B = annfs(F);  setring B;
  minIntRoot(BS[1],0);
  // So, the minimal integer root is -1
  setring r;
  def  A = SannfsBM(F);
  setring A;
  poly F =  x3+y3+z3;
  annfspecial(LD,F,-1,3/4); // generic root
  annfspecial(LD,F,-1,-2);  // integer but still generic root
  annfspecial(LD,F,-1,1);   // exceptional root
}
*/

proc minIntRoot(ideal P, int fact)
"USAGE:  minIntRoot(P, fact); P an ideal, fact an int
RETURN:  int
PURPOSE: minimal integer root of a maximal ideal P
NOTE:    if fact==1, P is the result of some 'factorize' call,
@*       else P is treated as the result of bernstein::gmssing.lib
@*       in both cases without constants and multiplicities
EXAMPLE: example minIntRoot; shows examples
"
{
  //    ideal P = factorize(p,1);
  // or ideal P = bernstein(F)[1];
  intvec vP;
  number nP;
  int i;
  if ( fact )
  {
    // the result comes from "factorize"
    P = normalize(P); // now leadcoef = 1
    // TODO: hunt for units and kill then !!!
    P = matrix(lead(P))-P;
    // nP = leadcoef(P[i]-lead(P[i])); // for 1 var only, extract the coeff
  }
  else
  {
    // bernstein deletes -1 usually
    P = P, -1;
  }
  // for both situations:
  // now we have an ideal of fractions of type "number"
  int sP = size(P);
  for(i=1; i<=sP; i++)
  {
    nP = leadcoef(P[i]);
    if ( (nP - int(nP)) == 0 )
    {
      vP = vP,int(nP);
    }
  }
  if ( size(vP)>=2 )
  {
    vP = vP[2..size(vP)];
  }
  sP = -Max(-vP);
  if (sP == 0)
  {
    "Warning: zero root present!";
  }
  return(sP);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r   = 0,(x,y),ds;
  poly f1  = x*y*(x+y);
  ideal I1 = bernstein(f1)[1]; // a local Bernstein poly
  I1;
  minIntRoot(I1,0);
  poly  f2  = x2-y3;
  ideal I2  = bernstein(f2)[1];
  I2;
  minIntRoot(I2,0);
  // now we illustrate the behaviour of factorize
  // together with a global ordering
  ring r2  = 0,x,dp;
  poly f3   = 9*(x+2/3)*(x+1)*(x+4/3); //global b-polynomial of f1=x*y*(x+y)
  ideal I3 = factorize(f3,1);
  I3;
  minIntRoot(I3,1);
  // and a more interesting situation
  ring  s  = 0,(x,y,z),ds;
  poly  f  = x3 + y3 + z3;
  ideal I  = bernstein(f)[1];
  I;
  minIntRoot(I,0);
}

proc isHolonomic(def M)
"USAGE:  isHolonomic(M); M an ideal/module/matrix
RETURN:  int, 1 if M is holonomic over the base ring, and 0 otherwise
ASSUME: basering is a Weyl algebra in characteristic 0
PURPOSE: check whether M is holonomic over the base ring
NOTE:    M is holonomic if 2*dim(M) = dim(R), where R is the
base ring; dim stands for Gelfand-Kirillov dimension
EXAMPLE: example isHolonomic; shows examples
"
{
  // char>0 or qring
  if (  (size(ideal(basering)) >0) || (char(basering) >0) )
  {
    ERROR("Basering is inappropriate: characteristic>0 or qring present");
  }
  if (!isWeyl(basering))
  {
    ERROR("Basering is not a Weyl algebra");
  }

  if ( (typeof(M) != "ideal") && (typeof(M) != "module") && (typeof(M) != "matrix") )
  {
  //  print(typeof(M));
    ERROR("an argument of type ideal/module/matrix expected");
  }
  if (attrib(M,"isSB")!=1)
  {
    M = engine(M,0);
  }
  int dimR = gkdim(std(0));
  int dimM = gkdim(M);
  return( (dimR==2*dimM) );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring R = 0,(x,y),dp;
  poly F = x*y*(x+y);
  def A = annfsBM(F,0);
  setring A;
  LD;
  isHolonomic(LD);
  ideal I = std(LD[1]);
  I;
  isHolonomic(I);
}

proc reiffen(int p, int q)
"USAGE:  reiffen(p, q);  int p, int q
RETURN:  ring
PURPOSE: set up the polynomial, describing a Reiffen curve
NOTE:  activate the output ring with the @code{setring} command and
@*       find the curve as a polynomial RC.
@*  A Reiffen curve is defined as RC = x^p + y^q + xy^{q-1}, q >= p+1 >= 5

EXAMPLE: example reiffen; shows examples
"
{
  // we allow also other numbers, p \geq 1, q\geq 1
// a Reiffen curve is defined as
// F = x^p + y^q +x*y^{q-1}, q \geq p+1 \geq 5

// ASSUME: q >= p+1 >= 5. Otherwise an error message is returned

//   if ( (p<4) || (q<5) || (q-p<1) )
//   {
//     ERROR("Some of conditions p>=4, q>=5 or q>=p+1 is not satisfied!");
//   }
  if ( (p<1) || (q<1) )
  {
    ERROR("Some of conditions p>=1, q>=1 is not satisfied!");
  }
  ring @r = 0,(x,y),dp;
  poly RC = y^q +x^p + x*y^(q-1);
  export RC;
  return(@r);
}
example
{
  "EXAMPLE:"; echo = 2;
  def r = reiffen(4,5);
  setring r;
  RC;
}

proc arrange(int p)
"USAGE:  arrange(p);  int p
RETURN:  ring
PURPOSE: set up the polynomial, describing a hyperplane arrangement
NOTE:  must be executed in a commutative ring
ASSUME: basering is present and it is commutative
EXAMPLE: example arrange; shows examples
"
{
  int UseBasering = 0 ;
  if (p<2)
  {
    ERROR("p>=2 is required!");
  }
  if ( nameof(basering)!="basering" )
  {
    if (p > nvars(basering))
    {
      ERROR("too big p");
    }
    else
    {
      def @r = basering;
      UseBasering = 1;
    }
  }
  else
  {
    // no basering found
    ERROR("no basering found!");
    //    ring @r = 0,(x(1..p)),dp;
  }
  int i,j,sI;
  poly  q = 1;
  list ar;
  ideal tmp;
  tmp = ideal(var(1));
  ar[1] = tmp;
  for (i = 2; i<=p; i++)
  {
    // add i-nary sums to the product
    ar = indAR(ar,i);
  }
  for (i = 1; i<=p; i++)
  {
    tmp = ar[i]; sI = size(tmp);
    for (j = 1; j<=sI; j++)
    {
      q = q*tmp[j];
    }
  }
  if (UseBasering)
  {
    return(q);
  }
  //  poly AR = q; export AR;
  //  return(@r);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring X = 0,(x,y,z,t),dp;
  poly q = arrange(3);
  factorize(q,1);
}

proc checkRoot(poly F, number a, list #)
"USAGE:  checkRoot(f,alpha [,S,eng]);  poly f, number alpha, string S, int eng
RETURN:  int
ASSUME: Basering is a commutative ring, alpha is a positive rational number.
PURPOSE: check whether a negative rational number -alpha is a root of the global
@* Bernstein-Sato polynomial of f and compute its multiplicity,
@* with the algorithm given by S and with the Groebner basis engine given by eng.
NOTE:  The annihilator of f^s in D[s] is needed and hence it is computed with the
@*       algorithm by Briancon and Maisonobe. The value of a string S can be
@*       'alg1' (default) - for the algorithm 1 of [LM08]
@*       'alg2' - for the algorithm 2 of [LM08]
@*       Depending on the value of S, the output of type int is:
@*       'alg1': 1 only if -alpha is a root of the global Bernstein-Sato polynomial
@*       'alg2':  the multiplicity of -alpha as a root of the global Bernstein-Sato
@*                  polynomial of f. If -alpha is not a root, the output is 0.
@*       If eng <>0, @code{std} is used for Groebner basis computations,
@*       otherwise (and by default) @code{slimgb} is used.
DISPLAY: If printlevel=1, progress debug messages will be printed,
@*          if printlevel>=2, all the debug messages will be printed.
EXAMPLE: example checkRoot; shows examples
"
{
  int eng = 0;
  int chs = 0; // choice
  string algo = "alg1";
  string st;
  if ( size(#)>0 )
  {
   if ( typeof(#[1]) == "string" )
   {
     st = string(#[1]);
     if ( (st == "alg1") || (st == "ALG1") || (st == "Alg1") ||(st == "aLG1"))
     {
       algo = "alg1";
       chs  = 1;
     }
     if ( (st == "alg2") || (st == "ALG2") || (st == "Alg2") || (st == "aLG2"))
     {
       algo = "alg2";
       chs  = 1;
     }
     if (chs != 1)
     {
       // incorrect value of S
       print("Incorrect algorithm given, proceed with the default alg1");
       algo = "alg1";
     }
     // second arg
     if (size(#)>1)
     {
       // exists 2nd arg
       if ( typeof(#[2]) == "int" )
       {
         // the case: given alg, given engine
         eng = int(#[2]);
       }
       else
       {
         eng = 0;
       }
     }
     else
     {
       // no second arg
       eng = 0;
     }
   }
   else
   {
     if ( typeof(#[1]) == "int" )
     {
       // the case: default alg, engine
       eng = int(#[1]);
       // algo = "alg1";  //is already set
     }
     else
     {
       // incorr. 1st arg
       algo = "alg1";
     }
   }
  }
  // size(#)=0, i.e. there is no algorithm and no engine given
  //  eng = 0; algo = "alg1";  //are already set
  // int ppl = printlevel-voice+2;
// check assume: a is positive rational number
  if (!isRational(a))
  {
    ERROR("rational root expected for checking");
  }
  if (numerator(a) < 0 )
  {
    ERROR("expected positive -alpha");
    // the following is more user-friendly but less correct
    //    print("proceeding with the negated root");
    //    a = -a;
  }
  printlevel=printlevel+1;
  def save = basering;
  def @A = SannfsBM(F);
  setring @A;
  poly F = imap(save,F);
  number a = imap(save,a);
  if ( algo=="alg1")
  {
    int output = checkRoot1(LD,F,a,eng);
  }
  else
  {
    if ( algo=="alg2")
    {
      int output = checkRoot2(LD,F,a,eng);
    }
  }
  printlevel=printlevel-1;
  return(output);
}
example
{
  "EXAMPLE:"; echo = 2;
  printlevel=0;
  ring r = 0,(x,y),Dp;
  poly F = x^4+y^5+x*y^4;
  checkRoot(F,11/20);    // -11/20 is a root of bf
  poly G = x*y;
  checkRoot(G,1,"alg2"); // -1 is a root of bg with multiplicity 2
}

proc checkRoot1(ideal I, poly F, number a, list #)
"USAGE:  checkRoot1(I,f,alpha [,eng]);  ideal I, poly f, number alpha, int eng
ASSUME:  Basering is D[s], I is the annihilator of f^s in D[s],
@* that is basering and I are the output of Sannfs-like procedure,
@* f is a polynomial in K[x] and alpha is a rational number.
RETURN:  int, 1 if -alpha is a root of the Bernstein-Sato polynomial of f
PURPOSE: check, whether alpha is a root of the global Bernstein-Sato polynomial of f
NOTE:  If eng <>0, @code{std} is used for Groebner basis 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 checkRoot1; shows examples
"
{
  // to check: alpha is rational (has char 0 check inside)
  if (!isRational(a))
  {
    "ERROR: alpha must be a rational number!";
  }
  //  no qring
  if ( size(ideal(basering)) >0 )
  {
    "ERROR: no qring is allowed";
  }
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  int ppl = printlevel-voice+2;
  dbprint(ppl,"// -0-1- starting the procedure checkRoot1");
  def save = basering;
  int N = nvars(basering);
  int Nnew = N-1;
  int n = Nnew div 2;
  int i;
  string s;
  list RL = ringlist(basering);
  list L, Lord;
  list tmp;
  intvec iv;
  L[1] = RL[1];  // char
  L[4] = RL[4];  // char, minpoly
  // check whether basering is D[s]=K(_x,_Dx,s)
  list Name = RL[2];
//   for (i=1; i<=n; i++)
//   {
//     if ( bracket(var(i+n),var(i))!=1 )
//     {
//       ERROR("basering should be D[s]=K(_x,_Dx,s)");
//     }
//   }
  if ( Name[N]!="s" )
  {
    ERROR("the last variable of basering should be s");
  }
  // now, create the new vars
  list NName;
  for (i=1; i<=Nnew; i++)
  {
    NName[i] = Name[i];
  }
  L[2] = NName;
  kill Name,NName;
  // block ord (dp);
  tmp[1] = "dp"; // string
  s = "iv=";
  for (i=1; i<=Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)]=";";
  execute(s);
  kill s;
  tmp[2]  = iv;
  Lord[1] = tmp;
  tmp[1]  = "C";
  iv      = 0;
  tmp[2]  = iv;
  Lord[2] = tmp;
  tmp     = 0;
  L[3]    = Lord;
  // we are done with the list
  def @R@ = ring(L);
  setring @R@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=n; i++)
  {
    @D[i,i+n]=1;
  }
  def @R = nc_algebra(1,@D);
  setring @R;
  kill @R@;
  dbprint(ppl,"// -1-1- the ring @R(_x,_Dx) is ready");
  dbprint(ppl-1, S);
  // create the ideal K = ann_D[s](f^s)_{s=-alpha} + < f >
  setring save;
  ideal K = subst(I,s,-a);
  dbprint(ppl,"// -1-2- the variable s has been substituted by "+string(-a));
  dbprint(ppl-1, K);
  K = NF(K,std(F));
  // make leadcoeffs positive
  for (i=1; i<=ncols(K); i++)
  {
    if ( leadcoef(K[i])<0 )
    {
      K[i] = -K[i];
    }
  }
  K = K,F;
  // ------------ the ideal K is ready ------------
  setring @R;
  ideal K = imap(save,K);
  dbprint(ppl,"// -1-3- starting the computation of a Groebner basis of K in @R");
  dbprint(ppl-1, K);
  ideal G = engine(K,eng);
  dbprint(ppl,"// -1-4- the Groebner basis has been computed");
  dbprint(ppl-1, G);
  return(G[1]!=1);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),Dp;
  poly F = x^4+y^5+x*y^4;
  printlevel = 0;
  def A = Sannfs(F);
  setring A;
  poly F = imap(r,F);
  checkRoot1(LD,F,31/20);   // -31/20 is not a root of bs
  checkRoot1(LD,F,11/20);   // -11/20 is a root of bs
}

proc checkRoot2 (ideal I, poly F, number a, list #)
"USAGE:  checkRoot2(I,f,a [,eng]);  I an ideal, f a poly, alpha a number, eng an optional int
ASSUME:  I is the annihilator of f^s in D[s], basering is D[s],
@* that is basering and I are the output os Sannfs-like procedure,
@* f is a polynomial in K[_x] and alpha is a rational number.
RETURN:  int, the multiplicity of -alpha as a root of the BS polynomial of f.
PURPOSE: check whether a rational number alpha is a root of the global Bernstein-
@* Sato polynomial of f and compute its multiplicity from the known Ann F^s in D[s]
NOTE:   If -alpha is not a root, the output is 0.
@*       If eng <>0, @code{std} is used for Groebner basis computations,
@*       otherwise (and by default) @code{slimgb} is used.
DISPLAY: If printlevel=1, progress debug messages will be printed,
@*          if printlevel>=2, all the debug messages will be printed.
EXAMPLE: example checkRoot2; shows examples
"
{


  // to check: alpha is rational (has char 0 check inside)
  if (!isRational(a))
  {
    "ERROR: alpha must be a rational number!";
  }
  //  no qring
  if ( size(ideal(basering)) >0 )
  {
    "ERROR: no qring is allowed";
  }

  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  int ppl = printlevel-voice+2;
  dbprint(ppl,"// -0-1- starting the procedure checkRoot2");
  def save = basering;
  int N = nvars(basering);
  int n = (N-1) div 2;
  int i;
  string s;
  list RL = ringlist(basering);
  list L, Lord;
  list tmp;
  intvec iv;
  L[1] = RL[1];  // char
  L[4] = RL[4];  // char, minpoly
  // check whether basering is D[s]=K(_x,_Dx,s)
  list Name = RL[2];
  for (i=1; i<=n; i++)
  {
    if ( bracket(var(i+n),var(i))!=1 )
    {
      ERROR("basering should be D[s]=K(_x,_Dx,s)");
    }
  }
  if ( Name[N]!="s" )
  {
    ERROR("the last variable of basering should be s");
  }
  // now, create the new vars
  L[2] = Name;
  kill Name;
  // block ord (dp);
  tmp[1] = "dp"; // string
  s = "iv=";
  for (i=1; i<=N; i++)
  {
    s = s+"1,";
  }
  s[size(s)]=";";
  execute(s);
  kill s;
  tmp[2]  = iv;
  Lord[1] = tmp;
  tmp[1]  = "C";
  iv      = 0;
  tmp[2]  = iv;
  Lord[2] = tmp;
  tmp     = 0;
  L[3]    = Lord;
  // we are done with the list
  def @R@ = ring(L);
  setring @R@;
  matrix @D[N][N];
  for (i=1; i<=n; i++)
  {
    @D[i,i+n]=1;
  }
  def @R = nc_algebra(1,@D);
  setring @R;
  kill @R@;
  dbprint(ppl,"// -1-1- the ring @R(_x,_Dx,s) is ready");
  dbprint(ppl-1, @R);
  // now, continue with the algorithm
  ideal I = imap(save,I);
  poly F = imap(save,F);
  number a = imap(save,a);
  ideal II = NF(I,std(F));
  // make leadcoeffs positive
  for (i=1; i<=ncols(II); i++)
  {
    if ( leadcoef(II[i])<0 )
    {
      II[i] = -II[i];
    }
  }
  ideal J,G;
  int m;  // the output (multiplicity)
  dbprint(ppl,"// -2- starting the bucle");
  for (i=0; i<=n; i++)  // the multiplicity has to be <= n
  {
    // create the ideal Ji = ann_D[s](f^s) + < f, (s+alpha)^{i+1} >
    // (s+alpha)^i in Ji <==> -alpha is a root with multiplicity >= i
    J = II,F,(s+a)^(i+1);
    // ------------ the ideal Ji is ready -----------
    dbprint(ppl,"// -2-"+string(i+1)+"-1- starting the computation of a Groebner basis of J"+string(i)+" in @R");
    dbprint(ppl-1, J);
    G = engine(J,eng);
    dbprint(ppl,"// -2-"+string(i+1)+"-2- the Groebner basis has been computed");
    dbprint(ppl-1, G);
    if ( NF((s+a)^i,G)==0 )
    {
      dbprint(ppl,"// -2-"+string(i+1)+"-3- the number "+string(-a)+" has not multiplicity "+string(i+1));
      m = i;
      break;
    }
    dbprint(ppl,"// -2-"+string(i+1)+"-3- the number "+string(-a)+" has multiplicity at least "+string(i+1));
  }
  dbprint(ppl,"// -3- the bucle has finished");
  return(m);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),Dp;
  poly F = x*y*z;
  printlevel = 0;
  def A = Sannfs(F);
  setring A;
  poly F = imap(r,F);
  checkRoot2(LD,F,1);    // -1 is a root of bs with multiplicity 3
  checkRoot2(LD,F,1/3);  // -1/3 is not a root
}

proc checkFactor(ideal I, poly F, poly q, list #)
"USAGE:  checkFactor(I,f,qs [,eng]);  I an ideal, f a poly, qs a poly, eng an optional int
ASSUME:  checkFactor is called from the basering, created by Sannfs-like proc,
@* that is, from the Weyl algebra in x1,..,xN,d1,..,dN tensored with K[s].
@* The ideal I is the annihilator of f^s in D[s], that is the ideal, computed
@* by Sannfs-like procedure (usually called LD there).
@* Moreover, f is a polynomial in K[x1,..,xN] and qs is a polynomial in K[s].
RETURN:  int, 1 if qs is a factor of the global Bernstein polynomial of f and 0 otherwise
PURPOSE: check whether a univariate polynomial qs is a factor of the
@*  Bernstein-Sato polynomial of f without explicit knowledge of the latter.
NOTE:    If eng <>0, @code{std} is used for Groebner basis 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 checkFactor; shows examples
"
{

  // ASSUME too complicated, cannot check it.

  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  int ppl = printlevel-voice+2;
  def @R2 = basering;
  int N = nvars(@R2);
  int i;
  // we're in D_n[s], where the elim ord for s is set
  dbprint(ppl,"// -0-1- starting the procedure checkFactor");
  dbprint(ppl,"// -1-1- the ring @R2(_x,_Dx,s) is ready");
  dbprint(ppl-1, @R2);
  // create the ideal J = ann_D[s](f^s) + < f,q >
  ideal J = NF(I,std(F));
  // make leadcoeffs positive
  for (i=1; i<=ncols(J); i++)
  {
    if ( leadcoef(J[i])<0 )
    {
      J[i] = -J[i];
    }
  }
  J = J,F,q;
  // ------------ the ideal J is ready -----------
  dbprint(ppl,"// -1-2- starting the elimination of _x,_Dx in @R2");
  dbprint(ppl-1, J);
  ideal G = engine(J,eng);
  ideal K = nselect(G,1..N-1);
  kill J,G;
  dbprint(ppl,"// -1-3- _x,_Dx are eliminated");
  dbprint(ppl-1, K);
  //q is a factor of bs if and only if K = < q >
  //K = normalize(K);
  //q = normalize(q);
  //return( (K[1]==q) );
  return( NF(K[1],std(q))==0 );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),Dp;
  poly F = x^4+y^5+x*y^4;
  printlevel = 0;
  def A = Sannfs(F);
  setring A;
  poly F = imap(r,F);
  checkFactor(LD,F,20*s+31);     // -31/20 is not a root of bs
  checkFactor(LD,F,20*s+11);     // -11/20 is a root of bs
  checkFactor(LD,F,(20*s+11)^2); // the multiplicity of -11/20 is 1
}

proc varNum(string s)
"USAGE:  varNum(s);  string s
RETURN:  int
PURPOSE: returns the number of the variable with the name s
@*  among the variables of basering or 0 if there is no such variable
EXAMPLE: example varNum; shows examples
"
{
  int i;
  for (i=1; i<= nvars(basering); i++)
  {
    if ( string(var(i)) == s )
    {
      return(i);
    }
  }
  return(0);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring X = 0,(x,y1,t,z(0),z,tTa),dp;
  varNum("z");
  varNum("t");
  varNum("xyz");
}

static proc indAR(list L, int n)
"USAGE:  indAR(L,n);  list L, int n
RETURN:  list
PURPOSE: computes arrangement inductively, using L and
@* var(n) as the next variable
ASSUME: L has a structure of an arrangement
EXAMPLE: example indAR; shows examples
"
{
  if ( (n<2) || (n>nvars(basering)) )
  {
    ERROR("incorrect n");
  }
  int sl = size(L);
  list K;
  ideal tmp;
  poly @t = L[sl][1] + var(n); //1 elt
  K[sl+1] = ideal(@t);
  tmp = L[1]+var(n);
  K[1] = tmp; tmp = 0;
  int i,j,sI;
  ideal I;
  for(i=sl; i>=2; i--)
  {
    I = L[i-1]; sI = size(I);
    for(j=1; j<=sI; j++)
    {
      I[j] = I[j] + var(n);
    }
    tmp  = L[i],I;
    K[i] = tmp;
    I = 0; tmp = 0;
  }
  kill I; kill tmp;
  return(K);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z,t,v),dp;
  list L;
  L[1] = ideal(x);
  list K = indAR(L,2);
  K;
  list M = indAR(K,3);
  M;
  M = indAR(M,4);
  M;
}

proc isRational(number n)
"USAGE:  isRational(n); n number
RETURN:  int
PURPOSE: determine whether n is a rational number,
@* that is it does not contain parameters.
ASSUME: ground field is of characteristic 0
EXAMPLE: example indAR; shows examples
"
{
  if (char(basering) != 0)
  {
    ERROR("The ground field must be of characteristic 0!");
  }
  number dn = denominator(n);
  number nn = numerator(n);
  return( ((int(dn)==dn) && (int(nn)==nn)) );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = (0,a),(x,y),dp;
  number n1 = 11/73;
  isRational(n1);
  number n2 = (11*a+3)/72;
  isRational(n2);
}

proc bernsteinLift(ideal I, poly F, list #)
"USAGE:  bernsteinLift(I, F [,eng]);  I an ideal, F a poly, eng an optional int
RETURN:  list
PURPOSE: compute the (multiple of) Bernstein-Sato polynomial with lift-like method,
@* based on the output of Sannfs-like procedure
NOTE:  the output list contains the roots with multiplicities of the candidate
@* for being Bernstein-Sato polynomial of f.
@*  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 bernsteinLift; shows examples
"
{
  // assume: s is the last variable! check in the code
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "int" )
    {
      eng = int(#[1]);
    }
  }
  def @R2 = basering;
  int Nnew = nvars(@R2);
  int N = Nnew div 2;
  int ppl = printlevel-voice+2;
  // we're in D_n[s], where the elim ord for s is set
  // create D_n(s)
  // create the ordinary Weyl algebra and put the result into it,
  // keep: N, i,j,s, tmp, RL
  Nnew = Nnew - 1; // former 2*N;
  list L = 0;
  list Lord, tmp;
  intvec iv; int i;
  list RL = ringlist(basering);
  // if we work over alg. extension => problem!
  if (size(RL[1]) > 1)
  {
    ERROR("cannot work over algebraic field extension");
  }
  tmp[1] = RL[1]; // char
  tmp[2] = list("s");
  tmp[3] = list(list("lp",int(1)));
  tmp[4] = ideal(0);
  L[1] = tmp; // field
  tmp = 0;
  L[4] = RL[4];  // factor ideal

  // check whether vars have admissible names -> done earlier
  // list Name = RL[2]M
  // DName is defined earlier
  list NName; // = RL[2]; // skip the last var 's'
  for (i=1; i<=Nnew; i++)
  {
    NName[i] =  RL[2][i];
  }
  L[2] = NName;
  // (c, ) ordering:
  tmp[1]  = "c";
  iv = 0;
  tmp[2]  = iv;
  Lord[1] = tmp;
  tmp=0;
  // dp ordering;
  string s = "iv=";
  for (i=1; i<=Nnew; i++)
  {
    s = s+"1,";
  }
  s[size(s)] = ";";
  execute(s);
  tmp     = 0;
  tmp[1]  = "dp";  // string
  tmp[2]  = iv;  // intvec
  Lord[2] = tmp;
  kill s;
  tmp     = 0;
  L[3]    = Lord;
  // we are done with the list
  // Add: Plural part
  def @R4@ = ring(L);
  setring @R4@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  def @R4 = nc_algebra(1,@D);
  setring @R4;
  kill @R4@;
  dbprint(ppl,"// -3-1- the ring K(s)<x,dx> is ready");
  dbprint(ppl-1, @R4);
  // map things correctly, using names
  ideal J = imap(@R2, I), imap(@R2,F);
  module M;
  // make leadcoeffs positive
  for (i=1; i<= ncols(J); i++)
  {
    if (J[i]!=0)
    {
      M[i] = J[i]*gen(1) + gen(1+i);
    }
  }
  dbprint(ppl,"// -3-2- starting GB of the assoc. module M");
  M = engine(M,eng);
  dbprint(ppl,"// -3-3- finished GB of the assoc. module M");
  dbprint(ppl-1, M);
  // now look for (1) entry with 1st comp nonzero
  // determine whether there are several 1st comps nonzero
  module M2;
  for (i=1; i<= ncols(M); i++)
  {
    if (M[1,i]!=0)
    {
      M2 = M2, M[i];
    }
  }
  M2 = simplify(M2,2); // skip 0s
  if (ncols(M2) > 1)
  {
    dbprint(ppl,"// -*- more than 1 element with nonzero leading component");
    option(redSB);    option(redTail); // set them back?
    M2 = interred(M2);
    if (ncols(M2) > 1)
    {
      ERROR("more than one leading component after interred: assume violation!");
    }
    if (leadexp(M2[1]) != 0)
    {
      ERROR("nonconstant entry after interred: assume violation!");
    }
  }
  // now there's only one el-t with leadcomp<>0
  vector V = M2[1];
  number bcand = leadcoef(V[1]); // 1st component
  V[1]=0;
  number ct = content(V); // content of the cofactors
  poly CF = ct*V[ncols(J)]; // polynomial in K[s]<x,dx>, cofactor to F
  dbprint(ppl,"// -3-4- the cofactor candidate found");
  dbprint(ppl-1,CF);
  dbprint(ppl,"// -3-5- the entry as it is");
  dbprint(ppl-1,bcand);
  bcand = bcand*ct; // a product of both
  dbprint(ppl,"// -3-6- the content of the rest vector");
  dbprint(ppl-1,ct);
  ring @R3 = 0,s,dp;
  dbprint(ppl,"// -4-1- the ring @R3 i.e. K[s] is ready");
  poly bcand = imap(@R4,bcand);
  dbprint(ppl,"// -4-2- factorization");
  list P = factorize(bcand);          //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); // to get roots only
  setring @R2; // the ring the story started with
  ideal bs = imap(@R3,bs); // intvec m is global
  intvec mm = m; m = 0;
  kill @R3; // kills m as well....
  list @L = list(bs, mm);
  // look for (2) return the GB of syzygies?
  return(@L);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),Dp;
  poly F = x^3+y^3+z^3;
  printlevel = 0;
  def A = Sannfs(F);   setring A;
  LD;
  poly F = imap(r,F);
  list L  = bernsteinLift(LD,F); L;
  poly bs = fl2poly(L,"s"); bs; // the candidate for Bernstein-Sato polynomial
}

/// ****** EXAMPLES ************

/*

//static proc exCheckGenericity()
{
  LIB "control.lib";
  ring r = (0,a,b,c),x,dp;
  poly p = (x-a)*(x-b)*(x-c);
  def A = annfsBM(p);
  setring A;
  ideal J = slimgb(LD);
  matrix T = lift(LD,J);
  T = normalize(T);
  genericity(T);
  // Ann =x^3*Dx+3*x^2*t*Dt+(-a-b-c)*x^2*Dx+(-2*a-2*b-2*c)*x*t*Dt+3*x^2+(a*b+a*c+b*c)*x*Dx+(a*b+a*c+b*c)*t*Dt+(-2*a-2*b-2*c)*x+(-a*b*c)*Dx+(a*b+a*c+b*c)
  // genericity: g = a2-ab-ac+b2-bc+c2 =0
  // g = (a -(b+c)/2)^2 + (3/4)*(b-c)^2;
  // g ==0 <=> a=b=c
  // indeed, Ann = (x-a)^2*(x*Dx+3*t*Dt+(-a)*Dx+3)
  // --------------------------------------------
  // BUT a direct computation shows
  // when a=b=c,
  // Ann = x*Dx+3*t*Dt+(-a)*Dx+3
}

//static proc exOT_17()
{
  // Oaku-Takayama, p.208
  ring R = 0,(x,y),dp;
  poly F = x^3-y^2; // x^2+x*y+y^2;
  option(prot);
  option(mem);
  //  option(redSB);
  def A = annfsOT(F,0);
  setring A;
  LD;
  gkdim(LD); // a holonomic check
  //  poly F = x^3-y^2; // = x^7 - y^5; // x^3-y^4; // x^5 - y^4;
}

//static proc exOT_16()
{
  // Oaku-Takayama, p.208
  ring R = 0,(x),dp;
  poly F = x*(1-x);
  option(prot);
  option(mem);
  //  option(redSB);
  def A = annfsOT(F,0);
  setring A;
  LD;
  gkdim(LD); // a holonomic check
}

//static proc ex_bcheck()
{
  ring R = 0,(x,y),dp;
  poly F = x*y*(x+y);
  option(prot);
  option(mem);
  int eng = 0;
  //  option(redSB);
  def A = annfsOT(F,eng);
  setring A;
  LD;
}

//static proc ex_bcheck2()
{
  ring R = 0,(x,y),dp;
  poly F = x*y*(x+y);
  int eng = 0;
  def A = annfsBM(F,eng);
  setring A;
  LD;
}

//static proc ex_BMI()
{
  // a hard example
  ring r = 0,(x,y),Dp;
  poly F1 = (x2-y3)*(x3-y2);
  poly F2 = (x2-y3)*(xy4+y5+x4);
  ideal F = F1,F2;
  def A = annfsBMI(F);
  setring A;
  LD;
  BS;
}

//static proc ex2_BMI()
{
  // this example was believed to be intractable in 2005 by Gago-Vargas, Castro and Ucha
  ring r = 0,(x,y),Dp;
  option(prot);
  option(mem);
  ideal F = x2+y3,x3+y2;
  printlevel = 2;
  def A = annfsBMI(F);
  setring A;
  LD;
  BS;
}

//static proc ex_operatorBM()
{
  ring r = 0,(x,y,z,w),Dp;
  poly F = x^3+y^3+z^2*w;
  printlevel = 0;
  def A = operatorBM(F);
  setring A;
  F; // the original polynomial itself
  LD; // generic annihilator
  LD0; // annihilator
  bs; // normalized Bernstein poly
  BS; // root and multiplicities of the Bernstein poly
  PS; // the operator, s.t. PS*F^{s+1} = bs*F^s mod LD
  reduce(PS*F-bs,LD); // check the property of PS
}

*/