source: git/Singular/LIB/dmod.lib @ 3cec5e

//////////////////////////////////////////////////////////////////////////////
version="$Id: dmod.lib,v 1.35 2009-02-12 20:25:22 levandov Exp $";
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

THEORY: 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 F^s
@*      for a natural number 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 and
@*   the list of all roots of bs, which are rational, with their multiplicities 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 PS*F^(s+1) = bs*F^s holds in K[x_1,...,x_n,1/F^s].

@* We provide the following implementations:
@* OT) the classical Ann F^s algorithm from Oaku and Takayama (J. Pure
        Applied Math., 1999),
@* LOT) Levandovskyy's modification of the Oaku-Takayama algorithm (ISSAC 2007)
@* BM) the Ann F^s algorithm by Briancon and Maisonobe (Remarques sur
        l'ideal de Bernstein associe a des polynomes, preprint, 2002)

GUIDE:
@* - Ann F^s = I = I(F^s) = LD in D(R)[s] can be computed by SannfsBM, SannfsOT, SannfsLOT
@* - Ann^(1) F^s in D(R)[s] can be computed by Sannfslog
@* - global Bernstein polynomial bs resp. BS in K[s] can be computed by bernsteinBM
@* - Ann F^s0 = I(F^s0) = LD0 in D(R) can be computed by annfs0, annfsBM, annfsOT, annfsLOT, annfs2
@* - all the relevant data (LD, LD0, bs, PS) are computed by operatorBM

MAIN PROCEDURES:

annfs(F[,S,eng]);       compute Ann F^s0 in D and Bernstein poly for a poly F
annfspecial(I, F, m, n);  compute Ann F^n from Ann F^s for a poly F and a number n
Sannfs(F[,S,eng]);      compute Ann F^s in D[s] for a poly F
Sannfslog(F[,eng]);     compute Ann^(1) F^s in D[s] for a poly F
bernsteinBM(F[,eng]);   compute global Bernstein poly for a poly F (algorithm of Briancon-Maisonobe)
operatorBM(F[,eng]);    compute Ann F^s, Ann F^s0, BS and PS for a poly F (algorithm of Briancon-Maisonobe)
annfsParamBM(F[,eng]);  compute the generic Ann F^s (algorithm by Briancon and Maisonobe) and exceptional parametric constellations for a poly F with parametric coefficients
annfsBMI(F[,eng]);      compute Ann F^s and Bernstein ideal for a poly F=f1*..*fP (multivariate algorithm of Briancon-Maisonobe)
checkRoot(F,a[,S,eng]); check if a given rational is a root of the global Bernstein polynomial of F and compute its multiplicity
annfsBM(F[,eng]);          compute Ann F^s0 in D and Bernstein poly for a poly F (algorithm of Briancon-Maisonobe)
annfsLOT(F[,eng]);         compute Ann F^s0 in D and Bernstein poly for a poly F (Levandovskyy modification of the Oaku-Takayama algorithm)
annfsOT(F[,eng]);          compute Ann F^s0 in D and Bernstein poly for a poly F (algorithm of Oaku-Takayama)
SannfsBM(F[,eng]);         compute Ann F^s in D[s] for a poly F (algorithm of Briancon-Maisonobe)
SannfsBFCT(F[,eng]);      compute Ann F^s in D[s] for a poly F (algorithm of Briancon-Maisonobe, other output ordering)
SannfsLOT(F[,eng]);        compute Ann F^s in D[s] for a poly F (Levandovskyy modification of the Oaku-Takayama algorithm)
SannfsOT(F[,eng]);         compute Ann F^s in D[s] for a poly F (algorithm of Oaku-Takayama)
annfs0(I,F [,eng]);          compute Ann F^s0 in D and Bernstein poly from the known Ann F^s in D[s]
annfs2(I,F [,eng]);           compute Ann F^s0 in D and Bernstein poly from the known Ann F^s in D[s] by using a trick of Noro
annfsRB(I,F [,eng]);          compute Ann F^s0 in D and Bernstein poly from the known Ann F^s in D[s] by reduceB strategy of Macaulay
checkRoot1(I,F,a[,eng]);   check whether a rational is a root of the global Bernstein polynomial of F from the known Ann F^s in D[s]
checkRoot2(I,F,a[,eng]);   check whether a rational is a root of the global Bernstein polynomial of F and compute its multiplicity from the known Ann F^s in D[s]
checkFactor(I,F,q[,eng]); check whether a polynomial q in K[s] is a factor of the global Bernstein polynomial of F from the known Ann F^s in D[s]

AUXILIARY PROCEDURES:

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


SEE ALSO: gmssing_lib, bfun_lib, dmodapp_lib
";

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

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;
  "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.
@*       The value of a string S can be
@*       '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 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 (which is a Groebner basis) is the needed D-module structure,
@*       - the list  BS is the list of roots and multiplicities of a Bernstein polynomial of f.
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;
  LD;
  BS;
}


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 Oaku and Takayama,
@*       '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, where D is the Weyl algebra
NOTE:    activate this ring with the @code{setring} command.
@*       In the 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 this 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 its 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 this 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);
  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;
  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
NOTE:    activate this ring with the @code{setring} command. In this ring,
@*       - the ideal LD (which is a Groebner basis) is the annihilator of f^s,
@*       - 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.
@*       annfs2 uses the shorter form of expressions in the variable s (the idea of Noro).
DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
@*       if @code{printlevel}>=2, all the debug messages will be printed.
EXAMPLE: example annfs2; shows examples
"
{
  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/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 Macaulay's reduceB strategy, that is add
// not F but <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 this 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 follows the 'reduceB' strategy, used in Macaulay2.
DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
@*       if @code{printlevel}>=2, all the debug messages will be printed.
EXAMPLE: example annfsRB; shows examples
"
{
  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)/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 wrt 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/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; // 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, according to the algorithm by Briancon and Maisonobe
NOTE:    activate this ring with the @code{setring} command. In this 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);
  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;
  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;

// TODO: 1 has to appear in the 1nd column of transp. matrix
// this does not happen automatically
// for this, do special modulo with emphasis on the 1comp (c,<)
// of the transp matrix, there must appear 1 in the GB
// then use lift to get the expression...
// need: (c,<) ordering for such comp's

/*
proc operatorModulo(poly F, ideal I, poly b)
"USAGE:  operatorModulo(f,I,b);  f a poly, I an ideal, b a poly
RETURN:  poly
PURPOSE: compute the B-operator from the poly f, ideal I = Ann f^s and Bernstein-Sato
polynomial b using modulo i.e. kernel of module homomorphism
NOTE:  The computations take place in the ring, similar to the one returned by Sannfs procedure.
@*       If printlevel=1, progress debug messages will be printed,
@*       if printlevel>=2, all the debug messages will be printed.
EXAMPLE: example operatorModulo; shows examples
"
{
  // with hom_kernel;
  matrix AA[1][2] = F,-b;
  //  matrix M = matrix(subst(I,s,s+1)*freemodule(2)); // ann f^{s+1}
  matrix N = matrix(I); // ann f^s
  //  matrix K = hom_kernel(AA,M,N);
  matrix K = modulo(AA,N);
  K = transpose(K);
  print(K);
  ideal J;
  int i;
  poly t; number n;
  for(i=1; i<=nrows(K); i++)
  {
    if (K[i,2]!=0)
    {
      if ( leadmonom(K[i,2]) == 1)
      {
        t = K[i,1];
        n = leadcoef(K[i,2]);
        t = t/n;
        //        J = J, K[i][2];
        break;
      }
    }
  }
  ideal J = groebner(subst(I,s,s+1)); // for NF
  t = NF(t,J);
  "candidate:"; t;
  J = subst(J,s,s-1);
  // test:
  if ( NF(t*F-b,J) !=0)
  {
    "Problem: PS does not work on F";
  }
  return(t);
}
example
{
  "EXAMPLE:"; echo = 2;
  //  ring r = 0,(x,y,z),Dp;
  //   poly F = x^3+y^3+z^3;
  LIB "dmod.lib"; option(prot); option(mem);
  ring r = 0,(x,y),Dp;
  //  poly F = x^3+y^3+x*y^2;
  poly F = x^4 + y^4 + x*y^4;
  def A = Sannfs(F); // here we get LD = ann f^s
  setring A;
  poly F = imap(r,F);
  def B = annfs0(LD,F); // to obtain BS polynomial
  list BS = imap(B,BS);
  poly b = fl2poly(BS,"s");
  LD = groebner(LD);
  poly PS = operatorModulo(F,LD,b);
  PS; // the operator, s.t. PS*F^{s+1} = bs*F^s mod LD
  reduce(PS*F-bs,LD); // check the property of PS
}
*/



proc 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, according to the algorithm by Briancon and Maisonobe
NOTE:    activate this 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 the list of the special parameters.
@*       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;
  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;
  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]);  F an ideal, eng an optional int
RETURN:  ring
PURPOSE: compute the D-module structure of basering[1/f]*f^s where f = F[1]*..*F[P],
according to the algorithm by Briancon and Maisonobe.
NOTE:    activate this ring with the @code{setring} command. In this ring,
@*       - the ideal LD is the needed D-mod structure,
@*       - the list BS is the Bernstein 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 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;
  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 P = size(F);  //if F has some generators which are zero, int P = ncols(I);
  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 ----------
  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); ?
  poly f=1;
  for (j=1; j<=P; j++)
  {
    f = f*F[j];
  }
  K = K,f;       // to compute B (Bernstein-Sato ideal)
  //j=2;         // for example
  //K = K,F[j];  // to compute Bj (see "On the computation of Bernstein-Sato ideals"; Castro, Ucha)
  //K = K,F;     // to compute Bsigma (see "On the computation of Bernstein-Sato ideals"; 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 (size(K3)==1)
  {
    poly p = K3[1];
    dbprint(ppl,"// -3-2- 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: 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);
  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;
  return(@R4);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),Dp;
  ideal F = x,y,x+y;
  printlevel = 0;
  def A = annfsBMI(F);
  setring A;
  LD;
  BS;
}

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 this 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);
  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;
  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], where D is the Weyl algebra
NOTE:    activate this 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 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], where D is the Weyl algebra
NOTE:    activate this 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 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;
}

proc SannfsBFCT(poly F, list #)
"USAGE:  SannfsBFCT(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 this 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.
@*       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 SannfsBFCT; shows examples
"
{
  // DEBUG INFO: ordering on the output ring = dp, engine on the sum of ideals is used
  // hence it's the situation for slimgb

  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);
  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);
}

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 this 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 the ring D[s], where D is the Weyl algebra
NOTE:    activate this 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 this 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 this 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 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 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 this 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.
@*       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/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;
  // 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 this 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 poly 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.
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 a rational number n
ASSUME:  the basering contains 's' as a variable
NOTE:    We assume that the basering is D[s],
@*          ideal I is the Ann F^s in D[s] (obtained with e.g. SannfsBM, SannfsLOT, SannfsOT)
@*          integer 'mir' is the minimal integer root of the Bernstein polynomial of F
@*          and the number n is rational.
@*       We compute the real annihilator for any rational value of n (both generic and exceptional).
@*       The implementation goes along the lines of Saito-Sturmfels-Takayama, Alg. 5.3.15
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
"
{
  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 = 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!";
  }
  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
  minIntRoot(I1,0);
  poly  f2  = x2-y3;
  ideal I2  = bernstein(f2)[1];
  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-poly of f1=x*y*(x+y)
  ideal I3 = factorize(f3,1);
  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];
  minIntRoot(I,0);
}

proc isHolonomic(def M)
"USAGE:  isHolonomic(M); M an ideal/module/matrix
RETURN:  int, 1 if M is holonomic and 0 otherwise
ASSUME: basering is a Weyl algebra
PURPOSE: check the modules for the property of holonomy
NOTE:    M is holonomic, if 2*dim(M) = dim(R), where R is a
ground ring; dim stays for Gelfand-Kirillov dimension
EXAMPLE: example isHolonomic; shows examples
"
{
  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 = std(M);
  }
  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 this ring with the @code{setring} command and find the
         curve as a polynomial RC
@*   a Reiffen curve is defined as F = x^p + y^q + xy^{q-1}, q >= p+1 >= 5
ASSUME: q >= p+1 >= 5. Otherwise an error message is returned
EXAMPLE: example reiffen; shows examples
"
{
// a Reiffen curve is defined as
// F = x^p + y^q +x*y^{q-1}, q \geq p+1 \geq 5

  if ( (p<4) || (q<5) || (q-p<1) )
  {
    ERROR("Some of conditions p>=4, q>=5 or q>=p+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 ring
ASSUME: basering is present
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]);  f a poly, alpha a number, S a string , eng an optional int
RETURN:  int
PURPOSE: check whether a rational is a root of the global Bernstein polynomial of f (and compute its multiplicity)
with the algorithm given in S and with the Groebner basis engine given in eng
NOTE:    The annihilator of f^s in D[s] is needed and it is computed according to the algorithm by Briancon and Maisonobe
@*       The value of a string S can be
@*       'alg1' (default) - for the algorithm 1 of J. Martin-Morales (unpublished)
@*       'alg2' - for the algorithm 2 of J. Martin-Morales (unpublished)
@*       The output of type int is:
@*       'alg1': 1 if -alpha is a root of the global Bernstein polynomial and 0 otherwise
@*       'alg2':  the multiplicity of -alpha as a root of the global Bernstein polynomial of f.
@*                (If -alpha is not a root, the output is 0)
@*       If eng <>0, @code{std} is used for Groebner basis computations,
@*       otherwise (and by default) @code{slimgb} is used.
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 of J. Martín-Morales");
       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;
  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]);  I an ideal, f a poly, alpha a number, eng an optional int
ASSUME:  I is the annihilator of f^s in D[s], f is a polynomial in K[_x]
RETURN:  int, 1 if -alpha is a root of the global Bernstein polynomial of f and 0 otherwise
PURPOSE: check whether a rational is a root of the global Bernstein 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
"
{
  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 / 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,alpha [,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], f is a polynomial in K[_x]
RETURN:  int, the multiplicity of -alpha as a root of the global Bernstein polynomial of f. If -alpha is not a root, the output is 0
PURPOSE: 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]
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 checkRoot2; shows examples
"
{
  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) / 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:  I is the output of Sannfs, SannfsBM, SannfsLOT or SannfsOT,
         f is a polynomial in K[_x], 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 rational is a root of the global Bernstein polynomial of f and compute its multiplicity from the known Ann F^s in D[s]
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
"
{
  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 iff 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,z(0),tTa),dp;
  varnum("z(0)");
  varnum("tTa");
  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;
}

/// ****** 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
}

*/