Index: Singular/LIB/general.lib
===================================================================
--- Singular/LIB/general.lib	(revision 3d124a782bfbdf2088f156d5180ed7cb8a0b056c)
+++ Singular/LIB/general.lib	(revision 3d124a782bfbdf2088f156d5180ed7cb8a0b056c)
@@ -0,0 +1,591 @@
+// $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);
+}
+///////////////////////////////////////////////////////////////////////////////
