singular-data 4.0.3+ds-1 / usr / share / singular / LIB / random.lib

///////////////////////////////////////////////////////////////////////////
version="version random.lib 4.0.0.0 Jun_2013 "; // $Id: 9f9b7db9c658deb04a5d368847264ad29af5405f $
category="General purpose";
info="
LIBRARY:  random.lib    Creating Random and Sparse Matrices, Ideals, Polys

PROCEDURES:
 genericid(i[,p,b]);     generic sparse linear combinations of generators of i
 randomid(id,[k,b]);     random linear combinations of generators of id
 randommat(n,m[,id,b]);  nxm matrix of random linear combinations of id
 sparseid(k,u[,o,p,b]);  ideal of k random sparse poly's of degree d [u<=d<=o]
 sparsematrix(n,m,o[,.]);nxm sparse matrix of polynomials of degree<=o
 sparsemat(n,m[,p,b]);   nxm sparse integer matrix with random coefficients
 sparsepoly(u[,o,p,b]);  random sparse polynomial with terms of degree in [u,o]
 sparsetriag(n,m[,.]);   nxm sparse lower-triag intmat with random coefficients
 sparseHomogIdeal(k,u,[,.]); ideal with k sparse homogeneous generators of degree in [u, o]
 triagmatrix(n,m,o[,.]); nxm sparse lower-triag matrix of poly's of degree<=o
 randomLast(b);          random transformation of the last variable
 randomBinomial(k,u,..); binomial ideal, k random generators of degree >=u
           (parameters in square brackets [] are optional)
";

LIB "inout.lib";
LIB "general.lib";
LIB "matrix.lib";
///////////////////////////////////////////////////////////////////////////////

proc genericid (def id, list #)
"USAGE:   genericid(id[,p,b]);  id ideal/module, p,b integers
RETURN:  system of generators of id which are generic, sparse, triagonal linear
         combinations of given generators with coefficients in [1,b] and
         sparsety p percent, bigger p being sparser (default: p=75, b=30000)
NOTE:    For performance reasons try small bound b in characteristic 0
EXAMPLE: example genericid; shows an example
"
{
//----------------------------- set defaults ----------------------------------
   if( size(#)>=2 ) { int p=#[1]; int b=#[2];}
   if( size(#)==1 ) { int p=#[1]; int b=30000;}
   if( size(#)==0 ) { int p=75; int b=30000;}
//---------------- use sparsetriag for creation of genericid ------------------
   def i = simplify(id,10);
   i = i*sparsetriag(ncols(i),ncols(i),p,b);
   return(i);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(t,x,y,z),ds;
   ideal i= x3+y4,z4+yx,t+x+y+z;
   genericid(i,0,10);
   module m=[x,0,0,0],[0,y2,0,0],[0,0,z3,0],[0,0,0,t4];
   print(genericid(m));
}
///////////////////////////////////////////////////////////////////////////////

proc randomid (def id, list #)
"USAGE:   randomid(id[,k,b]);  id ideal/module, b,k integers
RETURN:  ideal/module having k generators which are random linear combinations
         of generators of id with coefficients in the interval [-b,b]
         (default: b=30000, k=size(id))
NOTE:    For performance reasons try small bound b in characteristic 0
EXAMPLE: example randomid; shows an example
"
{
//----------------------------- set defaults ----------------------------------
   if( size(#)>=2 ) { int k=#[1]; int b=#[2]; }
   if( size(#)==1 ) { int k=#[1]; int b=30000; }
   if( size(#)==0 ) { int k=size(id); int b=30000; }
//--------------------------- create randomid ---------------------------------
   def i = id;
   i = matrix(id)*random(b,ncols(id),k);
   return(i);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),dp;
   randomid(maxideal(2),2,9);
   module m=[x,0,1],[0,y2,0],[y,0,z3];
   show(randomid(m));
}
///////////////////////////////////////////////////////////////////////////////

proc randommat (int n, int m, list #)
"USAGE:   randommat(n,m[,id,b]);  n,m,b integers, id ideal
RETURN:  nxm matrix, entries are random linear combinations of elements
         of id and coefficients in [-b,b]
         [default: (id,b) = (maxideal(1),30000)]
NOTE:    For performance reasons try small bound b in char 0
EXAMPLE:  example randommat; shows an example
"
{
//----------------------------- set defaults ----------------------------------
   if( size(#)>=2 ) { ideal id=#[1]; int b=#[2]; }
   if( size(#)==1 ) { ideal id=#[1]; int b=30000; }
   if( size(#)==0 ) { ideal id=maxideal(1); int b=30000; }
//--------------------------- create randommat --------------------------------
   id=simplify(id,2);
   int g=ncols(id);
   matrix rand[n][m]; matrix ra[1][m];
   for (int k=1; k<=n; k=k+1)
   {
      ra = id*random(b,g,m);
      rand[k,1..m]=ra[1,1..m];
   }
   return(rand);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),dp;
   matrix A=randommat(3,3,maxideal(2),9);
   print(A);
   A=randommat(2,3);
   print(A);
}
///////////////////////////////////////////////////////////////////////////////

proc sparseid (int k, int u, list #)
"USAGE:   sparseid(k,u[,o,p,b]);  k,u,o,p,b integers
RETURN:  ideal having k generators, each of degree d, u<=d<=o, p percent of
         terms in degree d are 0, the remaining have random coefficients
         in the interval [1,b], (default: o=u, p=75, b=30000)
EXAMPLE: example sparseid; shows an example
"
{
//----------------------------- set defaults ----------------------------------
   if( size(#)>=3 ) { int o=#[1]; int p=#[2]; int b=#[3]; }
   else {if( size(#)==2 ) { int o=#[1]; int p=#[2]; int b=30000; }
   else {if( size(#)==1 ) { int o=#[1]; int p=75; int b=30000; }
   else {if( size(#)==0 ) { int o=u; int p=75; int b=30000; }}}}
//------------------ use sparsemat for creation of sparseid -------------------
   int ii; matrix i[1][k]; intmat m;
   if( u <=0 )
   {
      m = sparsemat(1,k,p,b);
      i = m;
      u=1;
   }
   for ( ii=u; ii<=o; ii++)
   {
       m = sparsemat(size(maxideal(ii)),k,p,b);
       i = i+matrix(maxideal(ii))*m;
   }
   return(ideal(i));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(a,b,c,d),ds;
   sparseid(2,3);"";
   sparseid(3,0,4,90,9);
}

///////////////////////////////////////////////////////////////////////////////
proc sparseHomogIdeal (int k, int u, list #)
"USAGE:   sparseid(k,u[,o,p,b]);  k,u,o,p,b integers
RETURN:  ideal having k homogeneous generators, each of random degree in the
         interval [u,o], p percent of terms in degree d are 0, the remaining
         have random coefficients in the interval [1,b], (default: o=u, p=75,
         b=30000)
EXAMPLE: example sparseid; shows an example
"
{
//----------------------------- set defaults ----------------------------------
   if( size(#)>=3 ) { int o=#[1]; int p=#[2]; int b=#[3]; }
   if( size(#)==2 ) { int o=#[1]; int p=#[2]; int b=30000; }
   if( size(#)==1 ) { int o=#[1]; int p=75; int b=30000; }
   if( size(#)==0 ) { int o=u; int p=75; int b=30000; }
//------------------ use sparsemat for creation of sparseid -------------------
   int ii; ideal i; intmat m; ideal id;

   for ( ii=k; ii>0; ii--)
   {
       id = maxideal(random(u, o)); // monomial basis of some degree
       m = sparsemat(size(id),1,p,b); // random coefficients
       i[ii] = (matrix(id)*m)[1,1];
   }
   return(i);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(a,b,c,d),dp;
   sparseHomogIdeal(2,3);"";
   sparseHomogIdeal(3,0,4,90,9);
}


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

proc sparsemat (int n, int m, list #)
"USAGE:   sparsemat(n,m[,p,b]);  n,m,p,b integers
RETURN:  nxm integer matrix, p percent of the entries are 0, the remaining
         are random coefficients >=1 and <= b; [defaults: (p,b) = (75,1)]
EXAMPLE: example sparsemat; shows an example
"
{
   int r,h,ii;
   int t = n*m;
   intmat v[1][t];
//----------------------------- set defaults ----------------------------------
   if( size(#)>=2 ) { int p=#[1]; int b=#[2]; }
   if( size(#)==1 ) { int p=#[1]; int b=1; }
   if( size(#)==0 ) { int p=75; int b=1; }
//------------------------- check trivial cases ------------------------------
   if( p<0 ) { p = 0; }
   if(p>100) { p=100; }
//--------------- this is faster for not very sparse matrices ----------------
   if( p<40 )
   {
      for( ii=1; ii<=t; ii++ )
      { r=( random(1,100)>p ); v[1,ii]=r*random(1,b); h=h+r; }
   }
  int bb = t*(100-p);
  if( 100*h > bb )
  {
     while( 100*h > bb )
     { r=random(1,t); h=h-( v[1,r]>0 ); v[1,r]=0; }
  }
  else
  {
//------------------- this is faster for sparse matrices ---------------------
     while ( 100*h < bb )
     { r=random(1,t); h=h+(v[1,r]==0); v[1,r]=random(1,b); }
  }
  intmat M[n][m] = v[1,1..t];
  return(M);
}
example
{ "EXAMPLE:"; echo = 2;
   sparsemat(5,5);"";
   sparsemat(5,5,95);"";
   sparsemat(5,5,5);"";
   sparsemat(5,5,50,100);
}
///////////////////////////////////////////////////////////////////////////////
proc sparsematrix (int n, int m, int o, list #)
"USAGE:   sparsematrix(n,m,o[,u,pe,pp,b]);  n,m,o,u,pe,pp,b integers
RETURN:  nxm matrix, about pe percent of the entries are 0, the remaining
         are random polynomials of degree d, u<=d<=o, with  pp percent of
         the terms being 0, the remaining have random coefficients
         in the interval [1,b] [default: (pe,u,pp,b) = (0,50,75,100)]
EXAMPLE: example sparsematrix; shows an example
"
{
   int ii,jj;
   ideal id;
   matrix M[n][m];
//----------------------------- set defaults ----------------------------------
   int pe=50;int u=0;int pp=75;int b=100;
   if( size(#)==4 ) { u=#[1]; pe=#[2]; pp=#[3]; b=#[4]; }
   else { if( size(#)==3 ) { u=#[1]; pe=#[2]; pp=#[3]; }
   else { if( size(#)==2 ) { u=#[1]; pe=#[2]; }
   else {if( size(#)==1 ) { u=#[1]; }}}}
//------------------- use sparsemat and sparseid  -----------------------------
   intmat I = sparsemat(n,m,pe,1);
   for(ii=n; ii>0;ii--)
   {
     id = sparseid(m,u,o,pp,b);
     for(jj=m; jj>0; jj--)
     {
        if( I[ii,jj] !=0)
        {
          M[ii,jj]=id[jj];
        }
     }
   }
   return(M);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(a,b,c,d),dp;
   // sparse matrix of sparse polys of degree <=2:
   print(sparsematrix(3,4,2));"";
   // dense matrix of sparse linear forms:
   print(sparsematrix(3,3,1,1,0,55,9));
}
///////////////////////////////////////////////////////////////////////////////

proc sparsepoly (int u, list #)
"USAGE:   sparsepoly(u[,o,p,b]);  u,o,p,b integers
RETURN:  poly having only terms in degree d, u<=d<=o, p percentage of the terms
         in degree d are 0, the remaining have random coefficients in [1,b),
         (defaults: o=u, p=75, b=30000)
EXAMPLE:  example sparsepoly; shows an example
"
{
//----------------------------- set defaults ----------------------------------
   if( size(#)>=3 ) { int o=#[1]; int p=#[2]; int b=#[3]; }
   else {if( size(#)==2 ) { int o=#[1]; int p=#[2]; int b=30000; }
   else {if( size(#)==1 ) { int o=#[1]; int p=75; int b=30000; }
   else {if( size(#)==0 ) { int o=u; int p=75; int b=30000; }}}}
   int ii; poly f;
//----------------- use sparseid for creation of sparsepoly -------------------
   for( ii=u; ii<=o; ii++ ) { f=f+sparseid(1,ii,ii,p,b)[1]; }
   return(f);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),dp;
   sparsepoly(5);"";
   sparsepoly(3,5,90,9);
}
///////////////////////////////////////////////////////////////////////////////

proc sparsetriag (int n, int m, list #)
"USAGE:   sparsetriag(n,m[,p,b]);  n,m,p,b integers
RETURN:  nxm lower triagonal integer matrix, diagonal entries equal to 1, about
         p percent of lower diagonal entries are 0, the remaining are random
         integers >=1 and <= b; [defaults: (p,b) = (75,1)]
EXAMPLE: example sparsetriag; shows an example
"
{
   int ii,min,l,r; intmat M[n][m];
   int t=(n*(n-1)) div 2;
//----------------------------- set defaults ----------------------------------
   if( size(#)>=2 ) { int p=#[1]; int b=#[2]; }
   if( size(#)==1 ) { int p=#[1]; int b=1; }
   if( size(#)==0 ) { int p=75; int b=1; }
//---------------- use sparsemat for creation of sparsetriag ------------------
   intmat v[1][t]=sparsemat(1,t,p,b);
   if( n<=m ) { min=n-1; M[n,n]=1; }
   else { min=m; }
   for( ii=1; ii<=min; ii++ )
   {
      l=r+1; r=r+n-ii;
      M[ii..n,ii]=1,v[1,l..r];
   }
   return(M);
}
example
{ "EXAMPLE:"; echo = 2;
   sparsetriag(5,7);"";
   sparsetriag(7,5,90);"";
   sparsetriag(5,5,0);"";
   sparsetriag(5,5,50,100);
}
///////////////////////////////////////////////////////////////////////////////
proc triagmatrix (int n, int m, int o, list #)
"USAGE:   triagmatrix(n,m,o[,u,pe,pp,b]);  n,m,o,u,pe,pp,b integers
RETURN:  nxm lower triagonal matrix, diagonal entries equal to 1, about
         p percent of lower diagonal entries are 0, the remaining
         are random polynomials of degree d, u<=d<=o, with  pp percent of
         the terms being 0, the remaining have random coefficients
         in the interval [1,b] [default: (pe,u,pp,b) = (0,50,75,100)]
EXAMPLE: example triagmatrix; shows an example
"
{
   int ii,jj;
   ideal id;
   matrix M[n][m];
//----------------------------- set defaults ----------------------------------
   int pe=50;int u=0;int pp=75;int b=100;
   if( size(#)==4 ) { u=#[1]; pe=#[2]; pp=#[3]; b=#[4]; }
   if( size(#)==3 ) { u=#[1]; pe=#[2]; pp=#[3]; }
   if( size(#)==2 ) { u=#[1]; pe=#[2]; }
   if( size(#)==1 ) { u=#[1]; }
//------------------- use sparsemat and sparseid  -----------------------------
   intmat I = sparsetriag(n,m,pe,1);
   for(ii=1; ii<=n;ii++)
   {
     id = sparseid(m,u,o,pp,b);
     for(jj=1; jj<ii; jj++)
     {
        if( I[ii,jj] !=0)
        {
          M[ii,jj]=id[jj];
        }
     }
   }
   for(ii=1; ii<=n;ii++)
   {
      M[ii,ii]=1;
   }
   return(M);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r = 0,(a,b,c,d),dp;
   // sparse triagonal matrix of sparse polys of degree <=2:
   print(triagmatrix(3,4,2));"";
   // dense triagonal matrix of sparse linear forms:
   print(triagmatrix(3,3,1,1,0,55,9));
}
///////////////////////////////////////////////////////////////////////////////

proc randomLast(int b)
"USAGE:   randomLast(b);  b int
RETURN:  ideal = maxideal(1), but the last variable is exchanged by a random
         linear combination of all variables, with coefficients in the
         interval [-b,b], except for the last variable which always has
         coefficient 1
EXAMPLE: example randomLast; shows an example
"
{
  ideal i=maxideal(1);
  int k=size(i);
  if (k<=1) { return(i);}
  i[k]=0;
  i=randomid(i,size(i),b);
  ideal ires=maxideal(1);
  ires[k]=i[1]+var(k);
  return(ires);
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 0,(x,y,z),lp;
   ideal i = randomLast(10);
   i;
}
///////////////////////////////////////////////////////////////////////////////

proc randomBinomial(int k, int u, list #)
"USAGE:   randomBinomial(k,u[,o,b]);  k,u,o,b integers
RETURN:  binomial ideal, k homogeneous generators of degree d, u<=d<=o, with
         randomly chosen monomials and coefficients in the interval [-b,b]
         (default: u=o, b=10).
EXAMPLE: example randomBinomial; shows an example
"
{
//----------------------------- set defaults ----------------------------------
   if( size(#)>=2 ) { int o=#[1]; int b=#[2]; }
   if( size(#)==1 ) { int o=#[1]; int b=10; }
   if( size(#)==0 ) { int o=u; int b=10; }
//------------------ use sparsemat for creation of sparseid -------------------
   ideal i,m;
   int ii,jj,s,r1,r2;
   if ( o<u ) { o=u; }
   int a = k div (o-u+1);
   int c = k mod (o-u+1);
   for ( ii = u; ii<=o; ii++ )
   { m = maxideal(ii);
     s = size(m);
     for ( jj=1; jj<=a; jj++ )
     { r1 = random(-b,b);
       r1 = r1 + (r1==0)*random(1,b);
       r2 = random(-b,b);
       r2 = r2 + (r2==0)*random(-b,-1);
       i = i,r1*m[random(1,s div 2)] + r1*m[random(s div 2+1,s)];
       if ( ii < c+u )
       {  r1 = random(-b,b);
          r1 = r1 + (r1==0)*random(1,b);
          r2 = random(-b,b);
          r2 = r2 + (r2==0)*random(-b,-1);
          i = i,r1*m[random(1,s div 2)] + r2*m[random(s div 2+1,s)];
       }
     }
   }
   i = i[2..k+1];
   return(i);
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 0,(x,y,z),lp;
   ideal i = randomBinomial(4,5,6);
   i;
}
///////////////////////////////////////////////////////////////////////////////