///////////////////////////////////////////////////////////////////////////////
version="$Id: dmod.lib,v 1.8 2006-07-18 15:48:12 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 a 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 provide two implementations:
@* 1) the classical Ann F^s algorithm from Oaku and Takayama (J. Pure
        Applied Math., 1999) and
@* 2) the newer Ann F^s algorithm by Briancon and Maisonobe (Remarques sur
        l'ideal de Bernstein associe a des polynomes, preprint, 2002).

PROCEDURES:
annfsOT(F[,eng]); compute Ann F^s for a poly F (algorithm of Oaku-Takayama)
annfsBM(F[,eng]); compute Ann F^s for a poly F (algorithm of Briancon-Maisonobe)
minIntRoot(P,fact); minimal integer root of a maximal ideal P
reiffen(p,q);     create a poly, describing a Reiffen curve
arrange(p);       create a poly, describing a generic hyperplane arrangement
isHolonomic(M);   check whether a module is holonomic
convloc(L);       replace global orderings with local in the ringlist L
";

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

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

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 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 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);
  dbprint(ppl,"// -1-3- t is eliminated");
  dbprint(ppl-1, K);  // K is without t
  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);
  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);
  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 is ready");
  ideal K3 = imap(@R2,K2);
  poly p = K3[1];
  dbprint(ppl,"// -3-2- factorization");
  ideal P = factorize(p,1);  //without constants and multiplicities
  //  "--------- b-function factorizes into ---------";   P;
  int sP = minIntRoot(P,1);
  dbprint(ppl,"// -3-3- minimal interger 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;
  //TODO: sort BS!
  // --------- substitute s found in the ideal ---------
  // --------- going back to @R and substitute ---------
  setring @R;
  ideal K2 = subst(K,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]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);
  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;
  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 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;
  ideal 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^2+y^3+z^5;
  printlevel = 0;
  def A = annfsBM(F);
  setring A;
  LD;
  print(matrix(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 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 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
  //  "------ b-function factorizes into ----------";  P;
  int sP  = minIntRoot(P, 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;
  // 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);
  export BS;
  kill @R4;
  ideal LD = K5;
  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;
  print(matrix(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];
  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); // 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 generic 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;
  LIB "gkdim.lib";
  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;
  LIB "gkdim.lib";
  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);
  //  def A = annfsgms(F,0);
  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);
  //  def A = annfsgms(F,0);
  setring A;
  LD;
}