source: git/Singular/LIB/general.lib @ 3d124a7

// $Id: general.lib,v 1.1.1.1 1997-04-25 15:13:26 obachman Exp $
//system("random",787422842);
//(GMG, last modified 22.06.96)
///////////////////////////////////////////////////////////////////////////////

LIBRARY:  general.lib   PROCEDURES OF GENERAL TYPE

 A_Z("a",n);            string a,b,... of n comma seperated letters
 binomial(n,m[,../..]); n choose m (type int), [type string/type number]
 factorial(n[,../..]);  n factorial (=n!) (type int), [type string/number]
 fibonacci(n[,p]);      nth Fibonacci number [char p]
 kmemory();             int = active memory (kilobyte)
 killall();             kill all user-defined variables
 number_e(n);           compute exp(1) up to n decimal digits
 number_pi(n);          compute pi (area of unit circle) up to n digits
 primes(n,m);           intvec of primes p, n<=p<=m
 product(../..[,v]);    multiply components of vector/ideal/...[indices v]
 ringweights(r);        intvec of weights of ring variables of ring r
 sort(ideal/module);    sort generators according to monomial ordering
 sum(vector/id/..[,v]); add components of vector/ideal/...[with indices v]
           (parameters in square brackets [] are optional)

LIB "inout.lib";
///////////////////////////////////////////////////////////////////////////////

proc A_Z (string s,int n)
USAGE:   A_Z("a",n);  a any letter, n integer (-26<= n <=26, !=0)
RETURN:  string of n small (if a is small) or capital (if a is capital)
         letters, comma seperated, beginning with a, in alphabetical
         order (or revers alphabetical order if n<0)
EXAMPLE: example A_Z; shows an example
{
  if ( n>=-26 and n<=26 and n!=0 )
  {
    string alpha =
    "a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,"+
    "a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,"+
    "A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P,Q,R,S,T,U,V,W,X,Y,Z,"+
    "A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P,Q,R,S,T,U,V,W,X,Y,Z";
    int ii; int aa;
    for(ii=1; ii<=51; ii=ii+2)
    {
       if( alpha[ii]==s ) { aa=ii; }
    }
    if ( aa==0)
    {
      for(ii=105; ii<=155; ii=ii+2)
      {
        if( alpha[ii]==s ) { aa=ii; }
      }
    }
    if( aa!=0 )
    {
      string out;
      if (n > 0) { out = alpha[aa,2*(n)-1];  return (out); }
      if (n < 0)
      {
        string beta =
        "z,y,x,w,v,u,t,s,r,q,p,o,n,m,l,k,j,i,h,g,f,e,d,c,b,a,"+
        "z,y,x,w,v,u,t,s,r,q,p,o,n,m,l,k,j,i,h,g,f,e,d,c,b,a,"+
        "Z,Y,X,W,V,U,T,S,R,Q,P,O,N,M,L,K,J,I,H,G,F,E,D,C,B,A,"+
        "Z,Y,X,W,V,U,T,S,R,Q,P,O,N,M,L,K,J,I,H,G,F,E,D,C,B,A";
        if ( aa < 52 ) { aa=52-aa; }
        if ( aa > 104 ) { aa=260-aa; }
        out = beta[aa,2*(-n)-1]; return (out);
      }
    }
  }
}
example
{ "EXAMPLE:"; echo = 2;
   A_Z("c",5);
   A_Z("Z",-5);
   string sR = "ring R = (0,"+A_Z("A",6)+"),("+A_Z("a",10)+"),dp;";
   sR;
   execute sR;
   R;
}
///////////////////////////////////////////////////////////////////////////////

proc binomial (int n, int k, list #)
USAGE:   binomial(n,k[,p/s]); n,k,p integers, s string
RETURN:  binomial(n,k);    binomial coefficient n choose k of type int
                           (machine integer, limited size! )
         binomial(n,k,p);  n choose k in char p of type string
         binomial(n,k,s);  n choose k of type number (s any string), computed
                           in char of basering if a basering is defined
EXAMPLE: example binomial; shows an example
{
   if ( size(#)==0 ) { int rr=1; }
   if ( typeof(#[1])=="int") { ring bin = #[1],x,dp; number rr=1; }
   if ( typeof(#[1])=="string") { number rr=1; }
   if ( size(#)==0 or typeof(#[1])=="int" or typeof(#[1])=="string" )
   {
      def r = rr;
      if ( k<=0 or k>n ) { return((k==0)*r); }
      if ( k>n-k ) { k = n-k; }
      int l;
      for (l=1; l<=k; l=l+1 )
      {
         r=r*(n+1-l)/l;
      }
      if ( typeof(#[1])=="int" ) { return(string(r)); }
      return(r);
   }
}
example
{ "EXAMPLE:"; echo = 2;
   int b1 = binomial(10,7); b1;
   binomial(37,17,0);
   ring t = 31,x,dp;
   number b2 = binomial(37,17,""); b2;
}
///////////////////////////////////////////////////////////////////////////////

proc factorial (int n, list #)
USAGE:   factorial(n[,string]);  n integer
RETURN:  factorial(n); string of n! in char 0
         factorial(n,s);  n! of type number (s any string), computed in char of
         basering if a basering is defined
EXAMPLE: example factorial; shows an example
{
   if ( size(#)==0 ) { ring R = 0,x,dp; poly r=1; }
   if ( typeof(#[1])=="string" ) { number r=1; }
   if ( size(#)==0 or typeof(#[1])=="string" )
   {
      int l;
      for (l=2; l<=n; l=l+1)
      {
         r=r*l;
      }
      if ( size(#)==0 ) { return(string(r)); }
      return(r);
   }
}
example
{ "EXAMPLE:"; echo = 2;
   factorial(37);
   ring r1 = 32003,(x,y,z),ds;
   number p = factorial(37,""); p;
}
///////////////////////////////////////////////////////////////////////////////

proc fibonacci (int n, list #)
USAGE:   fibonacci(n[,string]);  (n integer)
RETURN:  fibonacci(n); string of nth Fibonacci number,
            f(0)=f(1)=1, f(i+1)=f(i-1)+f(i)
         fibonacci(n,s);  nth Fibonacci number of type number (s any string),
         computed in characteristic of basering if a basering is defined
EXAMPLE: example fibonacci; shows an example
{
   if ( size(#)==0 ) { ring fibo = 0,x,dp; number f=1; }
   if ( typeof(#[1])=="string" ) { number f=1; }
   if ( size(#)==0 or typeof(#[1])=="string" )
   {
      number g,h = 1,1; int ii;
      for (ii=3; ii<=n; ii=ii+1)
      {
         h = f+g; f = g; g = h;
      }
      if ( size(#)==0 ) { return(string(h)); }
      return(h);
   }
}
example
{ "EXAMPLE:"; echo = 2;
   fibonacci(37);
   ring r = 17,x,dp;
   number b = fibonacci(37,""); b;
}
///////////////////////////////////////////////////////////////////////////////

proc kmemory ()
USAGE:   kmemory();
RETURN:  memory used by active variables, of type int (in kilobyte)
EXAMPLE: example kmemory; shows an example
{
  if ( voice==2 ) { "// memory used by active variables (kilobyte):"; }
   return ((memory(0)+1023)/1024);
}
example
{ "EXAMPLE:"; echo = 2;
   kmemory();
}
///////////////////////////////////////////////////////////////////////////////

proc killall
USAGE:   killall(); (no parameter)
         killall("proc");
COMPUTE: killall(); kills all user-defined variables but not loaded procedures
         killall("proc"); kills all loaded procedures
RETURN:  no return value
NOTE:    killall should never be used inside a procedure
EXAMPLE: example killall; shows an example AND KILLS ALL YOUR VARIABLES
{
   list L=names(); int joni=size(L);
   if( size(#)==0 )
   {
      for ( ; joni>0; joni-- )
      {
         if( L[joni]!="LIB" and typeof(`L[joni]`)!="proc" ) { kill `L[joni]`; }
      }
      return();
   }
   if( #[1] == "proc" )
   {
      for ( joni=size(L); joni>0; joni-- )
      {
         if( L[joni]=="LIB" or typeof(`L[joni]`)=="proc" ) { kill `L[joni]`; }
      }
   }
}
example
{ "EXAMPLE:"; echo = 2;
   ring rtest; ideal i=x,y,z; number n=37; string str="hi";
   export rtest,i,n,str;     //this makes the local variables global
   listvar(all);
   killall();
   listvar(all);
}
///////////////////////////////////////////////////////////////////////////////

proc number_e (int n)
USAGE:   number_e(n);  n integer
COMPUTE: exp(1) up to n decimal digits (no rounding)
         by A.H.J. Sale's algorithm
RETURN:  - string of exp(1) if no basering of char 0 is defined;
         - exp(1), of type number, if a basering of char 0 is defined and
         display its decimal format
EXAMPLE: example number_e; shows an example
{
   int i,m,s,t;
   intvec u,e;
   u[n+2]=0; e[n+1]=0; e=e+1;
   if( defined(basering) )
   {
      if( char(basering)==0 ) { number r=2; t=1; }
   }
   string result = "2.";
   for( i=1; i<=n+1; i=i+1 )
   {
      e = e*10;
      for( m=n+1; m>=1; m=m-1 )
      {
         s    = e[m]+u[m+1];
         u[m] = s/(m+1);
         e[m] = s%(m+1);
      }
      result = result+string(u[1]);
      if( t==1 ) { r = r+number(u[1])/number(10)^i; }
   }
   if( t==1 ) { "//",result[1,n+1]; return(r); }
   return(result[1,n+1]);
}
example
{ "EXAMPLE:"; echo = 2;
   number_e(15);
   ring R = 0,t,lp;
   number e = number_e(10);
   e;
}
///////////////////////////////////////////////////////////////////////////////

proc number_pi (int n)
USAGE:   number_pi(n);  n positive integer
COMPUTE: pi (area of unit circle) up to n decimal digits (no rounding)
         by algorithm of S. Rabinowitz
RETURN:  - string of pi if no basering of char 0 is defined,
         - pi, of type number, if a basering of char 0 is defined and display
         its decimal format
EXAMPLE: example number_pi; shows an example
{
   int i,m,t,e,q,N;
   intvec r,p,B,Prelim;
   string result,prelim;
   N = (10*n)/3 + 2;
   p[N+1]=0; p=p+2; r=p;
   for( i=1; i<=N+1; i=i+1 ) { B[i]=2*i-1; }
   if( defined(basering) )
   {
      if( char(basering)==0 ) { number pi; number pri; t=1; }
   }
   for( i=0; i<=n; i=i+1 )
   {
      p = r*10;
      e = p[N+1];
      for( m=N+1; m>=2; m=m-1 )
      {
         r[m] = e%B[m];
         q    = e/B[m];
         e    = q*(m-1)+p[m-1];
      }
      r[1] = e%10;
      q    = e/10;
      if( q!=10 and q!=9 )
      {
         result = result+prelim;
         Prelim = q;
         prelim = string(q);
      }
      if( q==9 )
      {
         Prelim = Prelim,9;
         prelim = prelim+"9";
      }
      if( q==10 )
      {
         Prelim = (Prelim+1)-((Prelim+1)/10)*10;
         for( m=size(Prelim); m>0; m=m-1)
         {
            prelim[m] = string(Prelim[m]);
         }
         result = result+prelim;
         if( t==1 ) { pi=pi+pri; }
         Prelim = 0;
         prelim = "0";
      }
      if( t==1 ) { pi=pi+number(q)/number(10)^i; }
   }
   result = result,prelim[1];
   result = "3."+result[2,n-1];
   if( t==1 ) { "//",result; return(pi); }
   return(result);
}
example
{ "EXAMPLE:"; echo = 2;
   number_pi(5);
   ring r = 0,t,lp;
   number pi = number_pi(6);
   pi;
}
///////////////////////////////////////////////////////////////////////////////

proc primes (int n, int m)
USAGE:   primes(n,m);  n,m integers
RETURN:  intvec, consisting of all primes p, prime(n)<=p<=m, in increasing
         order if n<=m, resp. prime(m)<=p<=n, in decreasing order if m<n
NOTE:    prime(n); returns the biggest prime number <= n (if n>=2, else 2)
EXAMPLE: example primes; shows an example
{  int change;
   if ( n>m ) { change=n; n=m ; m=change; change=1; }
   int q,p = prime(m),prime(n); intvec v = q; q = q-1;
   while ( q>=p ) { q = prime(q); v = q,v; q = q-1; }
   if ( change==1 ) { v = v[size(v)..1]; }
   return(v);
}
example
{  "EXAMPLE:"; echo = 2;
   primes(50,100);
   intvec v = primes(37,1); v;
}
///////////////////////////////////////////////////////////////////////////////

proc product (id, list #)
USAGE:    product(id[,v]); id=ideal/vector/module/matrix
          resp.id=intvec/intmat, v=intvec (e.g. v=1..n, n=integer)
RETURN:   poly resp. int which is the product of all entries of id, with index
          given by v (default: v=1..number of entries of id)
NOTE:     id is treated as a list of polys resp. integers. A module m is
          identified with corresponding matrix M (columns of M generate m)
EXAMPLE:  example product; shows an example
{
   int n,j;
   if( typeof(id)=="poly" or typeof(id)=="ideal" or typeof(id)=="vector"
       or typeof(id)=="module" or typeof(id)=="matrix" )
   {
      ideal i = ideal(matrix(id));
      if( size(#)!=0 ) { i = i[#[1]]; }
      n = ncols(i); poly f=1;
   }
   if( typeof(id)=="int" or typeof(id)=="intvec" or typeof(id)=="intmat" )
   {
      if ( typeof(id) == "int" ) { intmat S =id; }
      else { intmat S = intmat(id); }
      intvec i = S[1..nrows(S),1..ncols(S)];
      if( size(#)!=0 ) { i = i[#[1]]; }
      n = size(i); int f=1;
   }
   for( j=1; j<=n; j=j+1 ) { f=f*i[j]; }
   return(f);
}
example
{  "EXAMPLE:"; echo = 2;
   ring r= 0,(x,y,z),dp;
   ideal m = maxideal(1);
   product(m);
   matrix M[2][3] = 1,x,2,y,3,z;
   product(M);
   intvec v=2,4,6;
   product(M,v);
   intvec iv = 1,2,3,4,5,6,7,8,9;
   v=1..5,7,9;
   product(iv,v);
   intmat A[2][3] = 1,1,1,2,2,2;
   product(A,3..5);
}
///////////////////////////////////////////////////////////////////////////////

proc ringweights (r)
USAGE:   ringweights(r); r ring
RETURN:  intvec of weights of ring variables. If, say, x(1),...,x(n) are the
         variables of the ring r, in this order, the resulting intvec is
         deg(x(1)),...,deg(x(n)) where deg denotes the weighted degree if
         the monomial ordering of r has only one block of type ws,Ws,wp or Wp.
NOTE:    In all other cases, in particular if there is more than one block,
         the resulting intvec is 1,...,1
EXAMPLE: example ringweights; shows an example
{
   int i; intvec v; setring r;
   for (i=1; i<=nvars(basering); i=i+1) { v[i] = deg(var(i)); }
   return(v);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r1=32003,(x,y,z),wp(1,2,3);
   ring r2=32003,(x,y,z),Ws(1,2,3);
   ring r=0,(x,y,u,v),lp;
   intvec vr=ringweights(r1); vr;
   ringweights(r2);
   ringweights(r);
}
///////////////////////////////////////////////////////////////////////////////

proc sort (id, list #)
USAGE:   sort(id[v,o,n]); id=ideal/module/intvec/list (of intvec's or int's)
         sort may be called with 1, 2 or 3 arguments in the following way:
         -  sort(id[v,n]); v=intvec, n=integer,
         -  sort(id[o,n]); o=string (any allowed ordstr of a ring), n=integer
RETURN:  a list of two elements:
         [1]: object of same type as input but sorted in the following manner:
           - if id=ideal/module: generators of id are sorted w.r.t. intvec v
             (id[v[1]] becomes 1-st, id[v[2]]  2-nd element, etc.). If no v is
             present, id is sorted w.r.t. ordering o (if o is given) or w.r.t.
             actual monomial ordering (if no o is given):
                    generators with smaller leading term come first
             (e.g. sort(id); sorts w.r.t actual monomial ordering)
           - if id=list of intvec's or int's: consider a list element, say
             id[1]=3,2,5, as exponent vector of the monomial x^3*y^2*z^5;
             the corresponding monomials are ordered w.r.t. intvec v (s.a.).
             If no v is present, the monomials are sorted w.r.t. ordering o
             (if o is given) or w.r.t. lexicographical ordering (if no o is
             given). The corresponding ordered list of exponent vectors is
             returned.
             (e.g. sort(id); sorts lexicographically, smaller int's come first)
             WARNING: Since negative exponents create the 0 plynomial in
             Singular, id should not contain negative integers: the result might
             not be as exspected
           - if id=intvec: id is treated as list of integers
           - if n!=0 the ordering is inverse, i.e. w.r.t. v(size(v)..1)
             default: n=0
         [2]: intvec, describing the permutation of the input (hence [2]=v if
             v is given)
NOTE:    If v is given, id may be any simply indexed object (e.g. any list);
         entries of v must be pairwise distinct to get a permutation if id.
         Zero generators of ideal/module are deleted
EXAMPLE: example sort; shows an example
{
   int ii,jj,s,n = 0,0,1,0;
   intvec v;
   if ( defined(basering) ) { def P = basering; }
   if ( size(#)==0 and (typeof(id)=="ideal" or typeof(id)=="module") )
   {
      id = simplify(id,2);
      for ( ii=1; ii<size(id); ii++ )
      {
         if ( id[ii]!=id[ii+1] ) { break;}
      }
      if ( ii != size(id) ) { v = sortvec(id); }
      else  { v = size(id)..1; }
      if ( v == 0 ) { v = 1; }
   }
   if ( size(#)>=1 and (typeof(id)=="ideal" or typeof(id)=="module") )
   {
      if ( typeof(#[1])=="string" )
      {
         execute "ring r1 =("+charstr(P)+"),("+varstr(P)+"),("+#[1]+");";
         def i = imap(P,id);
         v = sortvec(i);
         setring P;
         n=2;
      }
   }
   if ( typeof(id)=="intvec" or typeof(id)=="list" and n==0 )
   {
      string o;
      if ( size(#)==0 ) { o = "lp"; n=1; }
      if ( size(#)>=1 )
      {
         if ( typeof(#[1])=="string" ) { o = #[1]; n=1; }
      }
   }
   if ( typeof(id)=="intvec" or typeof(id)=="list" and n==1 )
   {
      if ( typeof(id)=="list" )
      {
         for (ii=1; ii<=size(id); ii++)
         {
            if (typeof(id[ii]) != "intvec" and typeof(id[ii]) != "int")
               { "// list elements must be intvec/int"; return(); }
            else
               { s=size(id[ii])*(s < size(id[ii])) + s*(s >= size(id[ii])); }
         }
      }
      execute "ring r=0,x(1..s),("+o+");";
      ideal i;
      poly f;
      for (ii=1; ii<=size(id); ii++)
      {
         f=1;
         for (jj=1; jj<=size(id[ii]); jj++)
         {
            f=f*x(jj)^(id[ii])[jj];
         }
         i[ii]=f;
      }
      v = sort(i)[2];
   }
   if ( size(#)!=0 and n==0 ) { v = #[1]; }
   if( size(#)==2 )
   {
      if ( #[2] != 0 ) { v = v[size(v)..1]; }
   }
   s = size(v);
   def m = id;
   for ( jj=1; jj<=s; jj=jj+1) { m[jj] = id[v[jj]]; }
   list L=m,v;
   return(L);
}
example
{  "EXAMPLE:"; echo = 2;
   ring r0 = 0,(x,y,z),lp;
   ideal i = x3,y3,z3,x2z,x2y,y2z,y2x,z2y,z2x,xyz;
   show(sort(i));"";
   show(sort(i,"wp(1,2,3)"));"";
   intvec v=10..1;
   show(sort(i,v));"";
   show(sort(i,v,1));"";   // should be the identity
   ring r1  = 0,t,ls;
   ideal j = t14,t6,t28,t20,t12,t34,t26,t18,t40,t32,t24,t38,t30,t36;
   show(sort(j)[1]);"";
   show(sort(j,"lp")[1]);"";
   list L =1,5..8,10,2,8..5,8,3..10;
   sort(L)[1];"";          // sort L lexicographically
   sort(L,"Dp",1)[1];      // sort L w.r.t (total sum, reverse lex)
}
///////////////////////////////////////////////////////////////////////////////

proc sum (id, list #)
USAGE:    sum(id[,v]); id=ideal/vector/module/matrix resp. id=intvec/intmat,
                       v=intvec (e.g. v=1..n, n=integer)
RETURN:   poly resp. int which is the sum of all entries of id, with index
          given by v (default: v=1..number of entries of id)
NOTE:     id is treated as a list of polys resp. integers. A module m is
          identified with corresponding matrix M (columns of M generate m)
EXAMPLE:  example sum; shows an example
{
   if( typeof(id)=="poly" or typeof(id)=="ideal" or typeof(id)=="vector"
       or typeof(id)=="module" or typeof(id)=="matrix" )
   {
      ideal i = ideal(matrix(id));
      if( size(#)!=0 ) { i = i[#[1]]; }
      matrix Z = matrix(i);
   }
   if( typeof(id)=="int" or typeof(id)=="intvec" or typeof(id)=="intmat" )
   {
      if ( typeof(id) == "int" ) { intmat S =id; }
      else { intmat S = intmat(id); }
      intvec i = S[1..nrows(S),1..ncols(S)];
      if( size(#)!=0 ) { i = i[#[1]]; }
      intmat Z=transpose(i);
   }
   intvec v; v[ncols(Z)]=0; v=v+1;
   return((Z*v)[1,1]);
}
example
{  "EXAMPLE:"; echo = 2;
   ring r= 0,(x,y,z),dp;
   vector pv = [xy,xz,yz,x2,y2,z2];
   sum(pv);
   //sum(pv,2..5);
   //matrix M[2][3] = 1,x,2,y,3,z;
   //sum(M);
   //intvec w=2,4,6;
   //sum(M,w);
   //intvec iv = 1,2,3,4,5,6,7,8,9;
   //w=1..5,7,9;
   //sum(iv,w);
   //intmat m[2][3] = 1,1,1,2,2,2;
   //sum(m,3..4);
}
///////////////////////////////////////////////////////////////////////////////