////////////////////////////////////////////////////////////////////////
version="$Id$";
category="Noncommutative";
info="
LIBRARY: dmodvar.lib     Algebraic D-modules for varieties

AUTHORS: Daniel Andres,          daniel.andres@math.rwth-aachen.de
       Viktor Levandovskyy,    levandov@math.rwth-aachen.de
       Jorge Martin-Morales,   jorge@unizar.es

THEORY: Let K be a field of characteristic 0. Given a polynomial ring R = K[x_1,...,x_n] and 
@* a set of polynomial f_1,..., f_r in R, define F = f_1 * ... * f_r and F^s:=f_1^s_1*...*f_r^s_r 
@* for symbolic discrete (that is shiftable) variables s_1,..., s_r.
@* The module R[1/F]*F^s has a structure of a D<S>-module, where D<S> := D(R) tensored with S over K, 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>
@*   - S is the universal enveloping algebra of gl_r, generated by s_{ij}, where s_{ii}=s_i.
@* One is interested in the following data:
@*   - the left ideal Ann F^s in D<S>, usually denoted by LD in the output
@*   - global Bernstein polynomial in one variable s = s_1 + ...+ s_r, denoted by bs,
@*   - its minimal integer root s0, the list of all roots of bs, which are known
@*     to be negative rational numbers, with their multiplicities, which is denoted by BS
@*   - an r-tuple of operators in D<S>, denoted by PS, such that the functional equality
@*     sum(k=1 to k=r) P_k*f_k*F^s = bs*F^s holds in R[1/F]*F^s.

REFERENCES:
  (BMS06) Budur, Mustata, Saito: Bernstein-Sato polynomials of arbitrary varieties (2006).
  (ALM09) Andres, Levandovskyy, Martin-Morales : Principal Intersection and Bernstein-Sato Polynomial of an Affine Variety (2009).

MAIN PROCEDURES:
bfctVarIn(F[,L]);     compute the roots of the Bernstein-Sato polynomial b(s) of the variety V(F) using initial ideal approach
bfctVarAnn(F[,L]);  compute the roots of the Bernstein-Sato polynomial b(s) of the variety V(F) using Sannfs approach
SannfsVar(F[,O,e]); compute the annihilator of F^s in the ring D<S>

AUXILIARY PROCEDURES:
makeIF(F[,ORD]);    create the Malgrange ideal, associated with F = F[1],..,F[P]

SEE ALSO: bfun_lib, dmod_lib, dmodapp_lib, gmssing_lib

KEYWORDS: D-module; D-module structure; Bernstein-Sato polynomial for variety; global Bernstein-Sato polynomial for variety;
Weyl algebra; parametric annihilator for variety; Budur-Mustata-Saito approach; initial ideal approach;
";

// Static procs:
//coDim(I);           compute the codimension of the leading ideal of I
// dmodvarAssumeViolation()
// ORDstr2list (ORD, NN)
// smallGenCoDim(I,k)


LIB "bfun.lib";    // for pIntersect
LIB "dmodapp.lib"; // for isCommutative etc.


///////////////////////////////////////////////////////////////////////////////

// testing for consistency of the library:
proc testdmodvarlib ()
{
  "AUXILIARY PROCEDURES:";
  example makeIF;
  "MAIN PROCEDURES:";
  example bfctVarIn;
  example bfctVarAnn;
  example SannfsVar;
}

//   example coDim;

///////////////////////////////////////////////////////////////////////////////

static proc dmodvarAssumeViolation()
{
  // returns Boolean : yes/no [for assume violation]
  // char K = 0
  // no qring
  if (  (size(ideal(basering)) >0) || (char(basering) >0) )
  {
    //    "ERROR: no qring is allowed";
    return(1);
  }
  return(0);
}

// da: in smallGenCoDim(), rewritten using mstd business
static proc coDim (ideal I)
"USAGE:  coDim (I);  I an ideal
RETURN:  int
PURPOSE: computes the codimension of the ideal generated by the leading monomials
@*       of the given generators of the ideal. This is also the codimension of
@*       the ideal if it is represented by a standard basis.
NOTE:    The codimension of an ideal I means the number of variables minus the
@*       Krull dimension of the basering modulo I.
EXAMPLE: example SannfsVar; shows examples
"
{
  int n = nvars(basering);
  int d = dim(I); // to insert: check whether I is in GB
  return(n-d);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring R = 0,(x,y,z),Dp;
  ideal I = x^2+y^3, z;
  coDim(std(I));
}

static proc ORDstr2list (string ORD, int NN)
{
  /* convert an ordering defined in NN variables in the  */
  /* string form into the same ordering in the list form */
  string st;
  st = "ring @Z = 0,z(1.." + string(NN) + "),";
  st = st + ORD + ";";
  execute(st); kill st;
  list L = ringlist(@Z)[3];
  kill @Z;
  return(L);
}

proc SannfsVar (ideal F, list #)
"USAGE:  SannfsVar(F [,ORD,eng]);  F an ideal, ORD an optional string, eng an optional int
RETURN:  ring
PURPOSE: compute the D<S>-module structure of D<S>*f^s where f = F[1]*..*F[P]
and D<S> is the Weyl algebra D tensored with K<S>=U(gl_P), according to the
generalized algorithm by Briancon and Maisonobe for affine varieties.
NOTE:    activate this ring with the @code{setring} command.
@*       In the ring D<S>, the ideal LD is the needed D<S>-module structure.
@*       The value of ORD must be an elimination ordering in D<Dt,S> for Dt
@*       written in the string form, otherwise the result may have no meaning.
@*       By default ORD = '(a(1..(P)..1),a(1..(P+P^2)..1),dp)'.
@*       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 SannfsVar; shows examples
"
{
  if (dmodvarAssumeViolation())
  {
    ERROR("Basering is inappropriate: characteristic>0 or qring present");
  }
  if (!isCommutative())
  {
    ERROR("Basering must be commutative");
  }
  def save = basering;
  int N = nvars(basering);
  int P = ncols(F);  //ncols better than size, since F[i] could be zero
  // P is needed for default ORD
  int i,j,k,l;
  // st = "(a(1..(P)..1),a(1..(P+P^2)..1),dp)";
  string st = "(a(" + string(1:P);
  st = st + "),a(" + string(1:(P+P^2));
  st = st + "),dp)";
  // default values
  string ORD = st;
  int eng = 0;
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "string" )
    {
      ORD = string(#[1]);
      // second arg
      if (size(#)>1)
      {
        // exists 2nd arg
        if ( typeof(#[2]) == "int" )
        {
          // the case: given ORD, given engine
          eng = int(#[2]);
        }
        else
        {
          eng = 0;
        }
      }
      else
      {
        // no second arg
        eng = 0;
      }
    }
    else
    {
      if ( typeof(#[1]) == "int" )
      {
        // the case: default ORD, engine given
        eng = int(#[1]);
        // ORD = "(a(1..(P)..1),a(1..(P+P^2)..1),dp)";  //is already set
      }
      else
      {
        // incorr. 1st arg
        ORD = string(st);
      }
    }
  }
  // size(#)=0, i.e. there is no elimination ordering and no engine given
  // eng = 0; ORD = "(a(1..(P)..1),a(1..(P^2)..1),dp)";  //are already set
  int ppl = printlevel-voice+2;
  // returns a list with a ring and an ideal LD in it
  // save, N, P and the indices are already defined
  int Nnew = 2*N+P+P^2;
  list RL = ringlist(basering);
  list L;
  L[1] = RL[1];  //char
  L[4] = RL[4];  //char, minpoly
  // check whether vars have admissible names
  list Name  = RL[2];
  list RName;
  // (i,j) <--> (i-1)*p+j
  for (i=1; i<=P; i++)
  {
    RName[i] = "Dt("+string(i)+")";
    for (j=1; j<=P; j++)
    {
      st = "s("+string(i)+")("+string(j)+")";
      RName[P+(i-1)*P+j] = st;
    }
  }
  for(i=1; i<=N; i++)
  {
    for(j=1; j<=size(RName); j++)
    {
      if (Name[i] == RName[j])
      {
        ERROR("Variable names should not include Dt(i), s(i)(j)");
      }
    }
  }
  // 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 (a(1..(P)..1),a(1..(P+P^2)..1),dp);
  //export Nnew;
  L[3] = ORDstr2list(ORD,Nnew);
  // we are done with the list
  def @R@ = ring(L);
  setring @R@;
  matrix @D[Nnew][Nnew];
  // kronecker(i,j) equals (i==j)
  // (i,j) <--> (i-1)*p+j
  for (i=1; i<=P; i++)
  {
    for (j=1; j<=P; j++)
    {
      for (k=1; k<=P; k++)
      {
        //[sij,Dtk] = djk*Dti
        @D[k,P+(i-1)*P+j] = (j==k)*Dt(i);
        for (l=1; l<=P; l++)
        {
          if ( (i-k)*P < l-j )
          {
            //[sij,skl] = djk*sil - dil*skj
            @D[P+(i-1)*P+j,P+(k-1)*P+l] = -(j==k)*s(i)(l) + (i==l)*s(k)(j);
          }
        }
      }
    }
  }
  for (i=1; i<=N; i++)
  {
    //[Dx,x]=1
    @D[P+P^2+i,P+P^2+N+i] = 1;
  }
  def @R = nc_algebra(1,@D);
  setring @R;
  //@R@ will be used further
  dbprint(ppl,"// -1-1- the ring @R(_Dt,_s,_x,_Dx) is ready");
  dbprint(ppl-1, @R);
  // create the ideal I
  // (i,j) <--> (i-1)*p+j
  ideal  F = imap(save,F);
  ideal I;
  for (i=1; i<=P; i++)
  {
    for (j=1; j<=P; j++)
    {
      I[(i-1)*P+j] = Dt(i)*F[j] + s(i)(j);
    }
  }
  poly p,q;
  for (i=1; i<=N; i++)
  {
    p=0;
    for (j=1; j<=P; j++)
    {
      q = Dt(j);
      q = q*diff(F[j],var(P+P^2+i));
      p = p + q;
    }
    I = I, p + var(P+P^2+N+i);
  }
  // -------- the ideal I is ready ----------
  dbprint(ppl,"// -1-2- starting the elimination of "+string(Dt(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 Dt(i) are eliminated");
  dbprint(ppl-1, K);  //K is without Dt(i)
  // ----------- the ring @R2(_s,_x,_Dx) ------------
  //come back to the ring save, recover L and remove all Dt(i)
  //L[1],L[4] do not change
  setring save;
  list Lord, tmp;
  // variables
  tmp = L[2];
  Lord = tmp[P+1..Nnew];
  L[2] = Lord;
  // ordering
  // st = "(a(1..(P^2)..1),dp)";
  st = "(a(" + string(1:P^2);
  st = st + "),dp)";
  tmp = ORDstr2list(st,Nnew-P);
  L[3] = tmp;
  def @R2@ = ring(L);
  kill L;
  // we are done with the list
  setring @R2@;
  matrix tmpM,LordM;
  // non-commutative relations
  intvec iv = P+1..Nnew;
  tmpM = imap(@R@,@D);
  kill @R@;
  LordM = submat(tmpM,iv,iv);
  matrix @D2 = LordM;
  def @R2 = nc_algebra(1,@D2);
  setring @R2;
  kill @R2@;
  dbprint(ppl,"// -2-1- the ring @R(_s,_x,_Dx) is ready");
  dbprint(ppl-1, @R2);
  ideal K = imap(@R,K);
  kill @R;
  dbprint(ppl,"// -2-2- starting cosmetic Groebner basis computation");
  dbprint(ppl-1, K);
  K = engine(K,eng);
  dbprint(ppl,"// -2-3- the cosmetic Groebner basis has been computed");
  dbprint(ppl-1,K);
  ideal LD = K;
  attrib(LD,"isSB",1);
  export LD;
  return(@R2);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring R = 0,(x,y),Dp;
  ideal F = x^3, y^5;
  //ORD = "(a(1,1),a(1,1,1,1,1,1),dp)";
  //eng = 0;
  def A = SannfsVar(F);
  setring A;
  LD;
}

proc bfctVarAnn (ideal F, list #)
"USAGE:  bfctVarAnn(F[,gid,eng]); F an ideal, gid,eng optional ints
RETURN:  list of an ideal and an intvec
PURPOSE: computes the roots of the Bernstein-Sato polynomial and their multiplicities 
@*       for an affine algebraic variety defined by F = F[1],..,F[r].
ASSUME:  The basering is a commutative polynomial ring in char 0.
BACKGROUND: In this proc, the annihilator of f^s in D[s] is computed and then a
@*       system of linear equations is solved by linear reductions in order to
@*       find the minimal polynomial of S = s(1)(1) + ... + s(P)(P)
NOTE:    In the output list, the ideal contains all the roots and the intvec their multiplicities.
@*       If gid<>0, the ideal is used as given. Otherwise, and by default, a
@*       heuristically better suited generating set is used.
@*       If eng<>0, @code{std} is used for GB computations,
@*       otherwise, and by default, @code{slimgb} is used.
DISPLAY: If printlevel=1, progress debug messages will be printed,
@*       if printlevel=2, all the debug messages will be printed.
COMPUTATIONAL REMARK: The time of computation can be very different depending
@*       on the chosen generators of F, although the result is always the same.
EXAMPLE: example bfctVarAnn; shows examples
"
{
  if (dmodvarAssumeViolation())
  {
    ERROR("Basering is inappropriate: characteristic>0 or qring present");
  }
  if (!isCommutative())
  {
    ERROR("Basering must be commutative");
  }
  int gid = 0; // default
  int eng = 0; // default
  if (size(#)>0)
  {
    if (typeof(#[1])=="int" || typeof(#[1])=="number")
    {
      gid = int(#[1]);
    }
    if (size(#)>1)
    {
      if (typeof(#[2])=="int" || typeof(#[2])=="number")
      {
        eng = int(#[2]);
      }
    }
  }
  def save = basering;
  int ppl = printlevel - voice + 2;
  printlevel = printlevel+1;
  list L = smallGenCoDim(F,gid);
  F = L[1];
  int cd = L[2];
  kill L;
  def @R2 = SannfsVar(F,eng);
  printlevel = printlevel-1;
  setring @R2;
  // we are in D[s] and LD is a std of SannfsVar(F)
  ideal F = imap(save,F);
  ideal GF = std(F);
  ideal J = NF(LD,GF);
  J = J, F;
  dbprint(ppl,"// -3-1- starting Groebner basis of ann F^s + F ");
  dbprint(ppl-1,J);
  ideal K = engine(J,eng);
  dbprint(ppl,"// -3-2- finished Groebner basis of ann F^s + F ");
  dbprint(ppl-1,K);
  poly S;
  int i;
  for (i=1; i<=size(F); i++)
  {
    S = S + s(i)(i);
  }
  dbprint(ppl,"// -4-1- computing the minimal polynomial of S");
  //dbprint(ppl-1,"S = "+string(S));
  module M = pIntersect(S,K);
  dbprint(ppl,"// -4-2- the minimal polynomial has been computed");
  //dbprint(ppl-1,M);
  ring @R3 = 0,s,dp;
  dbprint(ppl,"// -5-1- the ring @R3(s) is ready");
  dbprint(ppl-1,@R3);
  ideal M = imap(@R2,M);
  //kill @R2;
  poly p;
  for (i=1; i<=size(M); i++)
  {
    p = p + M[i]*s^(i-1);
  }
  dbprint(ppl,"// -5-2- factorization of the minimal polynomial");
  list P = factorize(p);          //with constants and multiplicities
  dbprint(ppl-1,P);               //the Bernstein polynomial is monic,
  ideal bs; intvec m;             //so we are not interested in constants
  for (i=2; i<=ncols(P[1]); i++)   //and that is why we delete P[1][1] and P[2][1]
  {
    bs[i-1] = P[1][i];
    m[i-1]  = P[2][i];
  }
  // convert factors to a list of their roots and multiplicities
  bs =  normalize(bs);
  bs = -subst(bs,s,0);
  setring save;
//   ideal GF = imap(@R2,GF);
//   attrib(GF,"isSB",1);
  kill @R2;
  dbprint(ppl,"// -5-3- codimension of the variety");
//   int cd = coDim(GF);
  dbprint(ppl-1,cd);
  ideal bs = imap(@R3,bs);
  dbprint(ppl,"// -5-4- shifting BS(s)=minpoly(s-codim+1)");
  for (i=1; i<=ncols(bs); i++)
  {
    bs[i] = bs[i] + cd - 1;
  }
  kill @R3;
  list BS = bs,m;
  return(BS);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring R = 0,(x,y,z),Dp;
  ideal F = x^2+y^3, z;
  bfctVarAnn(F);
}

proc makeIF (ideal F, list #)
"USAGE:  makeIF(F [,ORD]);  F an ideal, ORD an optional string
RETURN:  ring
PURPOSE: create the ideal by Malgrange associated with F = F[1],..,F[P].
NOTE:    activate this ring with the @code{setring} command. In this ring,
@*       - the ideal IF is the ideal by Malgrange corresponding to F.
@*       The value of ORD must be an arbitrary ordering in K<_t,_x,_Dt,_Dx>
@*       written in the string form. By default ORD = 'dp'.
@*       If printlevel=1, progress debug messages will be printed,
@*       if printlevel>=2, all the debug messages will be printed.
EXAMPLE: example makeIF; shows examples
"
{
  string ORD = "dp";
  if ( size(#)>0 )
  {
    if ( typeof(#[1]) == "string" )
    {
      ORD = string(#[1]);
    }
  }
  int ppl = printlevel-voice+2;
  def save = basering;
  int N = nvars(save);
  int P = ncols(F);
  int Nnew = 2*P+2*N;
  int i,j;
  string st;
  list RL = ringlist(save);
  list L,Lord;
  list tmp;
  intvec iv;
  L[1] = RL[1];
  L[4] = RL[4];
  //check whether vars have admissible names
  list Name = RL[2];
  list TName, DTName;
  for (i=1; i<=P; i++)
  {
    TName[i] = "t("+string(i)+")";
    DTName[i] = "Dt("+string(i)+")";
  }
  for (i=1; i<=N; i++)
  {
    for (j=1; j<=P; j++)
    {
      if (Name[i] == TName[j])
      {
        ERROR("Variable names should not include t(i)");
      }
    }
  }
  //now, create the names for new vars
  list DName;
  for (i=1; i<=N; i++)
  {
    DName[i] = "D"+Name[i];  //concat
  }
  list NName = TName + Name + DTName + DName;
  L[2]   = NName;
  // Name, Dname will be used further
  kill NName, TName, Name, DTName, DName;
  // ORD already set, default ord dp;
  L[3] = ORDstr2list(ORD,Nnew);
  // we are done with the list
  def @R@ = ring(L);
  setring @R@;
  matrix @D[Nnew][Nnew];
  for (i=1; i<=N+P; i++)
  {
    @D[i,i+N+P]=1;
  }
  def @R = nc_algebra(1,@D);
  setring @R;
  kill @R@;
  //dbprint(ppl,"// -1-1- the ring @R(_t,_x,_Dt,_Dx) is ready");
  //dbprint(ppl-1, @R);
  // create the ideal I
  ideal  F = imap(save,F);
  ideal I;
  for (j=1; j<=P; j++)
  {
    I[j] = t(j) - F[j];
  }
  poly p,q;
  for (i=1; i<=N; i++)
  {
    p=0;
    for (j=1; j<=P; j++)
    {
      q = Dt(j);
      q = diff(F[j],var(P+i))*q;
      p = p + q;
    }
    I = I, p + var(2*P+N+i);
  }
  // -------- the ideal I is ready ----------
  ideal IF = I;
  export IF;
  return(@R);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring R = 0,(x,y,z),Dp;
  ideal I = x^2+y^3, z;
  def W = makeIF(I);
  setring W;
  IF;
}

proc bfctVarIn (ideal I, list #)
"USAGE:  bfctVarIn(I [,a,b,c]);  I an ideal, a,b,c optional ints
RETURN:  list of ideal and intvec
PURPOSE: computes the roots of the Bernstein-Sato polynomial and their
@*       multiplicities for an affine algebraic variety defined by I.
ASSUME:  The basering is commutative and of characteristic 0.
@*       Varnames of the basering do not include t(1),...,t(r) and
@*       Dt(1),...,Dt(r), where r is the number of entries of the input ideal.
BACKGROUND:  In this proc, the initial ideal of the multivariate Malgrange ideal
@*       defined by I is computed and then a system of linear equations is solved
@*       by linear reductions following the ideas by Noro.
NOTE:    In the output list, say L,
@*       - L[1] of type ideal contains all the rational roots of a b-function,
@*       - L[2] of type intvec contains the multiplicities of above roots,
@*       - optional L[3] of type string is the part of b-function without
@*        rational roots.
@*       Note, that a b-function of degree 0 is encoded via L[1][1]=0, L[2]=0 and
@*       L[3] is 1 (for nonzero constant) or 0 (for zero b-function).
@*       If a<>0, the ideal is used as given. Otherwise, and by default, a
@*       heuristically better suited generating set is used to reduce computation
@*       time.
@*       If b<>0, @code{std} is used for GB computations in characteristic 0,
@*       otherwise, and by default, @code{slimgb} is used.
@*       If c<>0, a matrix ordering is used for GB computations, otherwise,
@*       and by default, a block ordering is used.
DISPLAY: If printlevel=1, progress debug messages will be printed,
@*       if printlevel>=2, all the debug messages will be printed.
EXAMPLE: example bfctVarIn; shows examples
"
{
  if (dmodvarAssumeViolation())
  {
    ERROR("Basering is inappropriate: characteristic>0 or qring present");
  }
  if (!isCommutative())
  {
    ERROR("Basering must be commutative");
  }
  int ppl = printlevel - voice + 2;
  int idealasgiven  = 0; // default
  int whicheng      = 0; // default
  int whichord      = 0; // default
  if (size(#)>0)
  {
    if (typeof(#[1])=="int" || typeof(#[1])=="number")
    {
      idealasgiven = int(#[1]);
    }
    if (size(#)>1)
    {
      if (typeof(#[2])=="int" || typeof(#[2])=="number")
      {
        whicheng = int(#[2]);
      }
      if (size(#)>2)
      {
        if (typeof(#[3])=="int" || typeof(#[3])=="number")
        {
          whichord = int(#[3]);
        }
      }
    }
  }
  def save = basering;
  int i;
  int n = nvars(basering);
  // step 0: get small generating set
  I = simplify(I,2);
  list L = smallGenCoDim(I,idealasgiven);
  I = L[1];
  int c = L[2];
  kill L;
  // step 1: setting up the multivariate Malgrange ideal
  int r = size(I);
  def D = makeIF(I);
  setring D;
  dbprint(ppl-1,"// Computing in " + string(n+r) + "-th Weyl algebra:", D);
  dbprint(ppl-1,"// The Malgrange ideal: ", IF);
  // step 2: compute the b-function of the Malgrange ideal w.r.t. approriate weights
  intvec w = 1:r;
  w[r+n] = 0;
  dbprint(ppl,"// Computing the b-function of the Malgrange ideal...");
  list L = bfctIdeal(IF,w,whicheng,whichord);
  dbprint(ppl,"// ... done.");
  dbprint(ppl-1,"// The b-function: ",L);
  // step 3: shift the result
  ring S = 0,s,dp;
  list L = imap(D,L);
  kill D;
  if (size(L)==2)
  {
    ideal B = L[1];
    for (i=1; i<=ncols(B); i++)
    {
      B[i] = -B[i]+c-r-1;
    }
    L[1] = B;
  }
  else // should never get here: BS poly has non-rational roots
  {
    string str = L[3];
    L = delete(L,3);
    str = "poly @b = (" + str + ")*(" + string(fl2poly(L,"s")) + ");";
    execute(str);
    poly b = subst(@b,s,-s+c-r-1);
    L = bFactor(b);
  }
  setring save;
  list L = imap(S,L);
  return(L);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring R = 0,(x,y,z),dp;
  ideal F = x^2+y^3, z;
  list L = bfctVarIn(F);
  L;
}

static proc smallGenCoDim (ideal I, int Iasgiven)
{
  // call from K[x]
  // returns list L
  // L[1]=I or L[1]=smaller generating set of I
  // L[2]=codimension(I)
  int ppl = printlevel - voice + 3;
  int n = nvars(basering);
  int c;
  if (attrib(I,"isSB") == 1)
  {
    c = n - dim(I);
    if (!Iasgiven)
    {
      list L = mstd(I);
    }
  }
  else
  {
    def save = basering;
    list RL = ringlist(save);
    list @ord;
    @ord[1] = list("dp", intvec(1:n));
    @ord[2] = list("C", intvec(0));
    RL[3] = @ord;
    kill @ord;
    if (size(RL)>4) // commutative G-algebra present
    {
      RL = RL[1..4];
    }
    def R = ring(RL);
    kill RL;
    setring R;
    ideal I = imap(save,I);
    if (!Iasgiven)
    {
      list L = mstd(I);
      c = n - dim(L[1]);
      setring save;
      list L = imap(R,L);
    }
    else
    {
      I = std(I);
      c = n - dim(I);
      setring save;
    }
    kill R;
  }
  if (!Iasgiven)
  {
    if (size(L[2]) < size(I))
    {
      I = L[2];
      dbprint(ppl,"// Found smaller generating set of the given variety: ", I);
    }
    else
    {
      dbprint(ppl,"// Have not found smaller generating set of the given variety.");
    }
  }
  dbprint(ppl-1,"// The codim of the given variety is " + string(c) + ".");
  if (!defined(L))
  {
    list L;
  }
  L[1] = I;
  L[2] = c;
  return(L);
}


// Some more examples

static proc TXcups()
{
  "EXAMPLE:"; echo = 2;
  //TX tangent space of X=V(x^2+y^3)
  ring R = 0,(x0,x1,y0,y1),Dp;
  ideal F = x0^2+y0^3, 2*x0*x1+3*y0^2*y1;
  printlevel = 0;
  //ORD = "(a(1,1),a(1,1,1,1,1,1),dp)";
  //eng = 0;
  def A = SannfsVar(F);
  setring A;
  LD;
}

static proc ex47()
{
  ring r7 = 0,(x0,x1,y0,y1),dp;
  ideal I = x0^2+y0^3, 2*x0*x1+3*y0^2*y1;
  bfctVarIn(I);
  // second ex - too big
  ideal J = x0^4+y0^5, 4*x0^3*x1+5*y0^4*y1;
  bfctVarIn(J);
}

static proc ex48()
{
  ring r8 = 0,(x1,x2,x3),dp;
  ideal I = x1^3-x2*x3, x2^2-x1*x3, x3^2-x1^2*x2;
  bfctVarIn(I);
}

static proc ex49 ()
{
  ring r9 = 0,(z1,z2,z3,z4),dp;
  ideal I = z3^2-z2*z4, z2^2*z3-z1*z4, z2^3-z1*z3;
  bfctVarIn(I);
}

static proc ex410()
{
  LIB "toric.lib";
  ring r = 0,(z(1..7)),dp;
  intmat A[3][7];
  A = 6,4,2,0,3,1,0,0,1,2,3,0,1,0,0,0,0,0,1,1,2;
  ideal I = toric_ideal(A,"pt");
  I = std(I);
//ideal I = z(6)^2-z(3)*z(7), z(5)*z(6)-z(2)*z(7), z(5)^2-z(1)*z(7),
//  z(4)*z(5)-z(3)*z(6), z(3)*z(5)-z(2)*z(6), z(2)*z(5)-z(1)*z(6),
//  z(3)^2-z(2)*z(4), z(2)*z(3)-z(1)*z(4), z(2)^2-z(1)*z(3);
  bfctVarIn(I,1); // no result yet
}