source: git/Singular/LIB/poly.lib @ 8942a5

///////////////////////////////////////////////////////////////////////////////
version="$Id: poly.lib,v 1.29 2000-12-22 14:22:23 greuel Exp $";
category="General purpose";
info="
LIBRARY:  poly.lib      Procedures for Manipulating Polys, Ideals, Modules
AUTHORS:  O. Bachmann, G.-M: Greuel, A. Fruehbis

PROCEDURES:
 cyclic(int);           ideal of cyclic n-roots
 katsura([i]);          katsura [i] ideal
 freerank(poly/...)     rank of coker(input) if coker is free else -1
 is_homog(poly/...);    int, =1 resp. =0 if input is homogeneous resp. not
 is_zero(poly/...);     int, =1 resp. =0 if coker(input) is 0 resp. not
 lcm(ideal);            lcm of given generators of ideal
 maxcoef(poly/...);     maximal length of coefficient occuring in poly/...
 maxdeg(poly/...);      int/intmat = degree/s of terms of maximal order
 maxdeg1(poly/...);     int = [weighted] maximal degree of input
 mindeg(poly/...);      int/intmat = degree/s of terms of minimal order
 mindeg1(poly/...);     int = [weighted] minimal degree of input
 normalize(poly/...);   normalize poly/... such that leading coefficient is 1
 rad_con(p,I);          check radical containment of poly p in ideal I
 content(f);            content of polynomial/vector f
 numerator(n);          numerator of number n
 denominator(n)         denominator of number n
 mod2id(M,iv);          conversion of a module M to an ideal
 id2mod(i,iv);          conversion inverse to mod2id
 subrInterred(i1,i2,iv);interred w.r.t. a subset of variables
          (parameters in square brackets [] are optional)
";

LIB "general.lib";
///////////////////////////////////////////////////////////////////////////////

proc cyclic (int n)
"USAGE:   cyclic(n);  n integer
RETURN:  ideal of cyclic n-roots from 1-st n variables of basering
EXAMPLE: example cyclic; shows examples
"
{
//----------------------------- procedure body --------------------------------
   ideal m = maxideal(1);
   m = m[1..n],m[1..n];
   int i,j;
   ideal s; poly t;
   for ( j=0; j<=n-2; j=j+1 )
   {
      t=0;
      for( i=1;i<=n;i=i+1 ) { t=t+product(m,i..i+j); }
      s=s+t;
   }
   s=s,product(m,1..n)-1;
   return (s);
}
//-------------------------------- examples -----------------------------------
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(u,v,w,x,y,z),lp;
   cyclic(nvars(basering));
   homog(cyclic(5),z);
}
///////////////////////////////////////////////////////////////////////////////

proc katsura
"USAGE: katsura([n]): n integer
RETURN: katsura(n) : n-th katsura ideal of
         (1) newly created and set ring (32003, x(0..n), dp), if
             nvars(basering) < n
         (2) basering, if nvars(basering) >= n
        katsura()  : katsura ideal of basering
EXAMPLE: example katsura; shows examples
"
{
  int n;
  if ( size(#) == 1 && typeof(#[1]) == "int")
  {
    n = #[1] - 1;
    while (1)
    {
      if (defined(basering))
      {
        if (nvars(basering) >= #[1]) {break;}
      }
      ring katsura_ring = 32003, x(0..#[1]), dp;
      keepring katsura_ring;
      break;
    }
  }
  else
  {
    n = nvars(basering) -1;
  }

  ideal s;
  int i, j;
  poly p;

  p = -1;
  for (i = -n; i <= n; i++)
  {
    p = p + kat_var(i, n);
  }
  s[1] = p;

  for (i = 0; i < n; i++)
  {
    p = -1 * kat_var(i,n);
    for (j = -n; j <= n; j++)
    {
      p = p + kat_var(j,n) * kat_var(i-j, n);
    }
    s = s,p;
  }
  return (s);
}
//-------------------------------- examples -----------------------------------
example
{
  "EXAMPLE:"; echo = 2;
  ring r; basering;
  katsura();
  katsura(4); basering;
}

proc kat_var(int i, int n)
{
  poly p;
  if (i < 0)  { i = -i;}
  if (i <= n) { p = var(i+1); }
  return (p);
}
///////////////////////////////////////////////////////////////////////////////

proc freerank
"USAGE:   freerank(M[,any]);  M=poly/ideal/vector/module/matrix
COMPUTE: rank of module presented by M in case it is free. By definition this
         is vdim(coker(M)/m*coker(M)) if coker(M) is free, where m = maximal
         ideal of basering and M is considered as matrix (the 0-module is
         free of rank 0)
RETURN:  rank of coker(M) if coker(M) is free and -1 else;
         in case of a second argument return a list:
                L[1] = rank of coker(M) or -1
                L[2] = minbase(M)
NOTE:    freerank(syz(M)); computes the rank of M if M is free (and -1 else)
EXAMPLE: example freerank; shows examples
"
{
  int rk;
  def M = simplify(#[1],10);
  resolution mre = res(M,2);
  intmat B = betti(mre);
  if ( ncols(B)>1 ) { rk = -1; }
  else { rk = sum(B[1..nrows(B),1]); }
  if (size(#) == 2) { list L=rk,mre[1]; return(L);}
  return(rk);
}
example
{"EXAMPLE";   echo=2;
  ring r;
  ideal i=x;
  module M=[x,0,1],[-x,0,-1];
  freerank(M);           // should be 2, coker(M) is not free
  freerank(syz (M),"");  // [1] should be 1, coker(syz(M))=M is free of rank 1
                         // [2] should be gen(2)+gen(1) (minimal relation of M)
  freerank(i);
  freerank(syz(i));      // should be 1, coker(syz(i))=i is free of rank 1
}
///////////////////////////////////////////////////////////////////////////////

proc is_homog (id)
"USAGE:   is_homog(id);  id  poly/ideal/vector/module/matrix
RETURN:  integer which is 1 if input is homogeneous (resp. weighted homogeneous
         if the monomial ordering consists of one block of type ws,Ws,wp or Wp,
         assuming that all weights are positive) and 0 otherwise
NOTE:    A vector is homogeneous, if the components are homogeneous of same
         degree, a module/matrix is homogeneous if all column vectors are
         homogeneous
         //*** ergaenzen, wenn Matrizen-Spalten Gewichte haben
EXAMPLE: example is_homog; shows examples
"
{
//----------------------------- procedure body --------------------------------
   module M = module(matrix(id));
   M = simplify(M,2);                        // remove 0-columns
   intvec v = ringweights(basering);
   int i,j=1,1;
   for (i=1; i<=ncols(M); i=i+1)
   {
      if( M[i]!=jet(M[i],deg(lead(M[i])),v)-jet(M[i],deg(lead(M[i]))-1,v))
      { return(0); }
   }
   return(1);
}
//-------------------------------- examples -----------------------------------
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y,z),wp(1,2,3);
   is_homog(x5-yz+y3);
   ideal i = x6+y3+z2, x9-z3;
   is_homog(i);
   ring s = 0,(a,b,c),ds;
   vector v = [a2,0,ac+bc];
   vector w = [a3,b3,c4];
   is_homog(v);
   is_homog(w);
}
///////////////////////////////////////////////////////////////////////////////

proc is_zero
"USAGE:   is_zero(M[,any]); M=poly/ideal/vector/module/matrix
RETURN:  integer, 1 if coker(M)=0 resp. 0 if coker(M)!=0, where M is considered
         as matrix
         if a second argument is given, return a list:
                L[1] = 1 if coker(M)=0 resp. 0 if coker(M)!=0
                L[2] = dim(M)
EXAMPLE: example is_zero; shows examples
"
{
  int d=dim(std(#[1]));
  int a = ( d==-1 );
  if( size(#) >1 ) { list L=a,d; return(L); }
  return(a);
}
example
{ "EXAMPLE:";   echo=2;
  ring r;
  module m = [x],[y],[1,z];
  is_zero(m,1);
  qring q = std(ideal(x2+y3+z2));
  ideal j = x2+y3+z2-37;
  is_zero(j);
}
////////////////////////////////////////////////////////////////////////////////

proc maxcoef (f)
"USAGE:   maxcoef(f);  f  poly/ideal/vector/module/matrix
RETURN:  maximal length of coefficient of f of type int (by counting the
         length of the string of each coefficient)
EXAMPLE: example maxcoef; shows examples
"
{
//----------------------------- procedure body --------------------------------
   int max,s,ii,jj; string t;
   ideal i = ideal(matrix(f));
   i = simplify(i,6);            // delete 0's and keep first of equal elements
   poly m = var(1); matrix C;
   for (ii=2;ii<=nvars(basering);ii=ii+1) { m = m*var(ii); }
   for (ii=1; ii<=size(i); ii=ii+1)
   {
      C = coef(i[ii],m);
      for (jj=1; jj<=ncols(C); jj=jj+1)
      {
         t = string(C[2,jj]);  s = size(t);
         if ( t[1] == "-" ) { s = s - 1; }
         if ( s > max ) { max = s; }
      }
   }
   return(max);
}
//-------------------------------- examples -----------------------------------
example
{ "EXAMPLE:"; echo = 2;
   ring r= 0,(x,y,z),ds;
   poly g = 345x2-1234567890y+7/4z;
   maxcoef(g);
   ideal i = g,10/1234567890;
   maxcoef(i);
   // since i[2]=1/123456789
}
///////////////////////////////////////////////////////////////////////////////

proc maxdeg (id)
"USAGE:   maxdeg(id);  id  poly/ideal/vector/module/matrix
RETURN:  int/intmat, each component equals maximal degree of monomials in the
         corresponding component of id, independent of ring ordering
         (maxdeg of each var is 1)
         of type int if id is of type poly, of type intmat else
NOTE:    proc maxdeg1 returns 1 integer, the absolut maximum; moreover, it has
         an option for computing weighted degrees
EXAMPLE: example maxdeg; shows examples
"
{
   //-------- subprocedure to find maximal degree of given component ----------
   proc findmaxdeg
   {
      poly c = #[1];
      if (c==0) { return(-1); }
   //--- guess upper 'o' and lower 'u' bound, in case of negative weights -----
      int d = (deg(c)>=0)*deg(c)-(deg(c)<0)*deg(c);
      int i = d;
      while ( c-jet(c,i) != 0 ) { i = 2*(i+1); }
      int o = i-1;
      int u = (d != i)*((i /  2)-1);
   //----------------------- "quick search" for maxdeg ------------------------
      while ( (c-jet(c,i)==0)*(c-jet(c,i-1)!=0) == 0)
      {
         i = (o+1+u) /  2;
         if (c-jet(c,i)!=0) { u = i+1; }
         else { o = i-1; }
      }
      return(i);
   }
//------------------------------ main program ---------------------------------
   matrix M = matrix(id);
   int r,c = nrows(M), ncols(M); int i,j;
   intmat m[r][c];
   for (i=r; i>0; i=i-1)
   {
      for (j=c; j>0; j=j-1) { m[i,j] = findmaxdeg(M[i,j]); }
   }
   if (typeof(id)=="poly") { return(m[1,1]); }
   return(m);
}
//-------------------------------- examples -----------------------------------
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y,z),wp(-1,-2,-3);
   poly f = x+y2+z3;
   deg(f);               //deg; returns weighted degree (in case of 1 block)!
   maxdeg(f);
   matrix m[2][2]=f+x10,1,0,f^2;
   maxdeg(m);
}
///////////////////////////////////////////////////////////////////////////////

proc maxdeg1 (id,list #)
"USAGE:   maxdeg1(id[,v]);  id=poly/ideal/vector/module/matrix, v=intvec
RETURN:  integer, maximal [weighted] degree of monomials of id independent of
         ring ordering, maxdeg1 of i-th variable is v[i] (default: v=1..1).
NOTE:    This proc returns one integer while maxdeg returns, in general,
         a matrix of integers. For one polynomial and if no intvec v is given
         maxdeg is faster
EXAMPLE: example maxdeg1; shows examples
"
{
   //-------- subprocedure to find maximal degree of given component ----------
   proc findmaxdeg
   {
      poly c = #[1];
      if (c==0) { return(-1); }
      intvec v = #[2];
   //--- guess upper 'o' and lower 'u' bound, in case of negative weights -----
      int d = (deg(c)>=0)*deg(c)-(deg(c)<0)*deg(c);
      int i = d;
      if ( c == jet(c,-1,v))      //case: maxdeg is negative
      {
         i = -d;
         while ( c == jet(c,i,v) ) { i = 2*(i-1); }
         int o = (d != -i)*((i /  2)+2) - 1;
         int u = i+1;
         int e = -1;
      }
      else                        //case: maxdeg is nonnegative
      {
         while ( c != jet(c,i,v) ) { i = 2*(i+1); }
         int o = i-1;
         int u = (d != i)*((i /  2)-1);
         int e = 1;
      }
   //----------------------- "quick search" for maxdeg ------------------------
      while ( ( c==jet(c,i,v) )*( c!=jet(c,i-1,v) ) == 0 )
      {
         i = (o+e+u) /  2;
         if ( c!=jet(c,i,v) ) { u = i+1; }
         else { o = i-1; }
      }
      return(i);
   }
//------------------------------ main program ---------------------------------
   ideal M = simplify(ideal(matrix(id)),8);   //delete scalar multiples from id
   int c = ncols(M);
   int i,n;
   if( size(#)==0 )
   {
      int m = maxdeg(M[c]);
      for (i=c-1; i>0; i=i-1)
      {
          n = maxdeg(M[i]);
          m = (m>=n)*m + (m<n)*n;             //let m be the maximum of m and n
      }
   }
   else
   {
      intvec v=#[1];                          //weight vector for the variables
      int m = findmaxdeg(M[c],v);
      for (i=c-1; i>0; i--)
      {
         n = findmaxdeg(M[i],v);
         if( n>m ) { m=n; }
      }
   }
   return(m);
}
//-------------------------------- examples -----------------------------------
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y,z),wp(-1,-2,-3);
   poly f = x+y2+z3;
   deg(f);                  //deg returns weighted degree (in case of 1 block)!
   maxdeg1(f);
   intvec v = ringweights(r);
   maxdeg1(f,v);                             //weighted maximal degree
   matrix m[2][2]=f+x10,1,0,f^2;
   maxdeg1(m,v);                             //absolut weighted maximal degree
}
///////////////////////////////////////////////////////////////////////////////

proc mindeg (id)
"USAGE:   mindeg(id);  id  poly/ideal/vector/module/matrix
RETURN:  minimal degree/s of monomials of id, independent of ring ordering
         (mindeg of each variable is 1) of type int if id of type poly, else
         of type intmat.
NOTE:    proc mindeg1 returns one integer, the absolut minimum; moreover it
         has an option for computing weighted degrees.
EXAMPLE: example mindeg; shows examples
"
{
   //--------- subprocedure to find minimal degree of given component ---------
   proc findmindeg
   {
      poly c = #[1];
      if (c==0) { return(-1); }
   //--- guess upper 'o' and lower 'u' bound, in case of negative weights -----
      int d = (ord(c)>=0)*ord(c)-(ord(c)<0)*ord(c);
      int i = d;
      while ( jet(c,i) == 0 ) { i = 2*(i+1); }
      int o = i-1;
      int u = (d != i)*((i /  2)-1);
   //----------------------- "quick search" for mindeg ------------------------
      while ( (jet(c,u)==0)*(jet(c,o)!=0) )
      {
         i = (o+u) / 2;
         if (jet(c,i)==0) { u = i+1; }
         else { o = i-1; }
      }
      if (jet(c,u)!=0) { return(u); }
      else { return(o+1); }
   }
//------------------------------ main program ---------------------------------
   matrix M = matrix(id);
   int r,c = nrows(M), ncols(M); int i,j;
   intmat m[r][c];
   for (i=r; i>0; i=i-1)
   {
      for (j=c; j>0; j=j-1) { m[i,j] = findmindeg(M[i,j]); }
   }
   if (typeof(id)=="poly") { return(m[1,1]); }
   return(m);
}
//-------------------------------- examples -----------------------------------
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y,z),ls;
   poly f = x5+y2+z3;
   ord(f);                      // ord returns weighted order of leading term!
   mindeg(f);                   // computes minimal degree
   matrix m[2][2]=x10,1,0,f^2;
   mindeg(m);                   // computes matrix of minimum degrees
}
///////////////////////////////////////////////////////////////////////////////

proc mindeg1 (id, list #)
"USAGE:   mindeg1(id[,v]);  id=poly/ideal/vector/module/matrix, v=intvec
RETURN:  integer, minimal [weighted] degree of monomials of id independent of
         ring ordering, mindeg1 of i-th variable is v[i] (default v=1..1).
NOTE:    This proc returns one integer while mindeg returns, in general,
         a matrix of integers. For one polynomial and if no intvec v is given
         mindeg is faster.
EXAMPLE: example mindeg1; shows examples
"
{
   //--------- subprocedure to find minimal degree of given component ---------
   proc findmindeg
   {
      poly c = #[1];
      intvec v = #[2];
      if (c==0) { return(-1); }
   //--- guess upper 'o' and lower 'u' bound, in case of negative weights -----
      int d = (ord(c)>=0)*ord(c)-(ord(c)<0)*ord(c);
      int i = d;
      if ( jet(c,-1,v) !=0 )      //case: mindeg is negative
      {
         i = -d;
         while ( jet(c,i,v) != 0 ) { i = 2*(i-1); }
         int o = (d != -i)*((i /  2)+2) - 1;
         int u = i+1;
         int e = -1; i=u;
      }
      else                        //case: inded is nonnegative
      {
         while ( jet(c,i,v) == 0 ) { i = 2*(i+1); }
         int o = i-1;
         int u = (d != i)*((i /  2)-1);
         int e = 1; i=u;
      }
   //----------------------- "quick search" for mindeg ------------------------
      while ( (jet(c,i-1,v)==0)*(jet(c,i,v)!=0) == 0 )
      {
         i = (o+e+u) /  2;
         if (jet(c,i,v)==0) { u = i+1; }
         else { o = i-1; }
      }
      return(i);
   }
//------------------------------ main program ---------------------------------
   ideal M = simplify(ideal(matrix(id)),8);   //delete scalar multiples from id
   int c = ncols(M);
   int i,n;
   if( size(#)==0 )
   {
      int m = mindeg(M[c]);
      for (i=c-1; i>0; i=i-1)
      {
          n = mindeg(M[i]);
          m = (m<=n)*m + (m>n)*n;             //let m be the maximum of m and n
      }
   }
   else
   {
      intvec v=#[1];                          //weight vector for the variables
      int m = findmindeg(M[c],v);
      for (i=c-1; i>0; i=i-1)
      {
         n = findmindeg(M[i],v);
         m = (m<=n)*m + (m>n)*n;              //let m be the maximum of m and n
      }
   }
   return(m);
}
//-------------------------------- examples -----------------------------------
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y,z),ls;
   poly f = x5+y2+z3;
   ord(f);                      // ord returns weighted order of leading term!
   intvec v = 1,-3,2;
   mindeg1(f,v);                // computes minimal weighted degree
   matrix m[2][2]=x10,1,0,f^2;
   mindeg1(m,1..3);             // computes absolut minimum of weighted degrees
}
///////////////////////////////////////////////////////////////////////////////

proc normalize (id)
"USAGE:   normalize(id);  id=poly/vector/ideal/module
RETURN:  object of same type with leading coefficient equal to 1
EXAMPLE: example normalize; shows an example
"
{
   return(simplify(id,1));
}
//-------------------------------- examples -----------------------------------
example
{  "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y,z),ls;
   poly f = 2x5+3y2+4z3;
   normalize(f);
   module m=[9xy,0,3z3],[4z,6y,2x];
   normalize(m);
   ring s = 0,(x,y,z),(c,ls);
   module m=[9xy,0,3z3],[4z,6y,2x];
   normalize(m);
   normalize(matrix(m));             // by automatic type conversion to module!
}
///////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////
// Input: <ideal>=<f1,f2,...,fm> and <polynomial> g
// Question: Does g lie in the radical of <ideal>?
// Solution: Compute a standard basis G for <f1,f2,...,fm,gz-1> where z is a new
//           variable. Then g is contained in the radical of <ideal> <=> 1 is
//           generator in G.
////////////////////////////////////////////////////////////////////////////////
proc rad_con (poly g,ideal I)
"  USAGE:   rad_con(<poly>,<ideal>);
  RETURNS: 1 (TRUE) (type <int>) if <poly> is contained in the radical of
           <ideal>, 0 (FALSE) (type <int>) otherwise
  EXAMPLE: example rad_con; shows an example
"
{ def br=basering;
  int n=nvars(br);
  int dB=degBound;
  degBound=0;
  string mp=string(minpoly);
  execute("ring R=("+charstr(br)+"),(x(1..n),z),dp;");
  execute("minpoly=number("+mp+");");
  ideal irrel=x(1..n);
  map f=br,irrel;
  poly p=f(g);
  ideal J=f(I)+ideal(p*z-1);
  J=std(J);
  degBound=dB;
  if (J[1]==1)
  { return(1);
  }
  else
  { return(0);
  }
}
example
{ "  EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7.";
  echo=2;
           ring R=0,(x,y,z),dp;
           ideal I=x2+y2,z2;
           poly f=x4+y4;
           rad_con(f,I);
           ideal J=x2+y2,z2,x4+y4;
           poly g=z;
           rad_con(g,I);
}

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

proc lcm (ideal i)
"USAGE:   lcm(i); i ideal
RETURN:  poly = lcm(i[1],...,i[size(i)])
NOTE:
EXAMPLE: example lcm; shows an example
"
{
  int k,j;
   poly p,q;
  i=simplify(i,10);
  for(j=1;j<=size(i);j++)
  {
    if(deg(i[j])>0)
    {
      p=i[j];
      break;
    }
  }
  if(deg(p)==-1)
  {
    return(1);
  }
  for (k=j+1;k<=size(i);k++)
  {
     if(deg(i[k])!=0)
     {
        q=gcd(p,i[k]);
        if(deg(q)==0)
        {
           p=p*i[k];
        }
        else
        {
           p=p/q;
           p=p*i[k];
        }
     }
   }
  return(p);
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 0,(x,y,z),lp;
   poly  p = (x+y)*(y+z);
   poly  q = (z4+2)*(y+z);
   ideal l=p,q;
   poly  pr= lcm(l);
   pr;
   l=1,-1,p,1,-1,q,1;
   pr=lcm(l);
   pr;
}

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

proc content(f)
"USAGE:   content(f); f polynomial/vector
RETURN:  number, the content (greatest common factor of coefficients)
         of the polynomial/vector f
EXAMPLE: example content; shows an example
"
{
  return(leadcoef(f)/leadcoef(cleardenom(f)));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),(c,lp);
   content(3x2+18xy-27xyz);
   vector v=[3x2+18xy-27xyz,15x2+12y4,3];
   content(v);
}
///////////////////////////////////////////////////////////////////////////////


proc numerator(number n)
"USAGE:    numerator(n); n number
RETURN:   number, the numerator of n
SEE ALSO: denominator, content, cleardenom
EXAMPLE:  example numerator; shows an example
"
{
  poly p = cleardenom(n+var(1));
  return (coeffs(p,var(1))[1,1]);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,x, dp;
  number n = 3/2;
  numerator(n);
}

proc denominator(number n)
"USAGE:    denominator(n); n number
RETURN:   number, the denominator of n
SEE ALSO: denominator, content, cleardenom
EXAMPLE:  example denominator; shows an example
"
{
  poly p = cleardenom(n+var(1));
  return (coeffs(p,var(1))[2,1]);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,x, dp;
  number n = 3/2;
  denominator(n);
}
////////////////////////////////////////////////////////////////////////

proc mod2id(matrix M,intvec vpos)
"USAGE:     mod2id(M,vpos); M matrix, vpos intvec
ASSUME:    vpos is an integer vector such that gen(i) corresponds
           to var(vpos[i])
           the basering contains variables var(vpos[i]) which do not occur
           in M
NOTE:      this procedure should be used in the following situation:
           one wants to pass to a ring with new variables, say e(1),..,e(s),
           which correspond to the components gen(1),..,gen(s) of the
           module M such that e(i)*e(j)=0 for all i,j
           the new ring should already exist and be the current ring
RETURN:    ideal i in which each gen(i) from the module is replaced by
           var(vpos[i]) and all monomials var(vpos[i])*var(vpos[j]) have
           been added to the generating set of i
EXAMPLE:   example mod2id; shows an example"
{
  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
//----------------------------------------------------------------------
// define the ideal generated by the var(vpos[i]) and set up the matrix
// for the multiplication
//----------------------------------------------------------------------
  ideal vars=var(vpos[1]);
  for(int i=2;i<=size(vpos);i++)
  {
    vars=vars,var(vpos[i]);
  }
  matrix varm[1][size(vpos)]=vars;
  if (size(vpos) > nrows(M))
  {
    matrix Mt[size(vpos)][ncols(M)];
    Mt[1..nrows(M),1..ncols(M)]=M;
    kill M;
    matrix M=Mt;
  }
//----------------------------------------------------------------------
// define the desired ideal
//----------------------------------------------------------------------
  ideal ret=vars^2,varm*M;
  return(ret);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(e(1),e(2),x,y,z),(dp(2),ds(3));
  module mo=x*gen(1)+y*gen(2);
  intvec iv=2,1;
  mod2id(mo,iv);
}
////////////////////////////////////////////////////////////////////////

proc id2mod(ideal i,intvec vpos)
"USAGE:     id2mod(I,vpos); I ideal, vpos intvec
NOTE:      * use this procedure only makes sense if the ideal contains
             all var(vpos[i])*var(vpos[j]) as monomial generators and
             all other generators are linear combinations of the
             var(vpos[i]) over the ring in the other variables
           * this is the inverse procedure to mod2id and should be applied
             only to ideals created by mod2id using the same intvec vpos
             (possibly after a standard basis computation)
RETURN:    module corresponding to the ideal by replacing var(vpos[i]) by
           gen(i) and omitting all generators var(vpos[i])*var(vpos[j])
EXAMPLE:   example id2mod; shows an example"
{
  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
//---------------------------------------------------------------------
// Initialization
//---------------------------------------------------------------------
  int n=size(i);
  int v=size(vpos);
  matrix tempmat;
  matrix mm[v][n];
//---------------------------------------------------------------------
// Conversion
//---------------------------------------------------------------------
  for(int j=1;j<=v;j++)
  {
    tempmat=coeffs(i,var(vpos[j]));
    mm[j,1..n]=tempmat[2,1..n];
  }
  for(j=1;j<=v;j++)
  {
    mm=subst(mm,var(vpos[j]),0);
  }
  module ret=simplify(mm,10);
  return(ret);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(e(1),e(2),x,y,z),(dp(2),ds(3));
  ideal i=e(2)^2,e(1)*e(2),e(1)^2,e(1)*y+e(2)*x;
  intvec iv=2,1;
  id2mod(i,iv);
}
///////////////////////////////////////////////////////////////////////

proc subrInterred(ideal mon, ideal sm, intvec iv)
"USAGE:   subrInterred(mon,sm,iv);
         sm:   ideal in a ring r with n + s variables,
               e.g. x_1,..,x_n and t_1,..,t_s
         mon:  ideal with monomial generators (not divisible by
               one of the t_i) such that sm is contained in the module
               k[t_1,..,t_s]*mon[1]+..+k[t_1,..,t_s]*mon[size(mon)]
         iv:   intvec listing the variables which are supposed to be used
               as x_i
RETURN:  list l:
          l[1]=the monomials from mon in the order used
          l[2]=their coefficients after interreduction
          l[3]=l[1]*l[2]
          (interreduced system of generators of sm seen as a submodule
          of k[t_1,..,t_s]*mon[1]+..+k[t_1,..,t_s]*mon[size(mon)])
EXAMPLE: example subrInterred; shows an example"
{
  int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
//-----------------------------------------------------------------------
// check that mon is really generated by monomials
// and sort its generators with respect to the monomial ordering
//-----------------------------------------------------------------------
  int err;
  for(int i=1;i<=ncols(mon);i++)
  {
    if ( size(mon[i]) > 1 )
    {
      err=1;
    }
  }
  if (err==1)
  {
    ERROR("mon has to be generated by monomials");
  }
  intvec sv=sortvec(mon);
  ideal mons;
  for(i=1;i<=size(sv);i++)
  {
    mons[i]=mon[sv[i]];
  }
  ideal itemp=maxideal(1);
  for(i=1;i<=size(iv);i++)
  {
    itemp=subst(itemp,var(iv[i]),0);
  }
  itemp=simplify(itemp,10);
  def r=basering;
  string tempstr="ring rtemp=" + charstr(basering) + ",(" + string(itemp)
                               + "),(C,dp);";
//-----------------------------------------------------------------------
// express m in terms of the generators of mon and check whether m
// can be considered as a submodule of k[t_1,..,t_n]*mon
//-----------------------------------------------------------------------
  module motemp;
  motemp[ncols(sm)]=0;
  poly ptemp;
  int j;
  for(i=1;i<=ncols(mons);i++)
  {
    for(j=1;j<=ncols(sm);j++)
    {
      ptemp=sm[j]/mons[i];
      motemp[j]=motemp[j]+ptemp*gen(i);
    }
  }
  for(i=1;i<=size(iv);i++)
  {
    motemp=subst(motemp,var(iv[i]),0);
  }
  matrix monmat[1][ncols(mons)]=mons;
  ideal dummy=monmat*motemp;
  for(i=1;i<=size(sm);i++)
  {
    if(sm[i]-dummy[i]!=0)
    {
      ERROR("the second argument is not a submodule of the assumed structure");
    }
  }
//----------------------------------------------------------------------
// do the interreduction and convert back
//----------------------------------------------------------------------
  execute(tempstr);
  def motemp=imap(r,motemp);
  motemp=interred(motemp);
  setring r;
  kill motemp;
  def motemp=imap(rtemp,motemp);
  list ret=monmat,motemp,monmat*motemp;
  for(i=1;i<=ncols(ret[2]);i++)
  {
    ret[2][i]=cleardenom(ret[2][i]);
  }
  for(i=1;i<=ncols(ret[3]);i++)
  {
    ret[3][i]=cleardenom(ret[3][i]);
  }
  return(ret);
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(x,y,z),dp;
  ideal i=x^2+x*y^2,x*y+x^2*y,z;
  ideal j=x^2+x*y^2,x*y,z;
  ideal mon=x^2,z,x*y;
  intvec iv=1,3;
  subrInterred(mon,i,iv);
  subrInterred(mon,j,iv);
}
////////////////////////////////////////////////////////////////////////