///////////////////////////////////////////////////////////////////////////////
version="$Id: algebra.lib,v 1.17 2008-10-06 17:04:26 Singular Exp $";
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 poly 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)
";

 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 poly 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 poly 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;
        ");
     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));
      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 poly 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;
      ");
      return(l);
    }
  }
  else
  { "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 iff 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 iff 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 with a different
         algorithm, which is often faster
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)),dp;");
 //  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));
  //---------- 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)
     { l[2]=l[2]+string(ker[i]);
     }
  }
  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 choosen.
         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;
        ");
    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, 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;
   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+");");
     }
    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, 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 choosen.
         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 iff all the variables of the basering are
         contained in the polynomial subalgebra generated by the polynomials
         defining phi. Hence, 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 means \"surjectivity\"
         in this sense.
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;
       }
    }
   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, 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" );
      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" );
     }
     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:
           if we define  map phi=basering,I;
           then k[var(1),...,var(n)]/phi(id) is finite over k[I].
         - J is generated by a subset of the variables with size(J) = dim(id)
         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
"
{
   if ( deg(i[1]) <= 0)
   {
     list l = ideal(0),i;
     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("dp");
   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 = ideal(0),maxideal(1);
     return( l );
   }
   // ------------------------ change of ordering ---------------------------
   //Now change the order to (dp(n-d),lp) creating a new basering @S
    string s ="dp("+string(n-d)+"),lp";
   def @S=changeord(s);
   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].
EXAMPLE: example finitenessTest; shows an example
"
{  int n = nvars(basering);
   intvec v,w;
   int j,z,ii;
   intvec V = 1..n;
   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]);
      if(size(ideal(w))==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, R, list #)
"USAGE:   mapIsFinite(phi,R[,J]); R a ring, phi: R ---> basering a map
         J an ideal in the basering, J = 0 if not given
RETURN:  1 if R ---> basering/J is finite and 0 else
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);
  return(l[1]);
}
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);
}
//////////////////////////////////////////////////////////////////////////////