source: git/Singular/LIB/poly.lib @ d2b2a7

// $Id: poly.lib,v 1.12 1998-05-05 11:55:33 krueger Exp $
//system("random",787422842);
//(GMG, last modified 22.06.96)
//(obachman: 17.12.97 -- added katsura)
///////////////////////////////////////////////////////////////////////////////

version="$Id: poly.lib,v 1.12 1998-05-05 11:55:33 krueger Exp $";
info="
LIBRARY:  poly.lib      PROCEDURES FOR MANIPULATING POLYS, IDEALS, MODULES

 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
          (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 newly created and set ring
                     (32003, x(0..n), dp)
        katsura()  : katsura ideal of basering 
EXAMPLE: example katsura; shows examples
"
{
  if ( size(#) == 1 && typeof(#[1]) == "int")
  {
    ring katsura_ring = 32003, x(0..#[1]), dp;
    keepring katsura_ring;
  }
  ideal s;
  int i, j;
  int n = nvars(basering) -1;
  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;
  katsura();
  katsura(3);
}

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)
         //* Zur Zeit noch ein Bug, da erste Bettizahl falsch berechnet wird:
         //betti(0) ist -1 statt 0
EXAMPLE: example freerank; shows examples
"
{
  int rk;
  def M = simplify(#[1],10);
  list mre = mres(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 -1, coker(M) is not free
                         // [1] should be 1, coker(syz(M))=M is free of rank 1
  freerank(syz (M),"");  // [2] should be gen(2)+gen(1) (minimal relation of M)
  freerank(i);
  freerank(syz(i));      //* bug, 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 div  2)-1);
   //----------------------- "quick search" for maxdeg ------------------------
      while ( (c-jet(c,i)==0)*(c-jet(c,i-1)!=0) == 0)
      {
         i = (o+1+u) div  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 div  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 div  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) div  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 div  2)-1);
   //----------------------- "quick search" for mindeg ------------------------
      while ( (jet(c,u)==0)*(jet(c,o)!=0) )
      {
         i = (o+u) div  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 div  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 div  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) div  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);
}