source: git/Singular/LIB/algebra.lib @ 92b2f3

///////////////////////////////////////////////////////////////////////////////
version="$Id: algebra.lib,v 1.1 1999-07-19 14:54:07 obachman Exp $";
info="
LIBRARY:  algebra.lib   PROCEDURES FOR COMPUTING WITH ALGBRAS AND MAPS
AUTHORS:  Gert-Martin Greuel, email: greuel@mathematik.uni-kl.de
          Agnes Eileen Heydtmann, email:agnes@math.uni-sb.de
          Gerhard Pfister, email: 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
 algDependence(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
";

 LIB "inout.lib";
 LIB "elim.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:  - 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] contains the poly check = h(y(1),...,y(m)) if p=h(A[1],...,A[m])
           l[1]=0 if p is in not 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
DISPLAY: if k=0, polynomial h(y(1),...,y(m)) if p=h(A[1],...,A[m]) is contained
         in the subalgebra, provided printlevel >= voice+1 (default)
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);
  // ---------- 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);
    }
    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(i*(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. 
// This 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:  - 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]> and then l[2] contains the
             polynomial 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
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. 
// This 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:  a list, say 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.
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;
  // ---------- 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 ("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;
  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 algDependence( ideal A, list # )
"USAGE:   algDependence(f[,c]); f ideal (say, f = f1,...,fm), c integer
RETURN:  a list, say 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 
         onIy on the y(i) and generating the algebraic relations between the 
         f[i], i.e. substituting y(i) by f[i] yields 0.
         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.
         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 algDependence; 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 ----------------------
    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 ----------------------
 // eliminate does not work in qrings
      if ( s == 0 and  tt == 2  )
      { ideal ker = eliminate(A,product(vars));
      }
      else
 // --------------------- 3rd variant of algorithm ----------------------
      { 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, say l, is a ring with variables x(1),...,x(n),
// y(1),...,y(m) if the basering has n variables and the ideal is f[1],...,f[m]
// It contains the ideal ker only depending on the y(i) and generating 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 = algDependence(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 pr, 
         the kernel of phi with name equal to s.
         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.
         (alogrithm 1 uses the preimage command)
NOTE:    The basering may be a quotient ring.
EXAMPLE: example 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 = algDependence(A,tt);
    // algDependence 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
   // setring preimage(phi);      ## wenn realisiert ##
   ker;
}
//////////////////////////////////////////////////////////////////////////////

proc is_injective( map phi, pr,list #)
"USAGE:   is_injective(phi[,c,s]); phi map, pr reimage ring, c int, s string
RETURN:  - 1 (type int) if phi is injective, 0 if not (if s is not given). 
         - If s is given, return a list, say l, of size 2, l[1] int, l[2] ring:
           l[1] = 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
         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.
NOTE:    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 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 = algDependence(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..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, say, 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; printlevel = 1;
   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;
   
   printlevel = p;
 }
///////////////////////////////////////////////////////////////////////////////

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 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); 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 or mixed orderings, type
         help is_surjective;
DISPLAY: 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..m] = maxideal(1);
      map psi = bsr,proj;
      t = size(NF(psi(ker),std(0))) == 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;
}

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