source: git/Singular/LIB/general.lib @ a9ca4a

//GMG, last modified 18.6.99
//anne, added deleteSublist and watchdog 12.12.2000
//eric, added absValue 11.04.2002
///////////////////////////////////////////////////////////////////////////////
version="$Id: general.lib,v 1.63 2009-07-10 18:52:05 motsak Exp $";
category="General purpose";
info="
LIBRARY:  general.lib   Elementary Computations of General Type

PROCEDURES:
 A_Z(\"a\",n);          string a,b,... of n comma separated letters
 ASCII([n,m]);          string of printable ASCII characters (number n to m)
 absValue(c);           absolute value of c
 binomial(n,m[,../..]); n choose m (type int), [type bigint]
 deleteSublist(iv,l);   delete entries given by iv from list l
 factorial(n[,../..]);  n factorial (=n!) (type int), [type bigint]
 fibonacci(n[,p]);      nth Fibonacci number [char p]
 kmemory([n[,v]]);      active [allocated] memory in 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]
 sort(ideal/module);    sort generators according to monomial ordering
 sum(vector/id/..[,v]); add components of vector/ideal/...[with indices v]
 watchdog(i,cmd);       only wait for result of command cmd for i seconds
 which(command);        search for command and return absolute path, if found
 primecoeffs(J[,q]);    primefactors <= min(p,32003) of coeffs of J
 primefactors(n[,p]);  primefactors <= min(p,32003) of n
 timeStd(i,d)           std(i) if the standard basis computation finished
                        after d-1 seconds and i otherwhise
 timeFactorize(p,d)     factorize(p) if the factorization finished after d-1
                        seconds otherwhise f is considered to be irreducible
 factorH(p)             changes variables to become the last variable the
                        principal one in the multivariate factorization and
                        factorizes then the polynomial

           (parameters in square brackets [] are optional)
";

LIB "inout.lib";
LIB "poly.lib";
LIB "matrix.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 separated, beginning with a, in alphabetical
         order (or reverse 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 ASCII (list #)
"USAGE:   ASCII([n,m]); n,m= integers (32 <= n <= m <= 126)
RETURN:   string of printable ASCII characters (no native language support)@*
          ASCII():    string of  all ASCII characters with its numbers,@*
          ASCII(n):   n-th ASCII character@*
          ASCII(n,m): n-th up to m-th ASCII character (inclusive)
EXAMPLE: example ASCII; shows an example
"
{
  string s1 =
 "     !    \"    #    $    %    &    '    (    )    *    +    ,    -    .
32   33   34   35   36   37   38   39   40   41   42   43   44   45   46
/    0    1    2    3    4    5    6    7    8    9    :    ;    <    =
47   48   49   50   51   52   53   54   55   56   57   58   59   60   61
>    ?    @    A    B    C    D    E    F    G    H    I    J    K    L
62   63   64   65   66   67   68   69   70   71   72   73   74   75   76
M    N    O    P    Q    R    S    T    U    V    W    X    Y    Z    [
77   78   79   80   81   82   83   84   85   86   87   88   89   90   91
\\    ]    ^    _    `    a    b    c    d    e    f    g    h    i    j
92   93   94   95   96   97   98   99  100  101  102  103  104  105  10
k    l    m    n    o    p    q    r    s    t    u    v    w    x    y
107  108  109  110  111  112  113  114  115  116  117  118  119  120  121
z    {    |    }    ~
122  123  124  125  126 ";

   string s2 =
 " !\"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\\]^_`abcdefghijklmnopqrstuvwxyz{|}~";

   if ( size(#) == 0 )
   {
      return(s1);
   }
   if ( size(#) == 1 )
   {
      return( s2[#[1]-31] );
   }
   if ( size(#) == 2 )
   {
      return( s2[#[1]-31,#[2]-#[1]+1] );
   }
}
example
{ "EXAMPLE:"; echo = 2;
   ASCII();"";
   ASCII(42);
   ASCII(32,126);
}
///////////////////////////////////////////////////////////////////////////////

proc absValue(def c)
"USAGE:  absValue(c); c int, number or poly
RETURN:  absValue(c); the absolute value of c
NOTE:    absValue(c)=c if c>=0; absValue=-c if c<=0.
@*       So the function can be applied to any type, for which comparison
@*       operators are defined.
EXAMPLE: example absValue; shows an example
"
{
  if (c>=0) { return(c); }
  else { return(-c); }
}
example
{ "EXAMPLE:"; echo = 2;
   ring r1 = 0,x,dp;
   absValue(-2002);

   poly f=-4;
   absValue(f);
}
///////////////////////////////////////////////////////////////////////////////

proc binomial (int n, int k)
"USAGE:   binomial(n,k); n,k integers
RETURN:  binomial(n,k); binomial coefficient n choose k
@*       - of type bigint (computed in characteristic 0)
NOTE:    In any characteristic, binomial(n,k) = coefficient of x^k in (1+x)^n
SEE ALSO: prime
EXAMPLE: example binomial; shows an example
"
{
   bigint l;
   bigint r=1;
   bigint kk=k;
   bigint nn=n;
   if ( k > n-k )
   { k = n-k; }
   if ( k<=0 or k>n )               //trivial cases
   { r = (k==0)*r; }
   for (l=1; l<=kk; l=l+1 )
   {
      r=r*(nn+1-l)/l;
   }
   return(r);
}
example
{ "EXAMPLE:"; echo = 2;
   binomial(200,100);"";                 //type bigint
   int n,k = 200,100;
   bigint b1 = binomial(n,k);
   ring r = 0,x,dp;
   poly b2 = coeffs((x+1)^n,x)[k+1,1];  //coefficient of x^k in (x+1)^n
   b1-b2;                               //b1 and b2 should coincide
}
///////////////////////////////////////////////////////////////////////////////

static proc binomp (int n, int k, int p)
 //computes binomial coefficient n choose k in basering of char p > 0
 //binomial(n,k) = coefficient of x^k in (1+x)^n.
 //Let n=q*p^j, gcd(q,p)=1, then (1+x)^n = (1 + x^(p^j))^q. We have
 //binomial(n,k)=0 if k!=l*p^j and binomial(n,l*p^j) = binomial(q,l).
 //Do this reduction first. Then, in denominator and numerator
 //of defining formula for binomial coefficient, reduce those factors
 //mod p which are not divisible by p and cancel common factors p. Hence,
 //if n = h*p+r, k=l*p+s, r,s<p, binomial(n,k) = binomial(r,s)*binomial(h,l)
{
   int l,q,i= 1,n,1;
   number r=1;
   if ( k > n-k )
   { k = n-k;
   }
   if ( k<=0 or k>n)               //trivial cases
   { r = (k==0)*r;
   }
   else
   {
      while ( q mod p == 0 )
      {  l = l*p;
         q = q div p;
      }                            //we have now n=q*l, l=p^j, gcd(q,p)=1;
      if (k mod l != 0 )
      { r = 0;
      }
      else
      {  l = k div l;
         n = q mod p;
         k = l mod p;              //now 0<= k,n <p, use binom0 for n,k
         q = q div p;              //recursion for q,l
         l = l div p;              //use binomp for q,l
         r = binom0(n,k)*binomp(q,l,p);
      }
   }
   return(r);
}
///////////////////////////////////////////////////////////////////////////////

proc factorial (int n)
"USAGE:    factorial(n);  n integer
RETURN:   factorial(n):   n! of type bigint.
SEE ALSO: prime
EXAMPLE:  example factorial; shows an example
"
{
   bigint r=1;
//------------------------------ computation --------------------------------
   for (int l=2; l<=n; l++)
   {
      r=r*l;
   }
   return(r);
}
example
{ "EXAMPLE:"; echo = 2;
   factorial(37);"";                 //37! (as long integer)
}
///////////////////////////////////////////////////////////////////////////////

proc fibonacci (int n)
"USAGE:    fibonacci(n);  n,p integers
RETURN:   fibonacci(n): nth Fibonacci number, f(0)=f(1)=1, f(i+1)=f(i-1)+f(i)
@*        - computed in characteristic 0, of type bigint
SEE ALSO: prime
EXAMPLE: example fibonacci; shows an example
"
{
  bigint f,g,h=1,1,1;
//------------------------------ computation --------------------------------
   for (ii=3; ii<=n; ii=ii+1)
   {
      h = f+g; f = g; g = h;
   }
   return(h);
}
example
{ "EXAMPLE:"; echo = 2;
   fibonacci(42); "";             //f(42) of type string (as long integer)
   ring r = 2,x,dp;
   number b = fibonacci(42,2);    //f(42) of type number, computed in r
   b;
}
///////////////////////////////////////////////////////////////////////////////

proc kmemory (list #)
"USAGE:   kmemory([n,[v]]); n,v integers
RETURN:  memory in kilobyte of type bigint
@*       n=0: memory used by active variables (same as no parameters)
@*       n=1: total memory allocated by Singular
DISPLAY: detailed information about allocated and used memory if v!=0
NOTE:    kmemory uses internal function 'memory' to compute kilobyte, and
         is the same as 'memory' for n!=0,1,2
EXAMPLE: example kmemory; shows an example
"
{
   int n;
   int verb;
   if (size(#) != 0)
   {
     n=#[1];
     if (size(#) >1)
     { verb=#[2]; }
   }

  if ( verb != 0)
  {
    if ( n==0)
    { dbprint(printlevel-voice+3,
      "// memory used, at the moment, by active variables (kilobyte):"); }
    if ( n==1 )
    { dbprint(printlevel-voice+3,
      "// total memory allocated, at the moment, by SINGULAR (kilobyte):"); }
  }
  return ((memory(n)+1023) div 1024);
}
example
{ "EXAMPLE:"; echo = 2;
   kmemory();
   kmemory(1,1);
}
///////////////////////////////////////////////////////////////////////////////

proc killall
"USAGE:   killall(); (no parameter)
         killall(\"type_name\");
         killall(\"not\", \"type_name\");
RETURN:  killall(); kills all user-defined variables except loaded procedures,
         no return value.
@*       - killall(\"type_name\"); kills all user-defined variables,
           of type \"type_name\"
@*       - killall(\"not\", \"type_name\"); kills all user-defined variables,
           except those of type \"type_name\" and except loaded procedures
@*       - killall(\"not\", \"name_1\", \"name_2\", ...);
           kills all user-defined variables, except those of name \"name_i\"
           and except loaded procedures
NOTE:    killall should never be used inside a procedure
EXAMPLE: example killall; shows an example AND KILLS ALL YOUR VARIABLES
"
{
  if (system("with","Namespaces"))
  {
    list @marie=names(Top);
  }
  else
  {
    list @marie=names();
  }
  int j, no_kill, @joni;
  for ( @joni=1; @joni<=size(#);  @joni++)
  {
    if (typeof(#[@joni]) != "string")
    {
      ERROR("Need string as " + string(@joni) + "th argument");
    }
  }

  // kills all user-defined variables but not loaded procedures
  if( size(#)==0 )
  {
    for ( @joni=size(@marie); @joni>0; @joni-- )
    {
      if( @marie[@joni]!="LIB" and typeof(`@marie[@joni]`)!="proc"
      and typeof(`@marie[@joni]`)!="package")
      { kill `@marie[@joni]`; }
    }
  }
  else
  {
    // kills all user-defined variables
    if( size(#)==1 )
    {
      // of type proc
      if( #[1] == "proc" )
      {
        for ( @joni=size(@marie); @joni>0; @joni-- )
        {
          if( (@marie[@joni]!="General")
          and (@marie[@joni]!="Top")
          and (@marie[@joni]!="killall")
          and (@marie[@joni]!="LIB") and
              ((typeof(`@marie[@joni]`)=="package") or
              (typeof(`@marie[@joni]`)=="proc")))
          {
            if (defined(`@marie[@joni]`)) {kill `@marie[@joni]`;}
          }
          if (!defined(@joni)) break;
        }
        if ((system("with","Namespaces")) && defined(General))
        {
          @marie=names(General);
          for ( @joni=size(@marie); @joni>0; @joni-- )
          {
            if(@marie[@joni]!="killall"
            and typeof(`@marie[@joni]`)=="proc")
            { kill General::`@marie[@joni]`; }
          }
          kill General::killall;
        }
      }
      else
      {
        // other types
        for ( @joni=size(@marie); @joni>2; @joni-- )
        {
          if(typeof(`@marie[@joni]`)==#[1] and @marie[@joni]!="LIB"
             and typeof(`@marie[@joni]`)!="proc")
          { kill `@marie[@joni]`; }
        }
      }
    }
    else
    {
      // kills all user-defined variables whose name or type is not #i
      for ( @joni=size(@marie); @joni>2; @joni-- )
      {
        if ( @marie[@joni] != "LIB" && @marie[@joni] != "Top"
             && typeof(`@marie[@joni]`) != "proc")
        {
          no_kill = 0;
          for (j=2; j<= size(#); j++)
          {
            if (typeof(`@marie[@joni]`)==#[j] or @marie[@joni] == #[j])
            {
              no_kill = 1;
              break;
            }
          }
          if (! no_kill)
          {
            kill `@marie[@joni]`;
          }
        }
        if (!defined(@joni)) break;
      }
    }
  }
}
example
{ "EXAMPLE:"; echo = 2;
   ring rtest; ideal i=x,y,z; string str="hi"; int j = 3;
   export rtest,i,str,j;       //this makes the local variables global
   listvar();
   killall("ring");            // kills all rings
   listvar();
   killall("not", "int");      // kills all variables except int's (and procs)
   listvar();
   killall();                  // kills all vars except loaded procs
   listvar();
}
///////////////////////////////////////////////////////////////////////////////

proc number_e (int n)
"USAGE:   number_e(n);  n integer
RETURN:  Euler number e=exp(1) up to n decimal digits (no rounding)
@*       - of type string if no basering of char 0 is defined
@*       - of type number if a basering of char 0 is defined
DISPLAY: decimal format of e if printlevel > 0 (default:printlevel=0 )
NOTE:    procedure uses algorithm of A.H.J. Sale
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 div (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 )
   { dbprint(printlevel-voice+2,"// "+result[1,n+1]);
     return(r);
   }
   return(result[1,n+1]);
}
example
{ "EXAMPLE:"; echo = 2;
   number_e(30);"";
   ring R = 0,t,lp;
   number e = number_e(30);
   e;
}
///////////////////////////////////////////////////////////////////////////////

proc number_pi (int n)
"USAGE:   number_pi(n);  n positive integer
RETURN:  pi (area of unit circle) up to n decimal digits (no rounding)
@*       - of type string if no basering of char 0 is defined,
@*       - of type number, if a basering of char 0 is defined
DISPLAY: decimal format of pi if printlevel > 0 (default:printlevel=0 )
NOTE:    procedure uses algorithm of S. Rabinowitz
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) div 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 div B[m];
         e    = q*(m-1)+p[m-1];
      }
      r[1] = e%10;
      q    = e div 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) div 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 )
   { dbprint(printlevel-voice+2,"// "+result);
     return(pi);
   }
   return(result);
}
example
{ "EXAMPLE:"; echo = 2;
   number_pi(11);"";
   ring r = (real,10),t,dp;
   number pi = number_pi(11); 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 <= min(n,32003)
         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/intvec/intmat/list,
          v intvec  (default: v=1..number of entries of id)
ASSUME:   list members can be multiplied.
RETURN:   The product of all entries of id [with index given by v] of type
          depending on the entries of id.
NOTE:     If id is not a list, id is treated as a list of polys resp. integers.
          A module m is identified with the corresponding matrix M (columns
          of M generate m).
@*        If v is outside the range of id, we have the empty product and the
          result will be 1 (of type int).
EXAMPLE:  example product; shows an example
"
{
//-------------------- initialization and special feature ---------------------
   int n,j,tt;
   string ty;                                //will become type of id
   list l;

// We wish to allow something like product(x(1..10)) if x(1),...,x(10) are
// variables. x(1..10) is a list of polys and enters the procedure with
// id=x(1) and # a list with 9 polys, #[1]= x(2),...,#[9]= x(10). Hence, in
// this case # is never empty. If an additional intvec v is given,
// it is added to #, so we have to separate it first and make
// the rest a list which has to be multiplied.

   int s = size(#);
   if( s!=0 )
   {  if ( typeof(#[s])=="intvec" or typeof(#[s])=="int")
      {
         intvec v = #[s];
         tt=1;
         s=s-1;
         if ( s>0 ) { # = #[1..s]; }
      }
   }
   if ( s>0 )
   {
      l = list(id)+#;
      kill id;
      list id = l;                                    //case: id = list
      ty = "list";
      n = size(id);
   }
   else
   {
      ty = typeof(id);
      if( ty == "list" )
      { n = size(id); }
   }
//------------------------------ reduce to 3 cases ---------------------------
   if( ty=="poly" or ty=="ideal" or ty=="vector"
       or ty=="module" or ty=="matrix" )
   {
      ideal i = ideal(matrix(id));
      kill id;
      ideal id = i;                                   //case: id = ideal
      n = ncols(id);
   }
   if( ty=="int" or ty=="intvec" or ty=="intmat" )
   {
      if ( ty == "int" ) { intmat S =id; }
      else { intmat S = intmat(id); }
      intvec i = S[1..nrows(S),1..ncols(S)];
      kill id;
      intvec id = i;                                  //case: id = intvec
      n = size(id);
   }
//--------------- consider intvec v and empty product  -----------------------
   if( tt!=0 )
   {
      for (j=1; j<=size(v); j++)
      {
         if ( v[j] <= 0 or v[j] > n )                 //v outside range of id
         {
            return(1);                                //empty product is 1
         }
      }
      id = id[v];                                     //consider part of id
   }                                                  //corresponding to v
//--------------------- special case: one factor is zero  ---------------------
   if ( typeof(id) == "ideal")
   {
      if( size(id) < ncols(id) )
      {
          poly f; return(f);
      }
   }
//-------------------------- finally, multiply objects  -----------------------
   n = size(id);
   def f(1) = id[1];
   for( j=2; j<=n; j=j+1 ) { def f(j)=f(j-1)*id[j]; }
   return(f(n));
}
example
{  "EXAMPLE:"; echo = 2;
   ring r= 0,(x,y,z),dp;
   ideal m = maxideal(1);
   product(m);
   product(m[2..3]);
   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 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 of positive integers, n=integer,
@*       sort(id[,o,n]); o=string (any allowed ordstr of a ring), n=integer
RETURN:  a list l of two elements:
@format
        l[1]: object of same type as input but sorted in the following way:
           - 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):
             NOTE: generators with SMALLER(!) leading term come FIRST
             (e.g. sort(id); sorts backwards to 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 polynomial in
             Singular, id should not contain negative integers: the result
             might not be as expected
           - 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
         l[2]: intvec, describing the permutation of the input (hence l[2]=v
             if v is given (with positive integers))
@end format
NOTE:    If v is given id may be any simply indexed object (e.g. any list or
         string); if v[i]<0 and i<=size(id) v[i] is set internally to i;
         entries of v must be pairwise distinct to get a permutation 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"
                        or typeof(id)=="matrix"))
   {
      id = simplify(id,2);
      for ( ii=1; ii<ncols(id); ii++ )
      {
         if ( id[ii]!=id[ii+1] ) { break;}
      }
      if ( ii != ncols(id) ) { v = sortvec(id); }
      else  { v = ncols(id)..1; }
   }
   if ( size(#)>=1 and (typeof(id)=="ideal" or typeof(id)=="module"
                        or typeof(id)=="matrix") )
   {
      if ( typeof(#[1])=="string" )
      {
         if( typeof(id)=="matrix" )
         {
           module ID = lead(module(id));
         } else
         {
           def ID = lead(id);
         }
         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")
               { ERROR("// list elements must be intvec/int"); }
            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);
   if( size(id) < s ) { s = size(id); }
   def m = id;
   if ( size(m) != 0 )
   {
      for ( jj=1; jj<=s; jj++)
      {
         if ( v[jj]<=0 ) { v[jj]=jj; }
         m[jj] = id[v[jj]];
      }
   }
   if ( v == 0 ) { v = 1; }
   list L=m,v;
   return(L);
}
example
{  "EXAMPLE:"; echo = 2;
   ring r0 = 0,(x,y,z,t),lp;
   ideal i = x3,z3,xyz;
   sort(i);            //sorts using lex ordering, smaller polys come first

   sort(i,3..1);

   sort(i,"ls")[1];     //sort w.r.t. negative lex ordering

   intvec v =1,10..5,2..4;v;
   sort(v)[1];          // sort v lexicographically

   sort(v,"Dp",1)[1];   // sort v w.r.t (total sum, reverse lex)
}
///////////////////////////////////////////////////////////////////////////////

static proc lsum (int n, list l)
{ if (n>10)
  { return( lsum(n/2,list(l[1..(n/2)])) + lsum(n-n/2, list(l[(n/2+1)..n])) );
  }
  else
  { def Summe=l[1];
    for (int i=2;i<=n;i++)
    { Summe=Summe+l[i];
    }
    return(Summe);
  }
}

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

proc sum (id, list #)
"USAGE:    sum(id[,v]); id ideal/vector/module/matrix/intvec/intmat/list,
          v intvec  (default: v=1..number of entries of id)
ASSUME:   list members can be added.
RETURN:   The sum of all entries of id [with index given by v] of type
          depending on the entries of id.
NOTE:     If id is not a list, id is treated as a list of polys resp. integers.
          A module m is identified with the corresponding matrix M (columns
          of M generate m).
@*        If v is outside the range of id, we have the empty sum and the
          result will be 0 (of type int).
EXAMPLE:  example sum; shows an example
"
{
//-------------------- initialization and special feature ---------------------
   int n,j,tt;
   string ty;                                  // will become type of id
   list l;

// We wish to allow something like sum(x(1..10)) if x(1),...,x(10) are
// variables. x(1..10) is a list of polys and enters the procedure with
// id=x(1) and # a list with 9 polys, #[1]= x(2),...,#[9]= x(10). Hence, in
// this case # is never empty. If an additional intvec v is given,
// it is added to #, so we have to separate it first and make
// the rest a list which has to be added.

   int s = size(#);
   if( s!=0 )
   {  if ( typeof(#[s])=="intvec" or typeof(#[s])=="int")
      {  intvec v = #[s];
         tt=1;
         s=s-1;
         if ( s>0 ) { # = #[1..s]; }
      }
   }
   if ( s>0 )
   {
      l = list(id)+#;
      kill id;
      list id = l;                                    //case: id = list
      ty = "list";
   }
   else
   {
      ty = typeof(id);
   }
//------------------------------ reduce to 3 cases ---------------------------
   if( ty=="poly" or ty=="ideal" or ty=="vector"
       or ty=="module" or ty=="matrix" )
   {                                                 //case: id = ideal
      ideal i = ideal(matrix(id));
      kill id;
      ideal id = simplify(i,2);                      //delete 0 entries
   }
   if( ty=="int" or ty=="intvec" or ty=="intmat" )
   {                                                 //case: id = intvec
      if ( ty == "int" ) { intmat S =id; }
      else { intmat S = intmat(id); }
      intvec i = S[1..nrows(S),1..ncols(S)];
      kill id;
      intvec id = i;
   }
//------------------- consider intvec v and empty sum  -----------------------
   if( tt!=0 )
   {
      for (j=1; j<=size(v); j++)
      {
         if ( v[j] <= 0 or v[j] > size(id) )         //v outside range of id
         {
            return(0);                               //empty sum is 0
         }
      }
      id = id[v];                                    //consider part of id
   }                                                 //corresponding to v

//-------------------------- finally, add objects  ---------------------------
   n = size(id);
   if (n>10)
   { return( lsum(n/2,list(id[1..(n/2)])) + lsum(n-n/2, list(id[(n/2+1)..n])) );
   }
   else
   { def Summe=id[1];
     for (int lauf=2;lauf<=n;lauf++)
     { Summe=Summe+id[lauf];
     }
     return(Summe);
   }
}
example
{  "EXAMPLE:"; echo = 2;
   ring r1 = 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;
   intvec w=2,4,6;
   sum(M,w);
   intvec iv = 1,2,3,4,5,6,7,8,9;
   sum(iv,2..4);
   iv = intvec(1..100);
   sum(iv);
   ring r2 = 0,(x(1..10)),dp;
   sum(x(3..7),intvec(1,3,5));
}
///////////////////////////////////////////////////////////////////////////////


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

proc which (command)
"USAGE:    which(command); command = string expression
RETURN:   Absolute pathname of command, if found in search path.
          Empty string, otherwise.
NOTE:     Based on the Unix command 'which'.
EXAMPLE:  example which; shows an example
"
{
   int rs;
   int i;
   string fn = "which_" + string(system("pid"));
   string pn;
   string cmd;
   if( typeof(command) != "string")
   {
     return (pn);
   }
   if (system("uname") != "ix86-Win")
   {
     cmd = "which ";
   }
   else
   {
     // unfortunately, it does not take -path
     cmd = "type ";
   }
   i = system("sh", cmd + command + " > " + fn);
   pn = read(fn);
   if (system("uname") != "ix86-Win")
   {
     // TBC: Hmm... should parse output to get rid of 'command is '
     pn[size(pn)] = "";
     i = 1;
     while ((pn[i] != " ") and (pn[i] != ""))
     {
       i = i+1;
     }
     if (pn[i] == " ") {pn[i] = "";}
     rs = system("sh", "ls " + pn + " > " + fn + " 2>&1 ");
   }
   else
   {
     rs = 0;
   }
   i = system("sh", "rm " + fn);
   if (rs == 0) {return (pn);}
   else
   {
     print (command + " not found ");
     return ("");
   }
}
example
{  "EXAMPLE:"; echo = 2;
    which("sh");
}
///////////////////////////////////////////////////////////////////////////////

proc watchdog(int i, string cmd)
"USAGE:   watchdog(i,cmd); i integer, cmd string
RETURN:  Result of cmd, if the result can be computed in i seconds.
         Otherwise the computation is interrupted after i seconds,
         the string "Killed" is returned and the global variable
         'watchdog_interrupt' is defined.
NOTE:    * the MP package must be enabled
         * the current basering should not be watchdog_rneu, since
           watchdog_rneu will be killed
         * if there are variable names of the structure x(i) all
           polynomials have to be put into eval(...) in order to be
           interpreted correctly
         * a second Singular process is started by this procedure
EXAMPLE: example watchdog; shows an example
"
{
  string rname=nameof(basering);
  def rsave=basering;
  if (defined(watchdog_rneu))
  {
    kill watchdog_rneu;
  }
// If we do not have MP-links, watchdog cannot be used
  if (system("with","MP"))
  {
    if ( i > 0 )
    {
      int j=10;
      int k=999999;
// fork, get the pid of the child and send it the command
      link l_fork="MPtcp:fork";
      open(l_fork);
      write(l_fork,quote(system("pid")));
      int pid=read(l_fork);
      execute("write(l_fork,quote(" + cmd + "));");


// sleep in small, but growing intervals for appr. 1 second
      while(j < k)
      {
        if (status(l_fork, "read", "ready", j)) {break;}
        j = j + j;
      }

// sleep in intervals of one second
      j = 1;
      if (!status(l_fork,"read","ready"))
      {
        while (j < i)
        {
          if (status(l_fork, "read", "ready", k)) {break;}
          j = j + 1;
        }
      }
// check, whether we have a result, and return it
      if (status(l_fork, "read", "ready"))
      {
        def result = read(l_fork);
        if (nameof(basering)!=rname)
        {
          def watchdog_rneu=basering;
          setring rsave;
          if (!defined(result))
          {
            def result=fetch(watchdog_rneu,result);
          }
        }
        if(defined(watchdog_interrupt))
        {
          kill watchdog_interrupt;
        }
        close(l_fork);
      }
      else
      {
        string result="Killed";
        if(!defined(watchdog_interrupt))
        {
          int watchdog_interrupt=1;
          export watchdog_interrupt;
        }
        close(l_fork);
        j = system("sh","kill " + string(pid));
      }
      return(result);
    }
    else
    {
      ERROR("First argument of watchdog has to be a positive integer.");
    }
  }
  else
  {
    ERROR("MP-support is not enabled in this version of Singular.");
  }
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(x,y,z),dp;
  poly f=x^30+y^30;
  watchdog(1,"factorize(eval("+string(f)+"))");
  watchdog(100,"factorize(eval("+string(f)+"))");
}
///////////////////////////////////////////////////////////////////////////////

proc deleteSublist(intvec v,list l)
"USAGE:   deleteSublist(v,l); intvec v; list l
         where the entries of the integer vector v correspond to the
         positions of the elements to be deleted
RETURN:  list without the deleted elements
EXAMPLE: example deleteSublist; shows an example"
{
  list k;
  int i,j,skip;
  j=1;
  skip=0;
  intvec vs=sort(v)[1];
  for ( i=1 ; i <=size(vs) ; i++)
  {
    while ((j+skip) < vs[i])
    {
      k[j] = l[j+skip];
      j++;
    }
    skip++;
  }
  if(vs[size(vs)]<size(l))
  {
    k=k+list(l[(vs[size(vs)]+1)..size(l)]);
  }
  return(k);
}
example
{ "EXAMPLE:"; echo=2;
   list l=1,2,3,4,5;
   intvec v=1,3,4;
   l=deleteSublist(v,l);
   l;
}
///////////////////////////////////////////////////////////////////////////////
proc primefactors (n, list #)
"USAGE:   primefactors(n [,p]); n = int or number, p = integer
COMPUTE: primefactors <= min(p,32003) of n (default p = 32003)
RETURN:  a list, say l,
         l[1] : primefactors <= min(p,32003) of n
         l[2] : l[2][i] = multiplicity of l[1][i]
         l[3] : remaining factor ( n=product{ (l[1][i]^l[2][i])*l[3]} )
                type(l[3])=typeof(n)
NOTE:    If n is a long integer (of type number) then the procedure
         finds primefactors <= min(p,32003) but n may be as larger as
         2147483647 (max. integer representation)
WARNING: the procedure works for small integers only, just by testing all
         primes (not to be considerd as serious prime factorization!)
EXAMPLE: example primefactors; shows an example
"
{
   int ii,jj,z,p,num,w3,q;
   intvec w1,w2,v;
   list l;
   if (size(#) == 0)
   {
      p=32003;
   }
   else
   {
     if( typeof(#[1]) != "int")
     {
       ERROR("2nd parameter must be of type int"+newline);
     }
     p=#[1];
   }
   if( n<0) { n=-n;};

// ----------------- case: 1st parameter is a number --------------------
   if (typeof(n) =="number")
   {
     kill w3;
     number w3;
     if( n > 2147483647 )           //2147483647 max. integer representation
     {
        v = primes(2,p);
        number m;
        for( ii=1; ii<=size(v); ii++)
        {
          jj=0;
          while(1)
          {
            q  = v[ii];
            jj = jj+1;
            m  = n/q;                  //divide n as often as possible
            if (denominator(m)!=1) { break;  }
            n=m;
          }
          if( jj>1 )
          {
             w1 = w1,v[ii];          //primes
             w2 = w2,jj-1;           //powers
          }
          if( n <= 2147483647 ) { break;  }
        }
     }

     if( n >  2147483647 )            //n is still too big
     {
       if( size(w1) >1 )               //at least 1 primefactor was found
       {
         w1 = w1[2..size(w1)];
         w2 = w2[2..size(w2)];
       }
       else                           //no primefactor was found
       {
         w1 = 1; w2 = 1;
       }
       l  = w1,w2,n;
       return(l);
     }

     if( n <= 2147483647 )          //n is in inter range
     {
       num = int(n);
       kill n;
       int n = num;
     }
   }

// --------------------------- trivial cases --------------------
   if( n==0 )
   {
     w1=1; w2=1; w3=0; l=w1,w2,w3;
     return(l);
   }

   if( n==1 )
   {
       w3=1;
       if( size(w1) >1 )               //at least 1 primefactor was found
       {
         w1 = w1[2..size(w1)];
         w2 = w2[2..size(w2)];
       }
       else                           //no primefactor was found
       {
         w1 = 1; w2 = 1;
       }
      l=w1,w2,w3;
      return(l);
   }
   if ( prime(n)==n )         //note: prime(n) <= 32003 in Singular
   {                          //case n is a prime
      if (p > n)
      {
        w1=w1,n; w2=w2,1; w3=1;
        w1 = w1[2..size(w1)];
        w2 = w2[2..size(w2)];
        l=w1,w2,w3;
        return(l);
      }
      else
      {
        w3=n;
        if( size(w1) >1 )               //at least 1 primefactor was found
        {
           w1 = w1[2..size(w1)];
           w2 = w2[2..size(w2)];
         }
         else                           //no primefactor was found
         {
           w1 = 1; w2 = 1;
         }
         l=w1,w2,w3;
        return(l);
      }
   }
   else
   {
      if ( p >= n)
      {
         v = primes(q,n div 2 + 1);
      }
      else
      {
         v = primes(q,p);
      }
//------------- search for primfactors <= last entry of v ------------
      for(ii=1; ii<=size(v); ii++)
      {
         z=0;
         while( (n mod v[ii]) == 0 )
         {
            z=z+1;
            n = n div v[ii];
         }
         if (z!=0)
         {
            w1 = w1,v[ii];        //primes
            w2 = w2,z;            //multiplicities
         }
      }
   }
//--------------- case:at least 1 primefactor was found ---------------
   if( size(w1) >1 )               //at least 1 primefactor was found
   {
      w1 = w1[2..size(w1)];
      w2 = w2[2..size(w2)];
   }
   else                           //no primefactor was found
   {
     w1 = 1; w2 = 1;
   }
   w3 = n;
   l  = w1,w2,w3;
   return(l);
}
example
{ "EXAMPLE:"; echo = 2;
   primefactors(7*8*121);
   ring r = 0,x,dp;
   primefactors(123456789100);
}

///////////////////////////////////////////////////////////////////////////////
proc primecoeffs(J, list #)
"USAGE:   primecoeffs(J[,p]); J any type which can be converted to a matrix
         e.g. ideal, matrix, vector, module, int, intvec
         p = integer
COMPUTE: primefactors <= min(p,32003) of coeffs of J (default p = 32003)
RETURN:  a list, say l, of two intvectors:@*
         l[1] : the different primefactors of all coefficients of J@*
         l[2] : the different remaining factors
NOTE:    the procedure works for small integers only, just by testing all
         primes (not to be considered as serious prime factorization!)
EXAMPLE: example primecoeffs; shows an example
"
{
   int q,ii,n,mark;;
   if (size(#) == 0)
   {
      q=32003;
   }
   else
   {
     if( typeof(#[1]) != "int")
     {
       ERROR("2nd parameter must be of type int"+newline);
     }
     q=#[1];
   }

   if (defined(basering) == 0)
   {
     mark=1;
     ring r = 0,x,dp;
   }
   def I = ideal(matrix(J));
   poly p = product(maxideal(1));
   matrix Coef=coef(I[1],p);
   ideal id, jd, rest;
   intvec v,re;
   list result,l;
   for(ii=2; ii<=ncols(I); ii++)
   {
     Coef=concat(Coef,coef(I[ii],p));
   }
   id = Coef[2,1..ncols(Coef)];
   id = simplify(id,6);
   for (ii=1; ii<=size(id); ii++)
   {
     l = primefactors(number(id[ii]),q);
     jd = jd,l[1];
     rest = rest,l[3];
   }
   jd = simplify(jd,6);
   for (ii=1; ii<=size(jd); ii++)
   {
     v[ii]=int(jd[ii]);
   }
   v = sort(v)[1];
   rest = simplify(rest,6);
   id = sort(id)[1];
   if (mark)
   {
     for (ii=1; ii<=size(rest); ii++)
     {
       re[ii] = int(rest[ii]);
     }
     result = v,re;
   }
   else
   {
      result = v,rest;
   }
   return(result);
}
example
{ "EXAMPLE:"; echo = 2;
   primecoeffs(intvec(7*8*121,7*8));"";
   ring r = 0,(b,c,t),dp;
   ideal I = -13b6c3t+4b5c4t,-10b4c2t-5b4ct2;
   primecoeffs(I);
}
///////////////////////////////////////////////////////////////////////////////
proc timeFactorize(poly i,list #)
"USAGE:  timeFactorize(p,d); poly p , integer d
RETURN:  factorize(p) if the factorization finished after d-1
         seconds otherwhise f is considered to be irreducible
EXAMPLE: example timeFactorize; shows an example
"
{
  def P=basering;
  if (size(#) > 0)
  {
    if (system("with", "MP"))
    {
      if ((typeof(#[1]) == "int")&&(#[1]))
      {
        int wait = #[1];
        int j = 10;

        string bs = nameof(basering);
        link l_fork = "MPtcp:fork";
        open(l_fork);
        write(l_fork, quote(system("pid")));
        int pid = read(l_fork);
        write(l_fork, quote(timeFactorize(eval(i))));

        // sleep in small intervalls for appr. one second
        if (wait > 0)
        {
          while(j < 1000000)
          {
            if (status(l_fork, "read", "ready", j)) {break;}
            j = j + j;
          }
        }

        // sleep in intervalls of one second from now on
        j = 1;
        while (j < wait)
        {
          if (status(l_fork, "read", "ready", 1000000)) {break;}
          j = j + 1;
        }

        if (status(l_fork, "read", "ready"))
        {
          def result = read(l_fork);
          if (bs != nameof(basering))
          {
            def PP = basering;
            setring P;
            def result = imap(PP, result);
            kill PP;
          }
          kill l_fork;
        }
        else
        {
          list result;
          intvec v=1,1;
          result[1]=list(1,i);
          result[2]=v;
          j = system("sh", "kill " + string(pid));
        }
        return (result);
      }
    }
  }
  return(factorH(i));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y),dp;
   poly p=((x2+y3)^2+xy6)*((x3+y2)^2+x10y);
   p=p^2;
   //timeFactorize(p,2);
   //timeFactorize(p,20);
}

proc timeStd(ideal i,list #)
"USAGE:  timeStd(i,d), i ideal, d integer
RETURN:  std(i) if the standard basis computation finished after
         d-1 seconds and i otherwhise
EXAMPLE: example timeStd; shows an example
"
{
  def P=basering;
  if (size(#) > 0)
  {
    if (system("with", "MP"))
    {
      if ((typeof(#[1]) == "int")&&(#[1]))
      {
        int wait = #[1];
        int j = 10;

        string bs = nameof(basering);
        link l_fork = "MPtcp:fork";
        open(l_fork);
        write(l_fork, quote(system("pid")));
        int pid = read(l_fork);
        write(l_fork, quote(timeStd(eval(i))));

        // sleep in small intervalls for appr. one second
        if (wait > 0)
        {
          while(j < 1000000)
          {
            if (status(l_fork, "read", "ready", j)) {break;}
            j = j + j;
          }
        }
        j = 1;
        while (j < wait)
        {
          if (status(l_fork, "read", "ready", 1000000)) {break;}
          j = j + 1;
        }
        if (status(l_fork, "read", "ready"))
        {
          def result = read(l_fork);
          if (bs != nameof(basering))
          {
            def PP = basering;
            setring P;
            def result = imap(PP, result);
            kill PP;
          }
          kill l_fork;
        }
        else
        {
          ideal result=i;
          j = system("sh", "kill " + string(pid));
        }
        return (result);
      }
    }
  }
  return(std(i));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(a,b,c,d,e),dp;
   int n=6;
   ideal i=
   a^n-b^n,
   b^n-c^n,
   c^n-d^n,
   d^n-e^n,
   a^(n-1)*b+b^(n-1)*c+c^(n-1)*d+d^(n-1)*e+e^(n-1)*a;
   timeStd(i,2);
   timeStd(i,20);
}

proc factorH(poly p)
"USAGE:  factorH(p)  p poly
RETURN:  factorize(p)
NOTE:    changes variables to make the last variable the principal
         one in the multivariate factorization and factorizes then
         the polynomial
EXAMPLE: example factorH; shows an example
"
{
   def R=basering;
   int i,j;
   int n=1;
   int d=nrows(coeffs(p,var(1)));
   for(i=1;i<=nvars(R);i++)
   {
      j=nrows(coeffs(p,var(i)));
      if(d>j)
      {
         n=i;
         d=j;
      }
   }
   ideal ma=maxideal(1);    //die letzte Variable ist die Hauptvariable
   ma[nvars(R)]=var(n);
   ma[n]=var(nvars(R));
   map phi=R,ma;
   list fac=factorize(phi(p));
   list re=phi(fac);
   return(re);
}
example
{ "EXAMPLE:"; echo = 2;
  system("random",992851144);
  ring r=32003,(x,y,z,w,t),lp;
  poly p=y2w9+yz7t-yz5w4-z2w4t4-w8t3;
  factorize(p);  //fast
  system("random",992851262);
  //factorize(p);  //slow
  system("random",992851262);
  factorH(p);
}