source: git/Singular/LIB/dmod.lib @ a2c2031

//////////////////////////////////////////////////////////////////////////////
version="$Id: dmod.lib,v 1.12 2007-01-23 15:03:30 Singular Exp $";
category="Noncommutative";
info="
LIBRARY: dmod.lib     Algorithms for algebraic D-modules
AUTHORS: Viktor Levandovskyy,     levandov@risc.uni-linz.ac.at
@*       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 ring R[1/F^s] for a natural number s.
@* In fact, the ring R[1/F^s] has a structure of a D(R)-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>.
@* Constructively, one needs to find a left ideal I = I(F^s) in D(R), such
@* that K[x_1,...,x_n,1/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 tensor over K with K[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))
@* - an operator in D(R)[s], denoted 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 (unpublished)
@* 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
@* - 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
@* - all the relevant data (LD, LD0, bs, PS) are computed by operatorBM

PROCEDURES:

annfs0(I,F [,eng]);    compute Ann F^s0 in D and Bernstein poly from the known Ann F^s in D[s]
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)
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)

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 and exceptional parametric constellations for a poly F with parametric coefficients (algorithm by Briancon and Maisonobe)

annfsBMI(F[,eng]);  compute Ann F^s and Bernstein ideal for a poly F=f1*..*fP (multivariate algorithm of Briancon-Maisonobe)


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
";

LIB "nctools.lib";
LIB "elim.lib";
LIB "qhmoduli.lib"; // for Max
LIB "gkdim.lib";
LIB "gmssing.lib";
LIB "control.lib";  // for genericity

static proc engine(ideal I, int i)
{
  /* std - slimgb mix */
  ideal J;
  if (i==0)
  {
    J = slimgb(I);
  }
  else
  {
    // without options -> strange! (ringlist?)
    option(redSB);
    option(redTail);
    J = std(I);
  }
  return(J);
}

// alternative code for SannfsBM, rename from annfsBM to ALTannfsBM
// is superfluos - will not be included in the official documentation
proc ALTannfsBM (poly F, list #)
"USAGE:  annfsBM(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.
@*       If printlevel=1, progress debug messages will be printed,
@*       if 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;
  ncalgebra(1,@D);
  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;
  }
  ncalgebra(1,@D);
  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.
@*       If printlevel=1, progress debug messages will be printed,
@*       if 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;
  ncalgebra(1,@D);
  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;
  }
  ncalgebra(1,@D);
  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 interger 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[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 is the needed D-mod structure,
@*       - 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.
@*       If printlevel=1, progress debug messages will be printed,
@*       if 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;
  ncalgebra(1,@D);
  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;
  }
  ncalgebra(1,@D);
  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 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 interger 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;
  }
  ncalgebra(1,@D);
  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 = x^3+y^3+z^3;
  printlevel = 0;
  def A = annfsBM(F);
  setring A;
  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.
@*       If printlevel=1, progress debug messages will be printed,
@*       if 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;
  ncalgebra(1,@D);
  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;
  }
  ncalgebra(1,@D);
  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 interger 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;
  }
  ncalgebra(1,@D);
  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,w),Dp;
  ring r = 0,(x,y,z),Dp;
  //  poly F = x^3+y^3+z^2*w;
  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; // 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
}

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.
@*       If printlevel=1, progress debug messages will be printed,
@*       if 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;
  ncalgebra(1,@D);
  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;
  }
  ncalgebra(1,@D);
  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[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;
  ncalgebra(1,@D);
  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;
  }
  ncalgebra(1,@D);
  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;
  }
  ncalgebra(1,@D);
  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[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 is the needed D-mod structure,
@*       - 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;
  ncalgebra(1,@D);
  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;
  }
  ncalgebra(1,@D);
  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;
  }
  ncalgebra(1,@D);
  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;
  }
  ncalgebra(1,@D);
  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[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 is the needed D-mod 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;
  ncalgebra(1,@D);
  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;
  }
  ncalgebra(1,@D);
  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;
  }
  ncalgebra(1,@D);
  dbprint(ppl,"//  -2-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[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 is the needed D-mod 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;
  ncalgebra(1,@D);
  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 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 @R2 = ring(L);
  setring @R2;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N; i++)
  {
    @D[i,N+i]=1;
  }
  ncalgebra(1,@D);
  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 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[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 is the needed D-mod 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;
  ncalgebra(1,@D);
  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;
  }
  ncalgebra(1,@D);
  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  = SannfsLOT(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[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 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 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 = x^3+y^3+z^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 procedures SannfsBM, SannfsOT or SannfsLOT
NOTE:    activate this ring with the @code{setring} command. In this ring,
@*       - the ideal LD 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
  int sP = minIntRoot(P[1],1);
  ideal bs; intvec m;             //the Bernstein polynomial is monic, so we are not interested in constants
  for (i=2; i<= size(P[1]); i++)  //we delete P[1][1] and P[2][1]
  {
    bs[i-1] = P[1][i];
    m[i-1] = P[2][i];
  }
  bs =  normalize(bs);
  bs = -subst(bs,s,0);
  dbprint(ppl,"// -2-3- minimal interger 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;
  }
  ncalgebra(1,@D);
  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);
  // alternatively, one can use SannfsOT or SannfsLOT
  setring A;
  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[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;
//   ncalgebra(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;
//   }
//   ncalgebra(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;
//   }
//   ncalgebra(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 minIntRoot; 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 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
    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
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);
}

static proc indAR(list L, int n)
"USAGE:  indAR(L,n);  list L, int n
RETURN:  list
PURPOSE: computes arrangement inductively, using L and varn(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;
}

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 = annfsOT(F,eng);
  setring A;
  LD;
}

static proc ex_BMI()
{
  ring r = 0,(x,y),Dp;
  poly F1 = (x2-y3)*(x3-y2);
  poly F2 = (x2-y3)*(xy4+y5+x4);
}