source: git/Singular/LIB/dmodvar.lib @ 513673

///////////////////////////////////////////////////////////////////////
version="version dmodvar.lib 4.0.0.0 Jun_2013 "; // $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

Support: DFG Graduiertenkolleg 1632 'Experimentelle und konstruktive Algebra'

OVERVIEW: Let K be a field of characteristic 0. Given a polynomial ring R = K[x_1,...,x_n]
and polynomials 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 the structure of a D<S>-module, where D<S> = D(R)
tensored with S over K, where
@* - D(R) is the 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_i = s_{ii}.
@* 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).


PROCEDURES:
bfctVarIn(F[,L]);       computes the roots of the Bernstein-Sato polynomial b(s) of the variety V(F) using initial ideal approach
bfctVarAnn(F[,L]);      computes the roots of the Bernstein-Sato polynomial b(s) of the variety V(F) using Sannfs approach
SannfsVar(F[,O,e]);     computes the annihilator of F^s in the ring D<S>
makeMalgrange(F[,ORD]); creates 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)
*/

/*
  CHANGELOG
  11.10.10 by DA:
  - reformated help strings
  - simplified code
  - add and use of safeVarName
  - renamed makeIF to makeMalgrange
*/


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


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

// testing for consistency of the library:
proc testdmodvarlib ()
{
  example makeMalgrange;
  example bfctVarIn;
  example bfctVarAnn;
  example SannfsVar;
}
//   example coDim;

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

static proc dmodvarAssumeViolation()
{
  // char K = 0, no qring
  if (  (size(ideal(basering)) >0) || (char(basering) >0) )
  {
    ERROR("Basering is inappropriate: characteristic>0 or qring present");
  }
  return();
}

static proc safeVarName (string s, string cv)
// assumes 's' to be a valid variable name
// returns valid var name string @@..@s
{
  string S;
  if (cv == "v")  { S = "," + "," + varstr(basering)  + ","; }
  if (cv == "c")  { S = "," + "," + charstr(basering) + ","; }
  if (cv == "cv") { S = "," + charstr(basering) + "," + varstr(basering) + ","; }
  s = "," + s + ",";
  while (find(S,s) <> 0)
  {
    s[1] = "@";
    s = "," + s;
  }
  s = s[2..size(s)-1];
  return(s)
    }

// 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 coDim; 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 (Weyl algebra tensored with U(gl_P)), containing an ideal LD
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
ASSUME:  The basering is commutative and over a field of characteristic 0.
NOTE:    Activate the output ring D<S> 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.
DISPLAY: If printlevel=1, progress debug messages will be printed,
@* if printlevel>=2, all the debug messages will be printed.
EXAMPLE: example SannfsVar; shows examples
"
{
  dmodvarAssumeViolation();
  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
    {
      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] = safeVarName("Dt("+string(i)+")","cv");
    for (j=1; j<=P; j++)
    {
      RName[P+(i-1)*P+j] = safeVarName("s("+string(i)+")("+string(j)+")","cv");
    }
  }
  // now, create the names for new vars
  list DName;
  for(i=1; i<=N; i++)
  {
    DName[i] = safeVarName("D"+Name[i],"cv");  //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);
        @D[k,P+(i-1)*P+j] = (j==k)*var(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);
            @D[P+(i-1)*P+j,P+(k-1)*P+l] = -(j==k)*var(i*P+l) + (i==l)*var(k*P+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);
      I[(i-1)*P+j] = var(i)*F[j] + var(i*P+j);
    }
  }
  poly p,q;
  for (i=1; i<=N; i++)
  {
    p=0;
    for (j=1; j<=P; j++)
    {
      //       q = Dt(j);
      q = var(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 Dt(i) 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;
  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 commutative and over a field in char 0.
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.
@* Computational remark: The time of computation can be very different depending
   on the chosen generators of F, although the result is always the same.
@* Further note that 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).
   The resulted is shifted by 1-codim(Var(F)) following (BMS06).
DISPLAY: If printlevel=1, progress debug messages will be printed,
@* if printlevel=2, all the debug messages will be printed.
EXAMPLE: example bfctVarAnn; shows examples
"
{
  dmodvarAssumeViolation();
  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;
  int sF = size(F); // no 0 in F
  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<=sF; i++)
  {
    //     S = S + s(i)(i);
    S = S + var((i-1)*sF+i);
  }
  dbprint(ppl,"// -4-1- computing the minimal polynomial of S");
  dbprint(ppl-1,"S = "+string(S));
  vector M = pIntersect(S,K);
  dbprint(ppl,"// -4-2- the minimal polynomial has been computed");
  ring @R3 = 0,s,dp;
  vector M = imap(@R2,M);
  poly p = vec2poly(M);
  dbprint(ppl-1,p);
  dbprint(ppl,"// -5-1- codimension of the variety");
  dbprint(ppl-1,cd);
  dbprint(ppl,"// -5-2- shifting BS(s)=minpoly(s-codim+1)");
  p = subst(p,var(1),var(1)-cd+1);
  dbprint(ppl-1,p);
  dbprint(ppl,"// -5-3- factorization of the minimal polynomial");
  list BS = bFactor(p);
  setring save;
  list BS = imap(@R3,BS);
  kill @R2,@R3;
  return(BS);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring R = 0,(x,y,z),Dp;
  ideal F = x^2+y^3, z;
  bfctVarAnn(F);
}

proc makeMalgrange (ideal F, list #)
  "USAGE:  makeMalgrange(F [,ORD]);  F an ideal, ORD an optional string
RETURN:  ring (Weyl algebra) containing an ideal IF
PURPOSE: create the ideal by Malgrange associated with F = F[1],...,F[P].
NOTE:    Activate the output 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'.
DISPLAY: If printlevel=1, progress debug messages will be printed,
@* if printlevel>=2, all the debug messages will be printed.
EXAMPLE: example makeMalgrange; 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] = safeVarName("t("+string(i)+")","cv");
    DTName[i] = safeVarName("Dt("+string(i)+")","cv");
  }
  //now, create the names for new vars
  list DName;
  for (i=1; i<=N; i++)
  {
    DName[i] = safeVarName("D"+Name[i],"cv");  //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@;
  def @R = Weyl();
  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];
    I[j] = var(j) - F[j];
  }
  poly p,q;
  for (i=1; i<=N; i++)
  {
    p=0;
    for (j=1; j<=P; j++)
    {
      //       q = Dt(j);
      q = var(P+N+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 = makeMalgrange(I);
  setring W;
  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 over a field 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.
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.
@* Further note, that 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.
   The result is shifted by 1-codim(Var(F)) following (BMS06).
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
"
{
  dmodvarAssumeViolation();
  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 = makeMalgrange(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];
    ideal BB;
    int nB = ncols(B);
    for (i=nB; i>0; i--)
    {
      BB[i] = -B[nB-i+1]+c-r-1;
    }
    L[1] = BB;
  }
  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 + 2;
  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
}
*/