source: git/Singular/LIB/algebra.lib @ 3963c1

///////////////////////////////////////////////////////////////////////////////
version="version algebra.lib 4.0.0.0 Jun_2013 "; // $Id$
category="Commutative Algebra";
info="
LIBRARY:  algebra.lib   Compute with Algbras and Algebra Maps
AUTHORS:  Gert-Martin Greuel,     greuel@mathematik.uni-kl.de,
@*        Agnes Eileen Heydtmann, agnes@math.uni-sb.de,
@*        Gerhard Pfister,        pfister@mathematik.uni-kl.de

PROCEDURES:
 algebra_containment(); query of algebra containment
 module_containment();  query of module containment over a subalgebra
 inSubring(p,I);        test whether polynomial p is in subring generated by I
 algDependent(I);       computes algebraic relations between generators of I
 alg_kernel(phi);       computes the kernel of the ringmap phi
 is_injective(phi);     test for injectivity of ringmap phi
 is_surjective(phi);    test for surjectivity of ringmap phi
 is_bijective(phi);     test for bijectivity of ring map phi
 noetherNormal(id);     noether normalization of ideal id
 mapIsFinite(R,phi,I);  query for finiteness of map phi:R --> basering/I

 finitenessTest(i,z);   find variables which occur as pure power in lead(i)
 nonZeroEntry(id);      list describing non-zero entries of an identifier
";

 LIB "inout.lib";
 LIB "elim.lib";
 LIB "ring.lib";
 LIB "matrix.lib";

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

proc algebra_containment (poly p, ideal A, list #)
"USAGE:   algebra_containment(p,A[,k]); p poly, A ideal, k integer.
@*       A = A[1],...,A[m] generators of subalgebra of the basering
RETURN:
@format
         - k=0 (or if k is not given) an integer:
           1  : if p is contained in the subalgebra K[A[1],...,A[m]]
           0  : if p is not contained in K[A[1],...,A[m]]
         - k=1 : a list, say l, of size 2, l[1] integer, l[2] ring, satisfying
           l[1]=1 if p is in the subalgebra K[A[1],...,A[m]] and then the ring
           l[2]: ring, contains poly check = h(y(1),...,y(m)) if p=h(A[1],...,A[m])
           l[1]=0 if p is not in the subalgebra K[A[1],...,A[m]] and then
           l[2] contains the poly check = h(x,y(1),...,y(m)) if p satisfies
           the nonlinear relation p = h(x,A[1],...,A[m]) where
           x = x(1),...,x(n) denote the variables of the basering
@end format
DISPLAY: if k=0 and printlevel >= voice+1 (default) display the polynomial check
NOTE:    The proc inSubring uses a different algorithm which is sometimes
         faster.
THEORY:  The ideal of algebraic relations of the algebra generators A[1],...,
         A[m] is computed introducing new variables y(i) and the product
         order with x(i) >> y(i).
         p reduces to a polynomial only in the y(i) <=> p is contained in the
         subring generated by the polynomials A[1],...,A[m].
EXAMPLE: example algebra_containment; shows an example
"
{ int DEGB = degBound;
  degBound = 0;
    if (size(#)==0)
    { #[1] = 0;
    }
    def br=basering;
    int n = nvars(br);
    int m = ncols(A);
    int i;
    string mp=string(minpoly);
    //-----------------
    // neu CL 10/05:
    int is_qring;
    if (size(ideal(br))>0) {
      is_qring=1;
      ideal IdQ = ideal(br);
    }
    //-----------------
    // ---------- create new ring with extra variables --------------------
    execute ("ring R=("+charstr(br)+"),(x(1..n),y(1..m)),(dp(n),dp(m));");
    execute ("minpoly=number("+mp+");");
    ideal vars=x(1..n);
    map emb=br,vars;
    ideal A=ideal(emb(A));
    poly check=emb(p);
    for (i=1;i<=m;i=i+1)
    { A[i]=A[i]-y(i);
    }
    //-----------------
    // neu CL 10/05:
    if (is_qring) { A = A,emb(IdQ); }
    //-----------------
    A=std(A);
    check=reduce(check,A);
    /*alternatively we could use reduce(check,A,1) which is a little faster
      but result is bigger since it is not tail-reduced
    */
  //--- checking whether all variables from old ring have disappeared ------
  // if so, then the sum of the first n leading exponents is 0, hence i=1
  // use i also to control the display
    i = (sum(leadexp(check),1..n)==0);
  degBound = DEGB;
    if( #[1] == 0 )
    { dbprint(printlevel-voice+3,"// "+string(check));
      return(i);
    }
    else
    { list l = i,R;
      kill A,vars,emb;
      export check;
      dbprint(printlevel-voice+3,"
// 'algebra_containment' created a ring as 2nd element of the list.
// The ring contains the polynomial check which defines the algebraic relation.
// To access to the ring and see check you must give the ring a name,
// e.g.:
               def S = l[2]; setring S; check;
        ");
     setring br;
     return(l);
    }
}
example
{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2;
   int p = printlevel; printlevel = 1;
   ring R = 0,(x,y,z),dp;
   ideal A=x2+y2,z2,x4+y4,1,x2z-1y2z,xyz,x3y-1xy3;
   poly p1=z;
   poly p2=
   x10z3-x8y2z3+2x6y4z3-2x4y6z3+x2y8z3-y10z3+x6z4+3x4y2z4+3x2y4z4+y6z4;
   algebra_containment(p1,A);
   algebra_containment(p2,A);
   list L = algebra_containment(p2,A,1);
   L[1];
   def S = L[2]; setring S;
   check;
   printlevel = p;
}
///////////////////////////////////////////////////////////////////////////////

proc module_containment(poly p, ideal P, ideal S, list #)
"USAGE:   module_containment(p,P,M[,k]); p poly, P ideal, M ideal, k int
@*       P = P[1],...,P[n] generators of a subalgebra of the basering,
@*       M = M[1],...,M[m] generators of a module over the subalgebra K[P]
ASSUME:  ncols(P) = nvars(basering), the P[i] are algebraically independent
RETURN:
@format
         - k=0 (or if k is not given), an integer:
           1    : if p is contained in the module <M[1],...,M[m]> over K[P]
           0    : if p is not contained in <M[1],...,M[m]>
         - k=1, a list, say l, of size 2, l[1] integer, l[2] ring:
           l[1]=1 : if p is in <M[1],...,M[m]> and then the ring l[2] contains
             the polynomial check = h(y(1),...,y(m),z(1),...,z(n)) if
             p = h(M[1],...,M[m],P[1],...,P[n])
           l[1]=0 : if p is in not in <M[1],...,M[m]>, then l[2] contains the
             poly check = h(x,y(1),...,y(m),z(1),...,z(n)) if p satisfies
             the nonlinear relation p = h(x,M[1],...,M[m],P[1],...,P[n]) where
             x = x(1),...,x(n) denote the variables of the basering
@end format
DISPLAY: the polynomial h(y(1),...,y(m),z(1),...,z(n)) if k=0, resp.
         a comment how to access the relation check if k=1, provided
         printlevel >= voice+1 (default).
THEORY:  The ideal of algebraic relations of all the generators p1,...,pn,
         s1,...,st given by P and S is computed introducing new variables y(j),
         z(i) and the product order: x^a*y^b*z^c > x^d*y^e*z^f if x^a > x^d
         with respect to the lp ordering or else if z^c > z^f with respect to
         the dp ordering or else if y^b > y^e with respect to the lp ordering
         again. p reduces to a polynomial only in the y(j) and z(i), linear in
         the z(i) <=> p is contained in the module.
EXAMPLE: example module_containment; shows an example
"
{ def br=basering;
  int DEGB = degBound;
  degBound=0;
  if (size(#)==0)
  { #[1] = 0;
  }
  int n=nvars(br);
  if ( ncols(P)==n )
  { int m=ncols(S);
    string mp=string(minpoly);
  // ---------- create new ring with extra variables --------------------
    execute
   ("ring R=("+charstr(br)+"),(x(1..n),y(1..m),z(1..n)),(lp(n),dp(m),lp(n));");
    execute ("minpoly=number("+mp+");");
    ideal vars = x(1..n);
    map emb = br,vars;
    ideal P = emb(P);
    ideal S  = emb(S);
    poly check = emb(p);
    ideal I;
    for (int i=1;i<=m;i=i+1)
    { I[i]=S[i]-y(i);
    }
    for (i=1;i<=n;i=i+1)
    { I[m+i]=P[i]-z(i);
    }
    I=std(I);
    check = reduce(check,I);
  //--- checking whether all variables from old ring have disappeared ------
  // if so, then the sum of the first n leading exponents is 0
    i = (sum(leadexp(check),1..n)==0);
    if( #[1] == 0 )
    { dbprint(i*(printlevel-voice+3),"// "+string(check));
      setring br;
      return(i);
    }
    else
    { list l = i,R;
      kill I,vars,emb,P,S;
      export check;
      dbprint(printlevel-voice+3,"
// 'module_containment' created a ring as 2nd element of the list. The
// ring contains the polynomial check which defines the algebraic relation
// for p. To access to the ring and see check you must give the ring
// a name, e.g.:
     def S = l[2]; setring S; check;
      ");
      setring br;
      return(l);
    }
  }
  else
  {
      setring br;
    "ERROR: the first ideal must have nvars(basering) entries";
    return();
  }
}
example
{ "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7:"; echo=2;
   int p = printlevel; printlevel = 1;
   ring R=0,(x,y,z),dp;
   ideal P = x2+y2,z2,x4+y4;           //algebra generators
   ideal M = 1,x2z-1y2z,xyz,x3y-1xy3;  //module generators
   poly p1=
   x10z3-x8y2z3+2x6y4z3-2x4y6z3+x2y8z3-y10z3+x6z4+3x4y2z4+3x2y4z4+y6z4;
   module_containment(p1,P,M);
   poly p2=z;
   list l = module_containment(p2,P,M,1);
   l[1];
   def S = l[2]; setring S; check;
   printlevel=p;
}
///////////////////////////////////////////////////////////////////////////////

proc inSubring(poly p, ideal I)
"USAGE:   inSubring(p,i); p poly, i ideal
RETURN:
@format
         a list l of size 2, l[1] integer, l[2] string
         l[1]=1 if and only if p is in the subring generated by i=i[1],...,i[k],
                and then l[2] = y(0)-h(y(1),...,y(k)) if p = h(i[1],...,i[k])
         l[1]=0 if and only if p is in not the subring generated by i,
                and then l[2] = h(y(0),y(1),...,y(k)) where p satisfies the
                nonlinear relation h(p,i[1],...,i[k])=0.
@end format
NOTE:    the proc algebra_containment tests the same using a different
         algorithm, which is often faster
         if l[1] == 0 then l[2] may contain more than one relation h(y(0),y(1),...,y(k)),
         separated by comma
EXAMPLE: example inSubring; shows an example
"
{int z=ncols(I);
  int i;
  def gnir=basering;
  int n = nvars(gnir);
  string mp=string(minpoly);
  list l;
  // neu CL 10/05:
  int is_qring;
  if (size(ideal(gnir))>0) {
    is_qring=1;
    ideal IdQ = ideal(gnir);
  }
  // ---------- create new ring with extra variables --------------------
  //the intersection of ideal nett=(p-y(0),I[1]-y(1),...)
  //with the ring k[y(0),...,y(n)] is computed, the result is ker
   execute ("ring r1= ("+charstr(basering)+"),(x(1..n),y(0..z)),lp;");
 //  execute ("ring r1= ("+charstr(basering)+"),(y(0..z),x(1..n)),dp;");
  execute ("minpoly=number("+mp+");");
  ideal va = x(1..n);
  map emb = gnir,va;
  ideal nett = emb(I);
  for (i=1;i<=z;i++)
  { nett[i]=nett[i]-y(i);
  }
  nett=emb(p)-y(0),nett;
  // neu CL 10/05:
  if (is_qring) { nett = nett,emb(IdQ); }
  //-----------------
  ideal ker=eliminate(nett,product(va));
  ker=std(ker);
  //---------- test wether y(0)-h(y(1),...,y(z)) is in ker --------------
  l[1]=0;
  l[2]="";
  for (i=1;i<=size(ker);i++)
  { if (deg(ker[i]/y(0))==0)
     { string str=string(ker[i]);
        setring gnir;
        l[1]=1;
        l[2]=str;
        return(l);
     }
     if (deg(ker[i]/y(0))>0)
     { if( l[2] != "" ){ l[2] = l[2] + ","; }
       l[2] = l[2] + string(ker[i]);
     }
  }
  setring gnir;
  return(l);
}
example
{ "EXAMPLE:"; echo = 2;
   ring q=0,(x,y,z,u,v,w),dp;
   poly p=xyzu2w-1yzu2w2+u4w2-1xu2vw+u2vw2+xyz-1yzw+2u2w-1xv+vw+2;
   ideal I =x-w,u2w+1,yz-v;
   inSubring(p,I);
}
//////////////////////////////////////////////////////////////////////////////

proc algDependent( ideal A, list # )
"USAGE:   algDependent(f[,c]); f ideal (say, f = f1,...,fm), c integer
RETURN:
@format
         a list l  of size 2, l[1] integer, l[2] ring:
         - l[1] = 1 if f1,...,fm are algebraic dependent, 0 if not
         - l[2] is a ring with variables x(1),...,x(n),y(1),...,y(m) if the
           basering has n variables. It contains the ideal 'ker', depending
           only on the y(i) and generating the algebraic relations between the
           f[i], i.e. substituting y(i) by fi yields 0. Of course, ker is
           nothing but the kernel of the ring map
              K[y(1),...,y(m)] ---> basering,  y(i) --> fi.
@end format
NOTE:    Three different algorithms are used depending on c = 1,2,3.
         If c is not given or c=0, a heuristically best method is chosen.
         The basering may be a quotient ring.
         To access to the ring l[2] and see ker you must give the ring a name,
         e.g. def S=l[2]; setring S; ker;
DISPLAY: The above comment is displayed if printlevel >= 0 (default).
EXAMPLE: example algDependent; shows an example
"
{
    int bestoption = 3;
    // bestoption is the default algorithm, it may be set to 1,2 or 3;
    // it should be changed, if another algorithm turns out to be faster
    // in general. Is perhaps dependent on the input (# vars, size ...)
    int tt;
    if( size(#) > 0 )
    { if( typeof(#[1]) == "int" )
      { tt = #[1];
      }
    }
    if( size(#) == 0 or tt == 0 )
    { tt = bestoption;
    }
    def br=basering;
    int n = nvars(br);
    ideal B = ideal(br);
    int m = ncols(A);
    int s = size(B);
    int i;
    string mp = string(minpoly);
 // --------------------- 1st variant of algorithm ----------------------
 // use internal preimage command (should be equivalent to 2nd variant)
    if ( tt == 1 )
    {
      execute ("ring R1=("+charstr(br)+"),y(1..m),dp;");
      execute ("minpoly=number("+mp+");");
      setring br;
      map phi = R1,A;
      setring R1;
      ideal ker = preimage(br,phi,B);
    }
 // ---------- create new ring with extra variables --------------------
    execute ("ring R2=("+charstr(br)+"),(x(1..n),y(1..m)),(dp);");
    execute ("minpoly=number("+mp+");");
    if( tt == 1 )
    {
      ideal ker = imap(R1,ker);
    }
    else
    {
      ideal vars = x(1..n);
      map emb = br,vars;
      ideal A = emb(A);
      for (i=1; i<=m; i=i+1)
      { A[i] = A[i]-y(i);
      }
 // --------------------- 2nd variant of algorithm ----------------------
 // use internal eliminate for eliminating m variables x(i) from
 // ideal A[i] - y(i) (uses extra eliminating 'first row', a-order)
      if ( s == 0 and  tt == 2  )
      { ideal ker = eliminate(A,product(vars));
      }
      else
 // eliminate does not work in qrings
 // --------------------- 3rd variant of algorithm ----------------------
 // eliminate m variables x(i) from ideal A[i] - y(i) by choosing product
 // order
       {execute ("ring R3=("+charstr(br)+"),(x(1..n),y(1..m)),(dp(n),dp(m));");
        execute ("minpoly=number("+mp+");");
        if ( s != 0 )
        { ideal vars = x(1..n);
          map emb = br,vars;
          ideal B = emb(B);
          attrib(B,"isSB",1);
          qring Q = B;
        }
        ideal A = imap(R2,A);
        A = std(A);
        ideal ker = nselect(A,1..n);
        setring R2;
        if ( defined(Q)==voice )
        { ideal ker = imap(Q,ker);
        }
        else
        { ideal ker = imap(R3,ker);
        }
        kill A,emb,vars;
      }
    }
 // --------------------------- return ----------------------------------
    s = size(ker);
    list L = (s!=0), R2;
    export(ker);
    dbprint(printlevel-voice+3,"
// The 2nd element of the list l is a ring with variables x(1),...,x(n),
// and y(1),...,y(m) if the basering has n variables and if the ideal
// is f[1],...,f[m]. The ring contains the ideal ker which depends only
// on the y(i) and generates the relations between the f[i].
// I.e. substituting y(i) by f[i] yields 0.
// To access to the ring and see ker you must give the ring a name,
// e.g.:
             def S = l[2]; setring S; ker;
        ");
    setring br;
    return (L);
}
example
{ "EXAMPLE:"; echo = 2;
   int p = printlevel; printlevel = 1;
   ring R = 0,(x,y,z,u,v,w),dp;
   ideal I = xyzu2w-1yzu2w2+u4w2-1xu2vw+u2vw2+xyz-1yzw+2u2w-1xv+vw+2,
             x-w, u2w+1, yz-v;
   list l = algDependent(I);
   l[1];
   def S = l[2]; setring S;
   ker;
   printlevel = p;
}
//////////////////////////////////////////////////////////////////////////////
proc alg_kernel( map phi, def pr, list #)
"USAGE:   alg_kernel(phi,pr[,s,c]); phi map to basering, pr preimage ring,
         s string (name of kernel in pr), c integer.
RETURN:  a string, the kernel of phi as string.
         If, moreover, a string s is given, the algorithm creates, in the
         preimage ring pr the kernel of phi with name s.
         Three different algorithms are used depending on c = 1,2,3.
         If c is not given or c=0, a heuristically best method is chosen.
         (algorithm 1 uses the preimage command)
NOTE:    Since the kernel of phi lives in pr, it cannot be returned to the
         basering. If s is given, the user has access to it in pr via s.
         The basering may be a quotient ring.
EXAMPLE: example alg_kernel; shows an example
"
{   int tt;
   def BAS = basering;
   if( size(#) >0 )
   { if( typeof(#[1]) == "int")
     { tt = #[1];
     }
     if( typeof(#[1]) == "string")
     { string nker=#[1];
     }
     if( size(#)>1 )
     {  if( typeof(#[2]) == "string")
        { string nker=#[2];
        }
        if( typeof(#[2]) == "int")
       {  tt = #[2];
       }
     }
   }
    int n = nvars(basering);
    string mp = string(minpoly);
    ideal A = ideal(phi);
    //def pr = preimage(phi);
    //folgendes Auffuellen oder Stutzen ist ev nicht mehr noetig
    //falls map das richtig macht
    int m = nvars(pr);
    ideal j;
    j[m]=0;
    A=A,j;
    A=A[1..m];
    list L = algDependent(A,tt);
    // algDependent is called with "bestoption" if tt = 0
    def S = L[2];
    execute ("ring R=("+charstr(basering)+"),(@(1..n),"+varstr(pr)+"),(dp);");
    execute ("minpoly=number("+mp+");");
    ideal ker = fetch(S,ker);       //in order to have variable names correct
    string sker = string(ker);
    if (defined(nker) == voice)
    { setring pr;
      execute("ideal "+nker+"="+sker+";");
      execute("export("+nker+");");
     }
    setring BAS;
    return(sker);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(a,b,c),ds;
   ring s = 0,(x,y,z,u,v,w),dp;
   ideal I = x-w,u2w+1,yz-v;
   map phi = r,I;                // a map from r to s:
   alg_kernel(phi,r);            // a,b,c ---> x-w,u2w+1,yz-v

   ring S = 0,(a,b,c),ds;
   ring R = 0,(x,y,z),dp;
   qring Q = std(x-y);
   ideal i = x, y, x2-y3;
   map phi = S,i;                 // a map to a quotient ring
   alg_kernel(phi,S,"ker",3);     // uses algorithm 3
   setring S;                     // you have access to kernel in preimage
   ker;
}
//////////////////////////////////////////////////////////////////////////////

proc is_injective( map phi,def pr,list #)
"USAGE:   is_injective(phi,pr[,c,s]); phi map, pr preimage ring, c int, s string
RETURN:
@format
         - 1 (type int) if phi is injective, 0 if not (if s is not given).
         - If s is given, return a list l of size 2, l[1] int, l[2] ring:
           l[1] is 1 if phi is injective, 0 if not
           l[2] is a ring with variables x(1),...,x(n),y(1),...,y(m) if the
           basering has n variables and the map m components, it contains the
           ideal 'ker', depending only on the y(i), the kernel of the given map
@end format
NOTE:    Three differnt algorithms are used depending on c = 1,2,3.
         If c is not given or c=0, a heuristically best method is chosen.
         The basering may be a quotient ring. However, if the preimage ring is
         a quotient ring, say pr = P/I, consider phi as a map from P and then
         the algorithm returns 1 if the kernel of phi is 0 mod I.
         To access to the ring l[2] and see ker you must give the ring a name,
         e.g. def S=l[2]; setring S; ker;
DISPLAY: The above comment is displayed if printlevel >= 0 (default).
EXAMPLE: example is_injective; shows an example
"
{  def bsr = basering;
   int tt;
   if( size(#) >0 )
   { if( typeof(#[1]) == "int")
     { tt = #[1];
     }
     if( typeof(#[1]) == "string")
     { string pau=#[1];
     }
     if( size(#)>1 )
     {  if( typeof(#[2]) == "string")
        { string pau=#[2];
        }
        if( typeof(#[2]) == "int")
       {  tt = #[2];
       }
     }
   }
    int n = nvars(basering);
    string mp = string(minpoly);
    ideal A = ideal(phi);
    //def pr = preimage(phi);
    //folgendes Auffuellen oder Stutzen ist ev nicht mehr noetig
    //falls map das richtig macht
    int m = nvars(pr);
    ideal j;
    j[m]=0;
    A=A,j;
    A=A[1..m];
    list L = algDependent(A,tt);
    L[1] = L[1]==0;
// the preimage ring may be a quotient ring, we still have to check whether
// the kernel is 0 mod ideal of the quotient ring
    setring pr;
    if ( size(ideal(pr)) != 0 )
    { def S = L[2];
      ideal proj;
      proj [n+1..n+m] = maxideal(1);
      map psi = S,proj;
      L[1] = size(NF(psi(ker),std(0))) == 0;
    }
    if ( defined(pau) != voice )
    {  return (L[1]);
    }
    else
    {
      dbprint(printlevel-voice+3,"
// The 2nd element of the list is a ring with variables x(1),...,x(n),
// y(1),...,y(m) if the basering has n variables and the map is
// F[1],...,F[m].
// It contains the ideal ker, the kernel of the given map y(i) --> F[i].
// To access to the ring and see ker you must give the ring a name,
// e.g.:
     def S = l[2]; setring S; ker;
        ");
      return(L);
    }
 }
example
{ "EXAMPLE:"; echo = 2;
   int p = printlevel;
   ring r = 0,(a,b,c),ds;
   ring s = 0,(x,y,z,u,v,w),dp;
   ideal I = x-w,u2w+1,yz-v;
   map phi = r,I;                    // a map from r to s:
   is_injective(phi,r);              // a,b,c ---> x-w,u2w+1,yz-v
   ring R = 0,(x,y,z),dp;
   ideal i = x, y, x2-y3;
   map phi = R,i;                    // a map from R to itself, z --> x2-y3
   list l = is_injective(phi,R,"");
   l[1];
   def S = l[2]; setring S;
   ker;
}
///////////////////////////////////////////////////////////////////////////////

proc is_surjective( map phi )
"USAGE:   is_surjective(phi); phi map to basering, or ideal defining it
RETURN:  an integer,  1 if phi is surjective, 0 if not
NOTE:    The algorithm returns 1 if and only if all the variables of the basering are
         contained in the polynomial subalgebra generated by the polynomials
         defining phi. Hence, it tests surjectivity in the case of a global odering.
         If the basering has local or mixed ordering or if the preimage ring is a
         quotient ring (in which case the map may not be well defined) then the return
         value 1 needs to be interpreted with care.
EXAMPLE: example is_surjective; shows an example
"
{
  def br=basering;
    ideal B = ideal(br);
    int s = size(B);
    int n = nvars(br);
    ideal A = ideal(phi);
    int m = ncols(A);
    int ii,t=1,1;
    string mp=string(minpoly);
  // ------------ create new ring with extra variables ---------------------
    execute ("ring R=("+charstr(br)+"),(x(1..n),y(1..m)),(dp(n),dp(m));");
    execute ("minpoly=number("+mp+");");
    ideal vars = x(1..n);
    map emb = br,vars;
    if ( s != 0 )
    {  ideal B = emb(B);
       attrib(B,"isSB",1);
       qring Q = B;
       ideal vars = x(1..n);
       map emb = br,vars;
    }
    ideal A = emb(A);
    for ( ii=1; ii<=m; ii++ )
    { A[ii] = A [ii]-y(ii);
    }
    A=std(A);
  // ------------- check whether the x(i) are in the image -----------------
    poly check;
    for (ii=1; ii<=n; ii++ )
    {  check=reduce(x(ii),A,1);
  // --- checking whether all variables from old ring have disappeared -----
  // if so, then the sum of the first n leading exponents is 0
       if( sum(leadexp(check),1..n)!=0 )
       { t=0;
         break;
       }
    }
   setring br;
   return(t);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R = 0,(x,y,z),dp;
   ideal i = x, y, x2-y3;
   map phi = R,i;                    // a map from R to itself, z->x2-y3
   is_surjective(phi);
   qring Q = std(ideal(z-x37));
   map psi = R, x,y,x2-y3;           // the same map to the quotient ring
   is_surjective(psi);

   ring S = 0,(a,b,c),dp;
   map psi = R,ideal(a,a+b,c-a2+b3); // a map from R to S,
   is_surjective(psi);               // x->a, y->a+b, z->c-a2+b3
}

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

proc is_bijective ( map phi,def pr )
"USAGE:   is_bijective(phi,pr); phi map to basering, pr preimage ring
RETURN:  an integer,  1 if phi is bijective, 0 if not
NOTE:    The algorithm checks first injectivity and then surjectivity.
         To interprete this for local/mixed orderings, or for quotient rings
         type help is_surjective; and help is_injective;
DISPLAY: A comment if printlevel >= voice-1 (default)
EXAMPLE: example is_bijective; shows an example
"
{
  def br = basering;
    int n = nvars(br);
    ideal B = ideal(br);
    int s = size(B);
    ideal A = ideal(phi);
    //folgendes Auffuellen oder Stutzen ist ev nicht mehr noetig
    //falls map das richtig macht
    int m = nvars(pr);
    ideal j;
    j[m]=0;
    A=A,j;
    A=A[1..m];
    int ii,t = 1,1;
    string mp=string(minpoly);
  // ------------ create new ring with extra variables ---------------------
    execute ("ring R=("+charstr(br)+"),(x(1..n),y(1..m)),(dp(n),dp(m));");
    execute ("minpoly=number("+mp+");");
    ideal vars = x(1..n);
    map emb = br,vars;
    if ( s != 0 )
    {  ideal B = emb(B);
       attrib(B,"isSB",1);
       qring Q = B;
       ideal vars = x(1..n);
       map emb = br,vars;
    }
    ideal A = emb(A);
    for ( ii=1; ii<=m; ii++ )
    { A[ii] = A[ii]-y(ii);
    }
    A=std(A);
    def bsr = basering;
 // ------- checking whether phi is injective by computing the kernel -------
    ideal ker = nselect(A,1..n);
    t = size(ker);
    setring pr;
    if ( size(ideal(pr)) != 0 )
    {
      ideal proj;
      proj[n+1..n+m] = maxideal(1);
      map psi = bsr,proj;
      t = size(NF(psi(ker),std(0)));
    }
    if ( t != 0 )
    {  dbprint(printlevel-voice+3,"// map not injective" );
      setring br;
      return(0);
    }
   else
 // -------------- checking whether phi is surjective ----------------------
   { t = 1;
     setring bsr;
     poly check;
     for (ii=1; ii<=n; ii++ )
     {  check=reduce(x(ii),A,1);
  // --- checking whether all variables from old ring have disappeared -----
  // if so, then the sum of the first n leading exponents is 0
        if( sum(leadexp(check),1..n)!=0 )
        { t=0;
          break;
        }
     }
     if ( t == 0 )
     {  dbprint(printlevel-voice+3,"// map injective, but not surjective" );
     }
     setring br;
     return(t);
   }
}
example
{ "EXAMPLE:"; echo = 2;
   int p = printlevel;  printlevel = 1;
   ring R = 0,(x,y,z),dp;
   ideal i = x, y, x2-y3;
   map phi = R,i;                      // a map from R to itself, z->x2-y3
   is_bijective(phi,R);
   qring Q = std(z-x2+y3);
   is_bijective(ideal(x,y,x2-y3),Q);

   ring S = 0,(a,b,c,d),dp;
   map psi = R,ideal(a,a+b,c-a2+b3,0); // a map from R to S,
   is_bijective(psi,R);                // x->a, y->a+b, z->c-a2+b3
   qring T = std(d,c-a2+b3);
   map chi = Q,a,b,a2-b3;              // amap between two quotient rings
   is_bijective(chi,Q);

   printlevel = p;
}
///////////////////////////////////////////////////////////////////////////////

proc noetherNormal(ideal i, list #)
"USAGE:   noetherNormal(id[,p]);  id ideal, p integer
RETURN:
@format
         a list l of two ideals, say I,J:
         - I defines a map (coordinate change in the basering), such that:
         - J is generated by a subset of the variables with size(J) = dim(id)
           if we define  map phi=basering,I;
           then k[var(1),...,var(n)]/phi(id) is finite over k[J].
         If p is given, 0<=p<=100, a sparse coordinate change with p percent
         of the matrix entries being 0 (default: p=0 i.e. dense)
@end format
NOTE:    Designed for characteristic 0.It works also in char k > 0 if it
         terminates,but may result in an infinite loop in small characteristic.
EXAMPLE: example noetherNormal; shows an example
"
{
   i=simplify(i,2);
   if ( i== 0)
   {
     list l = maxideal(1),maxideal(1);
     return( l );
   }
   int p;
   if( size(#) != 0 )
   {
     p = #[1];
   }
   def r = basering;
   int n = nvars(r);
   list good;
   // ------------------------ change of ordering ---------------------------
   //a procedure from ring.lib changing the order to dp creating a new
   //basering @R in order to compute the dimension d of i
   def @R=changeord(list(list("dp",1:nvars(basering))));
   setring @R;
   ideal i = imap(r,i);
   list j = mstd(i);
   i = j[2];
   int d = dim(j[1]);
   if ( d <= 0)
   {
     setring r;
     list l = maxideal(1),ideal(0);
     return( l );
   }
   // ------------------------ change of ordering ---------------------------
   //Now change the order to (dp(n-d),lp) creating a new basering @S
   def @S=changeord(list(list("dp",1:(n-d)),list("lp",1:d)));
   setring @S;
   ideal m;

   // ----------------- sparse-random coordinate change  --------------------
   //creating lower triangular random generators for the maximal ideal
   //a procedure from random.lib, as sparse as possible
   if(  char(@S) >  0 )
   {
      m=ideal(sparsetriag(n,n,p,char(@S)+1)*transpose(maxideal(1)));
   }
   if(  char(@S) == 0 )
   {
      if ( voice <= 6 )
      {
        m=ideal(sparsetriag(n,n,p,10)*transpose(maxideal(1)));
      }
     if( voice > 6 and voice <= 11)
     {
        m=ideal(sparsetriag(n,n,p,100)*transpose(maxideal(1)));
      }
      if ( voice > 11 )
      {
        m=ideal(sparsetriag(n,n,p,30000)*transpose(maxideal(1)));
      }
   }

   map phi=r,m;
   //map phi=@R,m;
   ideal i=std(phi(i));

   // ----------------------- test for finiteness ---------------------------
   //We need a test whether the coordinate change was random enough, if yes
   //we are ready, else call noetherNormal again
   list l=finitenessTest(i);

   setring r;
   list l=imap(@S,l);

   if(size(l[3]) == d)                    //the generic case
   {
      good = fetch(@S,m),l[3];
      kill @S,@R;
      return(good);
   }
   else                                   //the bad case
   { kill @S,@R;
      if ( voice >= 21 )
      {
       "// WARNING: In case of a finite ground field";
       "// the characteristic may be too small.";
       "// This could result in an infinte loop.";
       "// Loop in noetherNormal, voice:";, voice;"";
      }
     if ( voice >= 16 )
     {
       "// Switch to dense coordinate change";"";
       return(noetherNormal(i));
     }
     return(noetherNormal(i,p));
   }
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),dp;
   ideal i= xy,xz;
   noetherNormal(i);
}
///////////////////////////////////////////////////////////////////////////////

proc finitenessTest(ideal i, list #)
"USAGE:   finitenessTest(J[,v]); J ideal, v intvec (say v1,...,vr with vi>0)
RETURN:
@format
         a list l with l[1] integer, l[2], l[3], l[4] ideals
         - l[1] = 1 if var(v1),...,var(vr) are in l[2] and 0 else
         - l[2] (resp. l[3]) contains those variables which occur,
           (resp. do not occur) as pure power in the leading term of one of the
           generators of J,
         - l[4] contains those J[i] for which the leading term is a pure power
           of a variable (which is then in l[2])
         (default: v = [1,2,..,nvars(basering)])
@end format
THEORY:  If J is a standard basis of an ideal generated by x_1 - f_1(y),...,
         x_n - f_n with y_j ordered lexicographically and y_j >> x_i, then,
         if y_i appears as pure power in the leading term of J[k], J[k] defines
         an integral relation for y_i over the y_(i+1),... and the f's.
         Moreover, in this situation, if l[2] = y_1,...,y_r, then K[y_1,...y_r]
         is finite over K[f_1..f_n]. If J contains furthermore polynomials
         h_j(y), then K[y_1,...y_z]/<h_j> is finite over K[f_1..f_n].
         For a proof cf. Prop. 3.1.5, p. 214. in [G.-M. Greuel, G. Pfister:
         A SINGULAR Introduction to Commutative Algebra, 2nd Edition,
         Springer Verlag (2007)]
EXAMPLE: example finitenessTest; shows an example
"
{  int n = nvars(basering);
   intvec v,w;
   int j,z,ii;
   v[n]=0;                             //v should have size n
   intvec V = 1..n;
   list nze;                           //the non-zero entries of a leadexp
   if (size(#) != 0 )
   {
     V = #[1];
   }
   intmat W[1][n];                     //create intmat with 1 row, having 1 at
                                       //position V[j], i = 1..size(V), 0 else
   for( j=1; j<=size(V); j++ )
   {
     W[1,V[j]] = 1;
   }
   ideal relation,zero,nonzero;
   // ---------------------- check leading exponents -------------------------

   for(j=1;j<=ncols(i);j++)
   {
      w = leadexp(i[j]);
      nze = nonZeroEntry(w);
      if( nze[1] == 1 )               //the leading term of i[j] is a
      {                               //pure power of some variable
        if( W*w != 0 )                //case: variable has index appearing in V
        {
          relation[size(relation)+1] = i[j];
          v=v+w;
        }
      }
   }
   // ----------------- pick the corresponding variables ---------------------
   //the nonzero entries of v correspond to variables which occur as
   //pure power in the leading term of some polynomial in i

   for(j=1; j<=size(v); j++)
   {
      if(v[j]==0)
      {
         zero[size(zero)+1]=var(j);
      }
      else
      {
        nonzero[size(nonzero)+1]=var(j);
      }
   }
   // ---------------- do we have all pure powers we want? -------------------
   // test this by dividing the product of corresponding variables
   ideal max = maxideal(1);
   max = max[V];
   z = (product(nonzero)/product(max) != 0);
   return(list(z,nonzero,zero,relation));
}
example
{ "EXAMPLE:"; echo = 2;
   ring s = 0,(x,y,z,a,b,c),(lp(3),dp);
   ideal i= a -(xy)^3+x2-z, b -y2-1, c -z3;
   ideal j = a -(xy)^3+x2-z, b -y2-1, c -z3, xy;
   finitenessTest(std(i),1..3);
   finitenessTest(std(j),1..3);
}
///////////////////////////////////////////////////////////////////////////////

proc mapIsFinite(map phi,def R, list #)
"USAGE:   mapIsFinite(phi,R[,J]); R the preimage ring of the map
         phi: R ---> basering
         J an ideal in the basering, J = 0 if not given
RETURN:  1 if R ---> basering/J is finite and 0 else
NOTE:    R may be a quotient ring (this will be ignored since a map R/I-->S/J
         is finite if and only if the induced map R-->S/J is finite).
SEE ALSO: finitenessTest
EXAMPLE: example mapIsFinite; shows an example
"
{
  def bsr = basering;
  ideal J;
  if( size(#) != 0 )
  {
    J = #[1];
  }
  string os = ordstr(bsr);
  int m = nvars(bsr);
  int n = nvars(R);
  ideal PHI = ideal(phi);
  if ( ncols(PHI) < n )
  { PHI[n]=0;
  }
  // --------------------- change of variable names -------------------------
  execute("ring @bsr = ("+charstr(bsr)+"),y(1..m),("+os+");");
  ideal J = fetch(bsr,J);
  ideal PHI = fetch(bsr,PHI);

  // --------------------------- enlarging ring -----------------------------
  execute("ring @rr = ("+charstr(bsr)+"),(y(1..m),x(1..n)),(lp(m),dp);");
  ideal J = imap(@bsr,J);
  ideal PHI = imap(@bsr,PHI);
  ideal M;
  int i;

  for(i=1;i<=n;i++)
  {
    M[i]=x(i)-PHI[i];
  }
  M = std(M+J);
  // ----------------------- test for finiteness ---------------------------
  list l = finitenessTest(M,1..m);
  int result = l[1];
  setring bsr;
  return( result );
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(a,b,c),dp;
   ring s = 0,(x,y,z),dp;
   ideal i= xy;
   map phi= r,(xy)^3+x2+z,y2-1,z3;
   mapIsFinite(phi,r,i);
}
//////////////////////////////////////////////////////////////////////////////

proc nonZeroEntry(def id)
"USAGE:  nonZeroEntry(id); id=object for which the test 'id[i]!=0', i=1,..,N,
         N=size(id) (resp. ncols(id) for id of type ideal or module)
         is defined (e.g. ideal, vector, list of polynomials, intvec,...)
RETURN:  @format
         a list, say l, with l[1] an integer, l[2], l[3] integer vectors:
         - l[1] number of non-zero entries of id
         - l[2] intvec of size l[1] with l[2][i]=i if id[i] != 0
           in case l[1]!=0 (and l[2]=0 if l[1]=0)
         - l[3] intvec with l[3][i]=1 if id[i]!=0 and l[3][i]=0 else
@end format
NOTE:
EXAMPLE: example nonZeroEntry; shows an example
"
{
   int ii,jj,N,n;
   intvec v,V;

   if ( typeof(id) == "ideal" || typeof(id) == "module" )
   {
      N = ncols(id);
   }
   else
   {
     N = size(id);
   }
   for ( ii=1; ii<=N; ii++ )
   {
      V[ii] = 0;
      if ( id[ii] != 0 )
      {
         n++;
         v=v,ii;      //the first entry of v is always 0
         V[ii] = 1;
      }
   }
   if ( size(v) > 1 ) //if id[ii] != 0 for at least one ii delete the first 0
   {
      v = v[2..size(v)];
   }

   list l = n,v,V;
   return(l);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(a,b,c),dp;
   poly f = a3c+b3+c2+a;
   intvec v = leadexp(f);
   nonZeroEntry(v);

   intvec w;
   list L = 37,0,f,v,w;
   nonZeroEntry(L);
}
//////////////////////////////////////////////////////////////////////////////