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

spielwiese
Last change on this file since 341696 was 341696, checked in by Hans Schönemann <hannes@…>, 14 years ago
Adding Id property to all files git-svn-id: file:///usr/local/Singular/svn/trunk@12231 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 45.8 KB
Line 
1//GMG, last modified 18.6.99
2//anne, added deleteSublist and watchdog 12.12.2000
3//eric, added absValue 11.04.2002
4///////////////////////////////////////////////////////////////////////////////
5version="$Id$";
6category="General purpose";
7info="
8LIBRARY:  general.lib   Elementary Computations of General Type
9
10PROCEDURES:
11 A_Z(\"a\",n);          string a,b,... of n comma separated letters
12 ASCII([n,m]);          string of printable ASCII characters (number n to m)
13 absValue(c);           absolute value of c
14 binomial(n,m[,../..]); n choose m (type int), [type bigint]
15 deleteSublist(iv,l);   delete entries given by iv from list l
16 factorial(n[,../..]);  n factorial (=n!) (type int), [type bigint]
17 fibonacci(n[,p]);      nth Fibonacci number [char p]
18 kmemory([n[,v]]);      active [allocated] memory in kilobyte
19 killall();             kill all user-defined variables
20 number_e(n);           compute exp(1) up to n decimal digits
21 number_pi(n);          compute pi (area of unit circle) up to n digits
22 primes(n,m);           intvec of primes p, n<=p<=m
23 product(../..[,v]);    multiply components of vector/ideal/...[indices v]
24 sort(ideal/module);    sort generators according to monomial ordering
25 sum(vector/id/..[,v]); add components of vector/ideal/...[with indices v]
26 watchdog(i,cmd);       only wait for result of command cmd for i seconds
27 which(command);        search for command and return absolute path, if found
28 primecoeffs(J[,q]);    primefactors <= min(p,32003) of coeffs of J
29 primefactors(n[,p]);  primefactors <= min(p,32003) of n
30 timeStd(i,d)           std(i) if the standard basis computation finished
31                        after d-1 seconds and i otherwhise
32 timeFactorize(p,d)     factorize(p) if the factorization finished after d-1
33                        seconds otherwhise f is considered to be irreducible
34 factorH(p)             changes variables to become the last variable the
35                        principal one in the multivariate factorization and
36                        factorizes then the polynomial
37
38           (parameters in square brackets [] are optional)
39";
40
41LIB "inout.lib";
42LIB "poly.lib";
43LIB "matrix.lib";
44///////////////////////////////////////////////////////////////////////////////
45
46proc A_Z (string s,int n)
47"USAGE:   A_Z(\"a\",n);  a any letter, n integer (-26<= n <=26, !=0)
48RETURN:  string of n small (if a is small) or capital (if a is capital)
49         letters, comma separated, beginning with a, in alphabetical
50         order (or reverse alphabetical order if n<0)
51EXAMPLE: example A_Z; shows an example
52"
53{
54  if ( n>=-26 and n<=26 and n!=0 )
55  {
56    string alpha =
57    "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,"+
58    "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,"+
59    "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,"+
60    "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";
61    int ii; int aa;
62    for(ii=1; ii<=51; ii=ii+2)
63    {
64       if( alpha[ii]==s ) { aa=ii; }
65    }
66    if ( aa==0)
67    {
68      for(ii=105; ii<=155; ii=ii+2)
69      {
70        if( alpha[ii]==s ) { aa=ii; }
71      }
72    }
73    if( aa!=0 )
74    {
75      string out;
76      if (n > 0) { out = alpha[aa,2*(n)-1];  return (out); }
77      if (n < 0)
78      {
79        string beta =
80        "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,"+
81        "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,"+
82        "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,"+
83        "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";
84        if ( aa < 52 ) { aa=52-aa; }
85        if ( aa > 104 ) { aa=260-aa; }
86        out = beta[aa,2*(-n)-1]; return (out);
87      }
88    }
89  }
90}
91example
92{ "EXAMPLE:"; echo = 2;
93   A_Z("c",5);
94   A_Z("Z",-5);
95   string sR = "ring R = (0,"+A_Z("A",6)+"),("+A_Z("a",10)+"),dp;";
96   sR;
97   execute(sR);
98   R;
99}
100///////////////////////////////////////////////////////////////////////////////
101proc ASCII (list #)
102"USAGE:   ASCII([n,m]); n,m= integers (32 <= n <= m <= 126)
103RETURN:   string of printable ASCII characters (no native language support)@*
104          ASCII():    string of  all ASCII characters with its numbers,@*
105          ASCII(n):   n-th ASCII character@*
106          ASCII(n,m): n-th up to m-th ASCII character (inclusive)
107EXAMPLE: example ASCII; shows an example
108"
109{
110  string s1 =
111 "     !    \"    #    $    %    &    '    (    )    *    +    ,    -    .
11232   33   34   35   36   37   38   39   40   41   42   43   44   45   46
113/    0    1    2    3    4    5    6    7    8    9    :    ;    <    =
11447   48   49   50   51   52   53   54   55   56   57   58   59   60   61
115>    ?    @    A    B    C    D    E    F    G    H    I    J    K    L
11662   63   64   65   66   67   68   69   70   71   72   73   74   75   76
117M    N    O    P    Q    R    S    T    U    V    W    X    Y    Z    [
11877   78   79   80   81   82   83   84   85   86   87   88   89   90   91
119\\    ]    ^    _    `    a    b    c    d    e    f    g    h    i    j
12092   93   94   95   96   97   98   99  100  101  102  103  104  105  10
121k    l    m    n    o    p    q    r    s    t    u    v    w    x    y
122107  108  109  110  111  112  113  114  115  116  117  118  119  120  121
123z    {    |    }    ~
124122  123  124  125  126 ";
125
126   string s2 =
127 " !\"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\\]^_`abcdefghijklmnopqrstuvwxyz{|}~";
128
129   if ( size(#) == 0 )
130   {
131      return(s1);
132   }
133   if ( size(#) == 1 )
134   {
135      return( s2[#[1]-31] );
136   }
137   if ( size(#) == 2 )
138   {
139      return( s2[#[1]-31,#[2]-#[1]+1] );
140   }
141}
142example
143{ "EXAMPLE:"; echo = 2;
144   ASCII();"";
145   ASCII(42);
146   ASCII(32,126);
147}
148///////////////////////////////////////////////////////////////////////////////
149
150proc absValue(def c)
151"USAGE:  absValue(c); c int, number or poly
152RETURN:  absValue(c); the absolute value of c
153NOTE:    absValue(c)=c if c>=0; absValue=-c if c<=0.
154@*       So the function can be applied to any type, for which comparison
155@*       operators are defined.
156EXAMPLE: example absValue; shows an example
157"
158{
159  if (c>=0) { return(c); }
160  else { return(-c); }
161}
162example
163{ "EXAMPLE:"; echo = 2;
164   ring r1 = 0,x,dp;
165   absValue(-2002);
166
167   poly f=-4;
168   absValue(f);
169}
170///////////////////////////////////////////////////////////////////////////////
171
172proc binomial (int n, int k)
173"USAGE:   binomial(n,k); n,k integers
174RETURN:  binomial(n,k); binomial coefficient n choose k
175@*       - of type bigint (computed in characteristic 0)
176NOTE:    In any characteristic, binomial(n,k) = coefficient of x^k in (1+x)^n
177SEE ALSO: prime
178EXAMPLE: example binomial; shows an example
179"
180{
181   bigint l;
182   bigint r=1;
183   bigint kk=k;
184   bigint nn=n;
185   if ( k > n-k )
186   { k = n-k; }
187   if ( k<=0 or k>n )               //trivial cases
188   { r = (k==0)*r; }
189   for (l=1; l<=kk; l=l+1 )
190   {
191      r=r*(nn+1-l)/l;
192   }
193   return(r);
194}
195example
196{ "EXAMPLE:"; echo = 2;
197   binomial(200,100);"";                 //type bigint
198   int n,k = 200,100;
199   bigint b1 = binomial(n,k);
200   ring r = 0,x,dp;
201   poly b2 = coeffs((x+1)^n,x)[k+1,1];  //coefficient of x^k in (x+1)^n
202   b1-b2;                               //b1 and b2 should coincide
203}
204///////////////////////////////////////////////////////////////////////////////
205
206static proc binomp (int n, int k, int p)
207 //computes binomial coefficient n choose k in basering of char p > 0
208 //binomial(n,k) = coefficient of x^k in (1+x)^n.
209 //Let n=q*p^j, gcd(q,p)=1, then (1+x)^n = (1 + x^(p^j))^q. We have
210 //binomial(n,k)=0 if k!=l*p^j and binomial(n,l*p^j) = binomial(q,l).
211 //Do this reduction first. Then, in denominator and numerator
212 //of defining formula for binomial coefficient, reduce those factors
213 //mod p which are not divisible by p and cancel common factors p. Hence,
214 //if n = h*p+r, k=l*p+s, r,s<p, binomial(n,k) = binomial(r,s)*binomial(h,l)
215{
216   int l,q,i= 1,n,1;
217   number r=1;
218   if ( k > n-k )
219   { k = n-k;
220   }
221   if ( k<=0 or k>n)               //trivial cases
222   { r = (k==0)*r;
223   }
224   else
225   {
226      while ( q mod p == 0 )
227      {  l = l*p;
228         q = q div p;
229      }                            //we have now n=q*l, l=p^j, gcd(q,p)=1;
230      if (k mod l != 0 )
231      { r = 0;
232      }
233      else
234      {  l = k div l;
235         n = q mod p;
236         k = l mod p;              //now 0<= k,n <p, use binom0 for n,k
237         q = q div p;              //recursion for q,l
238         l = l div p;              //use binomp for q,l
239         r = binom0(n,k)*binomp(q,l,p);
240      }
241   }
242   return(r);
243}
244///////////////////////////////////////////////////////////////////////////////
245
246proc factorial (int n)
247"USAGE:    factorial(n);  n integer
248RETURN:   factorial(n):   n! of type bigint.
249SEE ALSO: prime
250EXAMPLE:  example factorial; shows an example
251"
252{
253   bigint r=1;
254//------------------------------ computation --------------------------------
255   for (int l=2; l<=n; l++)
256   {
257      r=r*l;
258   }
259   return(r);
260}
261example
262{ "EXAMPLE:"; echo = 2;
263   factorial(37);"";                 //37! (as long integer)
264}
265///////////////////////////////////////////////////////////////////////////////
266
267proc fibonacci (int n)
268"USAGE:    fibonacci(n);  n,p integers
269RETURN:   fibonacci(n): nth Fibonacci number, f(0)=f(1)=1, f(i+1)=f(i-1)+f(i)
270@*        - computed in characteristic 0, of type bigint
271SEE ALSO: prime
272EXAMPLE: example fibonacci; shows an example
273"
274{
275  bigint f,g,h=1,1,1;
276//------------------------------ computation --------------------------------
277   for (ii=3; ii<=n; ii=ii+1)
278   {
279      h = f+g; f = g; g = h;
280   }
281   return(h);
282}
283example
284{ "EXAMPLE:"; echo = 2;
285   fibonacci(42); "";             //f(42) of type string (as long integer)
286   ring r = 2,x,dp;
287   number b = fibonacci(42,2);    //f(42) of type number, computed in r
288   b;
289}
290///////////////////////////////////////////////////////////////////////////////
291
292proc kmemory (list #)
293"USAGE:   kmemory([n,[v]]); n,v integers
294RETURN:  memory in kilobyte of type bigint
295@*       n=0: memory used by active variables (same as no parameters)
296@*       n=1: total memory allocated by Singular
297DISPLAY: detailed information about allocated and used memory if v!=0
298NOTE:    kmemory uses internal function 'memory' to compute kilobyte, and
299         is the same as 'memory' for n!=0,1,2
300EXAMPLE: example kmemory; shows an example
301"
302{
303   int n;
304   int verb;
305   if (size(#) != 0)
306   {
307     n=#[1];
308     if (size(#) >1)
309     { verb=#[2]; }
310   }
311
312  if ( verb != 0)
313  {
314    if ( n==0)
315    { dbprint(printlevel-voice+3,
316      "// memory used, at the moment, by active variables (kilobyte):"); }
317    if ( n==1 )
318    { dbprint(printlevel-voice+3,
319      "// total memory allocated, at the moment, by SINGULAR (kilobyte):"); }
320  }
321  return ((memory(n)+1023) div 1024);
322}
323example
324{ "EXAMPLE:"; echo = 2;
325   kmemory();
326   kmemory(1,1);
327}
328///////////////////////////////////////////////////////////////////////////////
329
330proc killall
331"USAGE:   killall(); (no parameter)
332         killall(\"type_name\");
333         killall(\"not\", \"type_name\");
334RETURN:  killall(); kills all user-defined variables except loaded procedures,
335         no return value.
336@*       - killall(\"type_name\"); kills all user-defined variables,
337           of type \"type_name\"
338@*       - killall(\"not\", \"type_name\"); kills all user-defined variables,
339           except those of type \"type_name\" and except loaded procedures
340@*       - killall(\"not\", \"name_1\", \"name_2\", ...);
341           kills all user-defined variables, except those of name \"name_i\"
342           and except loaded procedures
343NOTE:    killall should never be used inside a procedure
344EXAMPLE: example killall; shows an example AND KILLS ALL YOUR VARIABLES
345"
346{
347  list @marie=names(Top);
348  int j, no_kill, @joni;
349  for ( @joni=1; @joni<=size(#);  @joni++)
350  {
351    if (typeof(#[@joni]) != "string")
352    {
353      ERROR("Need string as " + string(@joni) + "th argument");
354    }
355  }
356
357  // kills all user-defined variables but not loaded procedures
358  if( size(#)==0 )
359  {
360    for ( @joni=size(@marie); @joni>0; @joni-- )
361    {
362      if( @marie[@joni]!="LIB" and typeof(`@marie[@joni]`)!="proc"
363      and typeof(`@marie[@joni]`)!="package")
364      { kill `@marie[@joni]`; }
365    }
366  }
367  else
368  {
369    // kills all user-defined variables
370    if( size(#)==1 )
371    {
372      // of type proc
373      if( #[1] == "proc" )
374      {
375        for ( @joni=size(@marie); @joni>0; @joni-- )
376        {
377          if( (@marie[@joni]!="General")
378          and (@marie[@joni]!="Top")
379          and (@marie[@joni]!="killall")
380          and (@marie[@joni]!="LIB") and
381              ((typeof(`@marie[@joni]`)=="package") or
382              (typeof(`@marie[@joni]`)=="proc")))
383          {
384            if (defined(`@marie[@joni]`)) {kill `@marie[@joni]`;}
385          }
386          if (!defined(@joni)) break;
387        }
388        if (defined(General))
389        {
390          @marie=names(General);
391          for ( @joni=size(@marie); @joni>0; @joni-- )
392          {
393            if(@marie[@joni]!="killall"
394            and typeof(`@marie[@joni]`)=="proc")
395            { kill General::`@marie[@joni]`; }
396          }
397          kill General::killall;
398        }
399      }
400      else
401      {
402        // other types
403        for ( @joni=size(@marie); @joni>2; @joni-- )
404        {
405          if(typeof(`@marie[@joni]`)==#[1] and @marie[@joni]!="LIB"
406             and typeof(`@marie[@joni]`)!="proc")
407          { kill `@marie[@joni]`; }
408        }
409      }
410    }
411    else
412    {
413      // kills all user-defined variables whose name or type is not #i
414      for ( @joni=size(@marie); @joni>2; @joni-- )
415      {
416        if ( @marie[@joni] != "LIB" && @marie[@joni] != "Top"
417             && typeof(`@marie[@joni]`) != "proc")
418        {
419          no_kill = 0;
420          for (j=2; j<= size(#); j++)
421          {
422            if (typeof(`@marie[@joni]`)==#[j] or @marie[@joni] == #[j])
423            {
424              no_kill = 1;
425              break;
426            }
427          }
428          if (! no_kill)
429          {
430            kill `@marie[@joni]`;
431          }
432        }
433        if (!defined(@joni)) break;
434      }
435    }
436  }
437}
438example
439{ "EXAMPLE:"; echo = 2;
440   ring rtest; ideal i=x,y,z; string str="hi"; int j = 3;
441   export rtest,i,str,j;       //this makes the local variables global
442   listvar();
443   killall("ring");            // kills all rings
444   listvar();
445   killall("not", "int");      // kills all variables except int's (and procs)
446   listvar();
447   killall();                  // kills all vars except loaded procs
448   listvar();
449}
450///////////////////////////////////////////////////////////////////////////////
451
452proc number_e (int n)
453"USAGE:   number_e(n);  n integer
454RETURN:  Euler number e=exp(1) up to n decimal digits (no rounding)
455@*       - of type string if no basering of char 0 is defined
456@*       - of type number if a basering of char 0 is defined
457DISPLAY: decimal format of e if printlevel > 0 (default:printlevel=0 )
458NOTE:    procedure uses algorithm of A.H.J. Sale
459EXAMPLE: example number_e; shows an example
460"
461{
462   int i,m,s,t;
463   intvec u,e;
464   u[n+2]=0; e[n+1]=0; e=e+1;
465   if( defined(basering) )
466   {
467      if( char(basering)==0 ) { number r=2; t=1; }
468   }
469   string result = "2.";
470   for( i=1; i<=n+1; i=i+1 )
471   {
472      e = e*10;
473      for( m=n+1; m>=1; m=m-1 )
474      {
475         s    = e[m]+u[m+1];
476         u[m] = s div (m+1);
477         e[m] = s%(m+1);
478      }
479      result = result+string(u[1]);
480      if( t==1 ) { r = r+number(u[1])/number(10)^i; }
481   }
482   if( t==1 )
483   { dbprint(printlevel-voice+2,"// "+result[1,n+1]);
484     return(r);
485   }
486   return(result[1,n+1]);
487}
488example
489{ "EXAMPLE:"; echo = 2;
490   number_e(30);"";
491   ring R = 0,t,lp;
492   number e = number_e(30);
493   e;
494}
495///////////////////////////////////////////////////////////////////////////////
496
497proc number_pi (int n)
498"USAGE:   number_pi(n);  n positive integer
499RETURN:  pi (area of unit circle) up to n decimal digits (no rounding)
500@*       - of type string if no basering of char 0 is defined,
501@*       - of type number, if a basering of char 0 is defined
502DISPLAY: decimal format of pi if printlevel > 0 (default:printlevel=0 )
503NOTE:    procedure uses algorithm of S. Rabinowitz
504EXAMPLE: example number_pi; shows an example
505"
506{
507   int i,m,t,e,q,N;
508   intvec r,p,B,Prelim;
509   string result,prelim;
510   N = (10*n) div 3 + 2;
511   p[N+1]=0; p=p+2; r=p;
512   for( i=1; i<=N+1; i=i+1 ) { B[i]=2*i-1; }
513   if( defined(basering) )
514   {
515      if( char(basering)==0 ) { number pi; number pri; t=1; }
516   }
517   for( i=0; i<=n; i=i+1 )
518   {
519      p = r*10;
520      e = p[N+1];
521      for( m=N+1; m>=2; m=m-1 )
522      {
523         r[m] = e%B[m];
524         q    = e div B[m];
525         e    = q*(m-1)+p[m-1];
526      }
527      r[1] = e%10;
528      q    = e div 10;
529      if( q!=10 and q!=9 )
530      {
531         result = result+prelim;
532         Prelim = q;
533         prelim = string(q);
534      }
535      if( q==9 )
536      {
537         Prelim = Prelim,9;
538         prelim = prelim+"9";
539      }
540      if( q==10 )
541      {
542         Prelim = (Prelim+1)-((Prelim+1) div 10)*10;
543         for( m=size(Prelim); m>0; m=m-1)
544         {
545            prelim[m] = string(Prelim[m]);
546         }
547         result = result+prelim;
548         if( t==1 ) { pi=pi+pri; }
549         Prelim = 0;
550         prelim = "0";
551      }
552      if( t==1 ) { pi=pi+number(q)/number(10)^i; }
553   }
554   result = result,prelim[1];
555   result = "3."+result[2,n-1];
556   if( t==1 )
557   { dbprint(printlevel-voice+2,"// "+result);
558     return(pi);
559   }
560   return(result);
561}
562example
563{ "EXAMPLE:"; echo = 2;
564   number_pi(11);"";
565   ring r = (real,10),t,dp;
566   number pi = number_pi(11); pi;
567}
568///////////////////////////////////////////////////////////////////////////////
569
570proc primes (int n, int m)
571"USAGE:   primes(n,m);  n,m integers
572RETURN:  intvec, consisting of all primes p, prime(n)<=p<=m, in increasing
573         order if n<=m, resp. prime(m)<=p<=n, in decreasing order if m<n.
574NOTE:    prime(n); returns the biggest prime number <= min(n,32003)
575         if n>=2, else 2
576EXAMPLE: example primes; shows an example
577"
578{  int change;
579   if ( n>m ) { change=n; n=m ; m=change; change=1; }
580   int q,p = prime(m),prime(n); intvec v = q; q = q-1;
581   while ( q>=p ) { q = prime(q); v = q,v; q = q-1; }
582   if ( change==1 ) { v = v[size(v)..1]; }
583   return(v);
584}
585example
586{  "EXAMPLE:"; echo = 2;
587    primes(50,100);"";
588    intvec v = primes(37,1); v;
589}
590///////////////////////////////////////////////////////////////////////////////
591
592proc product (id, list #)
593"USAGE:    product(id[,v]); id ideal/vector/module/matrix/intvec/intmat/list,
594          v intvec  (default: v=1..number of entries of id)
595ASSUME:   list members can be multiplied.
596RETURN:   The product of all entries of id [with index given by v] of type
597          depending on the entries of id.
598NOTE:     If id is not a list, id is treated as a list of polys resp. integers.
599          A module m is identified with the corresponding matrix M (columns
600          of M generate m).
601@*        If v is outside the range of id, we have the empty product and the
602          result will be 1 (of type int).
603EXAMPLE:  example product; shows an example
604"
605{
606//-------------------- initialization and special feature ---------------------
607   int n,j,tt;
608   string ty;                                //will become type of id
609   list l;
610
611// We wish to allow something like product(x(1..10)) if x(1),...,x(10) are
612// variables. x(1..10) is a list of polys and enters the procedure with
613// id=x(1) and # a list with 9 polys, #[1]= x(2),...,#[9]= x(10). Hence, in
614// this case # is never empty. If an additional intvec v is given,
615// it is added to #, so we have to separate it first and make
616// the rest a list which has to be multiplied.
617
618   int s = size(#);
619   if( s!=0 )
620   {  if ( typeof(#[s])=="intvec" or typeof(#[s])=="int")
621      {
622         intvec v = #[s];
623         tt=1;
624         s=s-1;
625         if ( s>0 ) { # = #[1..s]; }
626      }
627   }
628   if ( s>0 )
629   {
630      l = list(id)+#;
631      kill id;
632      list id = l;                                    //case: id = list
633      ty = "list";
634      n = size(id);
635   }
636   else
637   {
638      ty = typeof(id);
639      if( ty == "list" )
640      { n = size(id); }
641   }
642//------------------------------ reduce to 3 cases ---------------------------
643   if( ty=="poly" or ty=="ideal" or ty=="vector"
644       or ty=="module" or ty=="matrix" )
645   {
646      ideal i = ideal(matrix(id));
647      kill id;
648      ideal id = i;                                   //case: id = ideal
649      n = ncols(id);
650   }
651   if( ty=="int" or ty=="intvec" or ty=="intmat" )
652   {
653      if ( ty == "int" ) { intmat S =id; }
654      else { intmat S = intmat(id); }
655      intvec i = S[1..nrows(S),1..ncols(S)];
656      kill id;
657      intvec id = i;                                  //case: id = intvec
658      n = size(id);
659   }
660//--------------- consider intvec v and empty product  -----------------------
661   if( tt!=0 )
662   {
663      for (j=1; j<=size(v); j++)
664      {
665         if ( v[j] <= 0 or v[j] > n )                 //v outside range of id
666         {
667            return(1);                                //empty product is 1
668         }
669      }
670      id = id[v];                                     //consider part of id
671   }                                                  //corresponding to v
672//--------------------- special case: one factor is zero  ---------------------
673   if ( typeof(id) == "ideal")
674   {
675      if( size(id) < ncols(id) )
676      {
677          poly f; return(f);
678      }
679   }
680//-------------------------- finally, multiply objects  -----------------------
681   n = size(id);
682   def f(1) = id[1];
683   for( j=2; j<=n; j=j+1 ) { def f(j)=f(j-1)*id[j]; }
684   return(f(n));
685}
686example
687{  "EXAMPLE:"; echo = 2;
688   ring r= 0,(x,y,z),dp;
689   ideal m = maxideal(1);
690   product(m);
691   product(m[2..3]);
692   matrix M[2][3] = 1,x,2,y,3,z;
693   product(M);
694   intvec v=2,4,6;
695   product(M,v);
696   intvec iv = 1,2,3,4,5,6,7,8,9;
697   v=1..5,7,9;
698   product(iv,v);
699   intmat A[2][3] = 1,1,1,2,2,2;
700   product(A,3..5);
701}
702///////////////////////////////////////////////////////////////////////////////
703
704proc sort (id, list #)
705"USAGE:   sort(id[,v,o,n]); id = ideal/module/intvec/list(of intvec's or int's)
706@*       sort may be called with 1, 2 or 3 arguments in the following way:
707@*       sort(id[,v,n]); v=intvec of positive integers, n=integer,
708@*       sort(id[,o,n]); o=string (any allowed ordstr of a ring), n=integer
709RETURN:  a list l of two elements:
710@format
711        l[1]: object of same type as input but sorted in the following way:
712           - if id=ideal/module: generators of id are sorted w.r.t. intvec v
713             (id[v[1]] becomes 1-st, id[v[2]]  2-nd element, etc.). If no v is
714             present, id is sorted w.r.t. ordering o (if o is given) or w.r.t.
715             actual monomial ordering (if no o is given):
716             NOTE: generators with SMALLER(!) leading term come FIRST
717             (e.g. sort(id); sorts backwards to actual monomial ordering)
718           - if id=list of intvec's or int's: consider a list element, say
719             id[1]=3,2,5, as exponent vector of the monomial x^3*y^2*z^5;
720             the corresponding monomials are ordered w.r.t. intvec v (s.a.).
721             If no v is present, the monomials are sorted w.r.t. ordering o
722             (if o is given) or w.r.t. lexicographical ordering (if no o is
723             given). The corresponding ordered list of exponent vectors is
724             returned.
725             (e.g. sort(id); sorts lexicographically, smaller int's come first)
726             WARNING: Since negative exponents create the 0 polynomial in
727             Singular, id should not contain negative integers: the result
728             might not be as expected
729           - if id=intvec: id is treated as list of integers
730           - if n!=0 the ordering is inverse, i.e. w.r.t. v(size(v)..1)
731             default: n=0
732         l[2]: intvec, describing the permutation of the input (hence l[2]=v
733             if v is given (with positive integers))
734@end format
735NOTE:    If v is given id may be any simply indexed object (e.g. any list or
736         string); if v[i]<0 and i<=size(id) v[i] is set internally to i;
737         entries of v must be pairwise distinct to get a permutation id.
738         Zero generators of ideal/module are deleted
739         If 'o' is given, the input is sorted by considering leading terms
740         w.r.t. the new ring ordering given by 'o'
741EXAMPLE: example sort; shows an example
742"
743{  int ii,jj,s,n = 0,0,1,0;
744   intvec v;
745   if ( defined(basering) ) { def P = basering; }
746   if ( size(#)==0 and (typeof(id)=="ideal" or typeof(id)=="module"
747                        or typeof(id)=="matrix"))
748   {
749      id = simplify(id,2);
750      for ( ii=1; ii<ncols(id); ii++ )
751      {
752         if ( id[ii]!=id[ii+1] ) { break;}
753      }
754      if ( ii != ncols(id) ) { v = sortvec(id); }
755      else  { v = ncols(id)..1; }
756   }
757   if ( size(#)>=1 and (typeof(id)=="ideal" or typeof(id)=="module"
758                        or typeof(id)=="matrix") )
759   {
760      if ( typeof(#[1])=="string" )
761      {
762         execute("ring r1 =("+charstr(P)+"),("+varstr(P)+"),("+#[1]+");");
763         def i = imap(P,id);
764         v = sortvec(i);
765         setring P;
766         n=2;
767      }
768   }
769   if ( typeof(id)=="intvec" or typeof(id)=="list" and n==0 )
770   {
771      string o;
772      if ( size(#)==0 ) { o = "lp"; n=1; }
773      if ( size(#)>=1 )
774      {
775         if ( typeof(#[1])=="string" ) { o = #[1]; n=1; }
776      }
777   }
778   if ( typeof(id)=="intvec" or typeof(id)=="list" and n==1 )
779   {
780      if ( typeof(id)=="list" )
781      {
782         for (ii=1; ii<=size(id); ii++)
783         {
784            if (typeof(id[ii]) != "intvec" and typeof(id[ii]) != "int")
785               { ERROR("// list elements must be intvec/int"); }
786            else
787               { s=size(id[ii])*(s < size(id[ii])) + s*(s >= size(id[ii])); }
788         }
789      }
790      execute("ring r=0,x(1..s),("+o+");");
791      ideal i;
792      poly f;
793      for (ii=1; ii<=size(id); ii++)
794      {
795         f=1;
796         for (jj=1; jj<=size(id[ii]); jj++)
797         {
798            f=f*x(jj)^(id[ii])[jj];
799         }
800         i[ii]=f;
801      }
802      v = sort(i)[2];
803   }
804   if ( size(#)!=0 and n==0 ) { v = #[1]; }
805   if( size(#)==2 )
806   {
807      if ( #[2] != 0 ) { v = v[size(v)..1]; }
808   }
809   s = size(v);
810   if( size(id) < s ) { s = size(id); }
811   def m = id;
812   if ( size(m) != 0 )
813   {
814      for ( jj=1; jj<=s; jj++)
815      {
816         if ( v[jj]<=0 ) { v[jj]=jj; }
817         m[jj] = id[v[jj]];
818      }
819   }
820   if ( v == 0 ) { v = 1; }
821   list L=m,v;
822   return(L);
823}
824example
825{  "EXAMPLE:"; echo = 2;
826   ring r0 = 0,(x,y,z,t),lp;
827   ideal i = x3,z3,xyz;
828   sort(i);            //sorts using lex ordering, smaller polys come first
829
830   sort(i,3..1);
831
832   sort(i,"ls")[1];     //sort w.r.t. negative lex ordering
833
834   intvec v =1,10..5,2..4;v;
835   sort(v)[1];          // sort v lexicographically
836
837   sort(v,"Dp",1)[1];   // sort v w.r.t (total sum, reverse lex)
838
839   // Note that in general: lead(sort(M)) != sort(lead(M)), e.g:
840   module M = [0, 1, 1, 0], [1, 0, 0, 1]; M;
841   sort(lead(M), "c, dp")[1];
842   lead(sort(M, "c, dp")[1]);
843
844   // In order to sort M wrt a NEW ordering by considering OLD leading
845   // terms use one of the following equivalent commands:
846   module( M[ sort(lead(M), "c,dp")[2] ] );
847   sort( M, sort(lead(M), "c,dp")[2] )[1];
848}
849///////////////////////////////////////////////////////////////////////////////
850
851static proc lsum (int n, list l)
852{ if (n>10)
853  { return( lsum(n/2,list(l[1..(n/2)])) + lsum(n-n/2, list(l[(n/2+1)..n])) );
854  }
855  else
856  { def Summe=l[1];
857    for (int i=2;i<=n;i++)
858    { Summe=Summe+l[i];
859    }
860    return(Summe);
861  }
862}
863
864///////////////////////////////////////////////////////////////////////////////
865
866proc sum (id, list #)
867"USAGE:    sum(id[,v]); id ideal/vector/module/matrix/intvec/intmat/list,
868          v intvec  (default: v=1..number of entries of id)
869ASSUME:   list members can be added.
870RETURN:   The sum of all entries of id [with index given by v] of type
871          depending on the entries of id.
872NOTE:     If id is not a list, id is treated as a list of polys resp. integers.
873          A module m is identified with the corresponding matrix M (columns
874          of M generate m).
875@*        If v is outside the range of id, we have the empty sum and the
876          result will be 0 (of type int).
877EXAMPLE:  example sum; shows an example
878"
879{
880//-------------------- initialization and special feature ---------------------
881   int n,j,tt;
882   string ty;                                  // will become type of id
883   list l;
884
885// We wish to allow something like sum(x(1..10)) if x(1),...,x(10) are
886// variables. x(1..10) is a list of polys and enters the procedure with
887// id=x(1) and # a list with 9 polys, #[1]= x(2),...,#[9]= x(10). Hence, in
888// this case # is never empty. If an additional intvec v is given,
889// it is added to #, so we have to separate it first and make
890// the rest a list which has to be added.
891
892   int s = size(#);
893   if( s!=0 )
894   {  if ( typeof(#[s])=="intvec" or typeof(#[s])=="int")
895      {  intvec v = #[s];
896         tt=1;
897         s=s-1;
898         if ( s>0 ) { # = #[1..s]; }
899      }
900   }
901   if ( s>0 )
902   {
903      l = list(id)+#;
904      kill id;
905      list id = l;                                    //case: id = list
906      ty = "list";
907   }
908   else
909   {
910      ty = typeof(id);
911   }
912//------------------------------ reduce to 3 cases ---------------------------
913   if( ty=="poly" or ty=="ideal" or ty=="vector"
914       or ty=="module" or ty=="matrix" )
915   {                                                 //case: id = ideal
916      ideal i = ideal(matrix(id));
917      kill id;
918      ideal id = simplify(i,2);                      //delete 0 entries
919   }
920   if( ty=="int" or ty=="intvec" or ty=="intmat" )
921   {                                                 //case: id = intvec
922      if ( ty == "int" ) { intmat S =id; }
923      else { intmat S = intmat(id); }
924      intvec i = S[1..nrows(S),1..ncols(S)];
925      kill id;
926      intvec id = i;
927   }
928//------------------- consider intvec v and empty sum  -----------------------
929   if( tt!=0 )
930   {
931      for (j=1; j<=size(v); j++)
932      {
933         if ( v[j] <= 0 or v[j] > size(id) )         //v outside range of id
934         {
935            return(0);                               //empty sum is 0
936         }
937      }
938      id = id[v];                                    //consider part of id
939   }                                                 //corresponding to v
940
941//-------------------------- finally, add objects  ---------------------------
942   n = size(id);
943   if (n>10)
944   { return( lsum(n/2,list(id[1..(n/2)])) + lsum(n-n/2, list(id[(n/2+1)..n])) );
945   }
946   else
947   { def Summe=id[1];
948     for (int lauf=2;lauf<=n;lauf++)
949     { Summe=Summe+id[lauf];
950     }
951     return(Summe);
952   }
953}
954example
955{  "EXAMPLE:"; echo = 2;
956   ring r1 = 0,(x,y,z),dp;
957   vector pv = [xy,xz,yz,x2,y2,z2];
958   sum(pv);
959   sum(pv,2..5);
960   matrix M[2][3] = 1,x,2,y,3,z;
961   intvec w=2,4,6;
962   sum(M,w);
963   intvec iv = 1,2,3,4,5,6,7,8,9;
964   sum(iv,2..4);
965   iv = intvec(1..100);
966   sum(iv);
967   ring r2 = 0,(x(1..10)),dp;
968   sum(x(3..7),intvec(1,3,5));
969}
970///////////////////////////////////////////////////////////////////////////////
971
972
973///////////////////////////////////////////////////////////////////////////////
974
975proc which (command)
976"USAGE:    which(command); command = string expression
977RETURN:   Absolute pathname of command, if found in search path.
978          Empty string, otherwise.
979NOTE:     Based on the Unix command 'which'.
980EXAMPLE:  example which; shows an example
981"
982{
983   int rs;
984   int i;
985   string fn = "which_" + string(system("pid"));
986   string pn;
987   string cmd;
988   if( typeof(command) != "string")
989   {
990     return (pn);
991   }
992   if (system("uname") != "ix86-Win")
993   {
994     cmd = "which ";
995   }
996   else
997   {
998     // unfortunately, it does not take -path
999     cmd = "type ";
1000   }
1001   i = system("sh", cmd + command + " > " + fn);
1002   pn = read(fn);
1003   if (system("uname") != "ix86-Win")
1004   {
1005     // TBC: Hmm... should parse output to get rid of 'command is '
1006     pn[size(pn)] = "";
1007     i = 1;
1008     while ((pn[i] != " ") and (pn[i] != ""))
1009     {
1010       i = i+1;
1011     }
1012     if (pn[i] == " ") {pn[i] = "";}
1013     rs = system("sh", "ls " + pn + " > " + fn + " 2>&1 ");
1014   }
1015   else
1016   {
1017     rs = 0;
1018   }
1019   i = system("sh", "rm " + fn);
1020   if (rs == 0) {return (pn);}
1021   else
1022   {
1023     print (command + " not found ");
1024     return ("");
1025   }
1026}
1027example
1028{  "EXAMPLE:"; echo = 2;
1029    which("sh");
1030}
1031///////////////////////////////////////////////////////////////////////////////
1032
1033proc watchdog(int i, string cmd)
1034"USAGE:   watchdog(i,cmd); i integer, cmd string
1035RETURN:  Result of cmd, if the result can be computed in i seconds.
1036         Otherwise the computation is interrupted after i seconds,
1037         the string "Killed" is returned and the global variable
1038         'watchdog_interrupt' is defined.
1039NOTE:    * the MP package must be enabled
1040         * the current basering should not be watchdog_rneu, since
1041           watchdog_rneu will be killed
1042         * if there are variable names of the structure x(i) all
1043           polynomials have to be put into eval(...) in order to be
1044           interpreted correctly
1045         * a second Singular process is started by this procedure
1046EXAMPLE: example watchdog; shows an example
1047"
1048{
1049  string rname=nameof(basering);
1050  def rsave=basering;
1051  if (defined(watchdog_rneu))
1052  {
1053    kill watchdog_rneu;
1054  }
1055// If we do not have MP-links, watchdog cannot be used
1056  if (system("with","MP"))
1057  {
1058    if ( i > 0 )
1059    {
1060      int j=10;
1061      int k=999999;
1062// fork, get the pid of the child and send it the command
1063      link l_fork="MPtcp:fork";
1064      open(l_fork);
1065      write(l_fork,quote(system("pid")));
1066      int pid=read(l_fork);
1067      execute("write(l_fork,quote(" + cmd + "));");
1068
1069
1070// sleep in small, but growing intervals for appr. 1 second
1071      while(j < k)
1072      {
1073        if (status(l_fork, "read", "ready", j)) {break;}
1074        j = j + j;
1075      }
1076
1077// sleep in intervals of one second
1078      j = 1;
1079      if (!status(l_fork,"read","ready"))
1080      {
1081        while (j < i)
1082        {
1083          if (status(l_fork, "read", "ready", k)) {break;}
1084          j = j + 1;
1085        }
1086      }
1087// check, whether we have a result, and return it
1088      if (status(l_fork, "read", "ready"))
1089      {
1090        def result = read(l_fork);
1091        if (nameof(basering)!=rname)
1092        {
1093          def watchdog_rneu=basering;
1094          setring rsave;
1095          if (!defined(result))
1096          {
1097            def result=fetch(watchdog_rneu,result);
1098          }
1099        }
1100        if(defined(watchdog_interrupt))
1101        {
1102          kill watchdog_interrupt;
1103        }
1104        close(l_fork);
1105      }
1106      else
1107      {
1108        string result="Killed";
1109        if(!defined(watchdog_interrupt))
1110        {
1111          int watchdog_interrupt=1;
1112          export watchdog_interrupt;
1113        }
1114        close(l_fork);
1115        j = system("sh","kill " + string(pid));
1116      }
1117      return(result);
1118    }
1119    else
1120    {
1121      ERROR("First argument of watchdog has to be a positive integer.");
1122    }
1123  }
1124  else
1125  {
1126    ERROR("MP-support is not enabled in this version of Singular.");
1127  }
1128}
1129example
1130{ "EXAMPLE:"; echo=2;
1131  ring r=0,(x,y,z),dp;
1132  poly f=x^30+y^30;
1133  watchdog(1,"factorize(eval("+string(f)+"))");
1134  watchdog(100,"factorize(eval("+string(f)+"))");
1135}
1136///////////////////////////////////////////////////////////////////////////////
1137
1138proc deleteSublist(intvec v,list l)
1139"USAGE:   deleteSublist(v,l); intvec v; list l
1140         where the entries of the integer vector v correspond to the
1141         positions of the elements to be deleted
1142RETURN:  list without the deleted elements
1143EXAMPLE: example deleteSublist; shows an example"
1144{
1145  list k;
1146  int i,j,skip;
1147  j=1;
1148  skip=0;
1149  intvec vs=sort(v)[1];
1150  for ( i=1 ; i <=size(vs) ; i++)
1151  {
1152    while ((j+skip) < vs[i])
1153    {
1154      k[j] = l[j+skip];
1155      j++;
1156    }
1157    skip++;
1158  }
1159  if(vs[size(vs)]<size(l))
1160  {
1161    k=k+list(l[(vs[size(vs)]+1)..size(l)]);
1162  }
1163  return(k);
1164}
1165example
1166{ "EXAMPLE:"; echo=2;
1167   list l=1,2,3,4,5;
1168   intvec v=1,3,4;
1169   l=deleteSublist(v,l);
1170   l;
1171}
1172///////////////////////////////////////////////////////////////////////////////
1173proc primefactors (n, list #)
1174"USAGE:   primefactors(n [,p]); n = int or number, p = integer
1175COMPUTE: primefactors <= min(p,32003) of n (default p = 32003)
1176RETURN:  a list, say l,
1177         l[1] : primefactors <= min(p,32003) of n
1178         l[2] : l[2][i] = multiplicity of l[1][i]
1179         l[3] : remaining factor ( n=product{ (l[1][i]^l[2][i])*l[3]} )
1180                type(l[3])=typeof(n)
1181NOTE:    If n is a long integer (of type number) then the procedure
1182         finds primefactors <= min(p,32003) but n may be as larger as
1183         2147483647 (max. integer representation)
1184WARNING: the procedure works for small integers only, just by testing all
1185         primes (not to be considerd as serious prime factorization!)
1186EXAMPLE: example primefactors; shows an example
1187"
1188{
1189   int ii,jj,z,p,num,w3,q;
1190   intvec w1,w2,v;
1191   list l;
1192   if (size(#) == 0)
1193   {
1194      p=32003;
1195   }
1196   else
1197   {
1198     if( typeof(#[1]) != "int")
1199     {
1200       ERROR("2nd parameter must be of type int"+newline);
1201     }
1202     p=#[1];
1203   }
1204   if( n<0) { n=-n;};
1205
1206// ----------------- case: 1st parameter is a number --------------------
1207   if (typeof(n) =="number")
1208   {
1209     kill w3;
1210     number w3;
1211     if( n > 2147483647 )           //2147483647 max. integer representation
1212     {
1213        v = primes(2,p);
1214        number m;
1215        for( ii=1; ii<=size(v); ii++)
1216        {
1217          jj=0;
1218          while(1)
1219          {
1220            q  = v[ii];
1221            jj = jj+1;
1222            m  = n/q;                  //divide n as often as possible
1223            if (denominator(m)!=1) { break;  }
1224            n=m;
1225          }
1226          if( jj>1 )
1227          {
1228             w1 = w1,v[ii];          //primes
1229             w2 = w2,jj-1;           //powers
1230          }
1231          if( n <= 2147483647 ) { break;  }
1232        }
1233     }
1234
1235     if( n >  2147483647 )            //n is still too big
1236     {
1237       if( size(w1) >1 )               //at least 1 primefactor was found
1238       {
1239         w1 = w1[2..size(w1)];
1240         w2 = w2[2..size(w2)];
1241       }
1242       else                           //no primefactor was found
1243       {
1244         w1 = 1; w2 = 1;
1245       }
1246       l  = w1,w2,n;
1247       return(l);
1248     }
1249
1250     if( n <= 2147483647 )          //n is in inter range
1251     {
1252       num = int(n);
1253       kill n;
1254       int n = num;
1255     }
1256   }
1257
1258// --------------------------- trivial cases --------------------
1259   if( n==0 )
1260   {
1261     w1=1; w2=1; w3=0; l=w1,w2,w3;
1262     return(l);
1263   }
1264
1265   if( n==1 )
1266   {
1267       w3=1;
1268       if( size(w1) >1 )               //at least 1 primefactor was found
1269       {
1270         w1 = w1[2..size(w1)];
1271         w2 = w2[2..size(w2)];
1272       }
1273       else                           //no primefactor was found
1274       {
1275         w1 = 1; w2 = 1;
1276       }
1277      l=w1,w2,w3;
1278      return(l);
1279   }
1280   if ( prime(n)==n )         //note: prime(n) <= 32003 in Singular
1281   {                          //case n is a prime
1282      if (p > n)
1283      {
1284        w1=w1,n; w2=w2,1; w3=1;
1285        w1 = w1[2..size(w1)];
1286        w2 = w2[2..size(w2)];
1287        l=w1,w2,w3;
1288        return(l);
1289      }
1290      else
1291      {
1292        w3=n;
1293        if( size(w1) >1 )               //at least 1 primefactor was found
1294        {
1295           w1 = w1[2..size(w1)];
1296           w2 = w2[2..size(w2)];
1297         }
1298         else                           //no primefactor was found
1299         {
1300           w1 = 1; w2 = 1;
1301         }
1302         l=w1,w2,w3;
1303        return(l);
1304      }
1305   }
1306   else
1307   {
1308      if ( p >= n)
1309      {
1310         v = primes(q,n div 2 + 1);
1311      }
1312      else
1313      {
1314         v = primes(q,p);
1315      }
1316//------------- search for primfactors <= last entry of v ------------
1317      for(ii=1; ii<=size(v); ii++)
1318      {
1319         z=0;
1320         while( (n mod v[ii]) == 0 )
1321         {
1322            z=z+1;
1323            n = n div v[ii];
1324         }
1325         if (z!=0)
1326         {
1327            w1 = w1,v[ii];        //primes
1328            w2 = w2,z;            //multiplicities
1329         }
1330      }
1331   }
1332//--------------- case:at least 1 primefactor was found ---------------
1333   if( size(w1) >1 )               //at least 1 primefactor was found
1334   {
1335      w1 = w1[2..size(w1)];
1336      w2 = w2[2..size(w2)];
1337   }
1338   else                           //no primefactor was found
1339   {
1340     w1 = 1; w2 = 1;
1341   }
1342   w3 = n;
1343   l  = w1,w2,w3;
1344   return(l);
1345}
1346example
1347{ "EXAMPLE:"; echo = 2;
1348   primefactors(7*8*121);
1349   ring r = 0,x,dp;
1350   primefactors(123456789100);
1351}
1352
1353///////////////////////////////////////////////////////////////////////////////
1354proc primecoeffs(J, list #)
1355"USAGE:   primecoeffs(J[,p]); J any type which can be converted to a matrix
1356         e.g. ideal, matrix, vector, module, int, intvec
1357         p = integer
1358COMPUTE: primefactors <= min(p,32003) of coeffs of J (default p = 32003)
1359RETURN:  a list, say l, of two intvectors:@*
1360         l[1] : the different primefactors of all coefficients of J@*
1361         l[2] : the different remaining factors
1362NOTE:    the procedure works for small integers only, just by testing all
1363         primes (not to be considered as serious prime factorization!)
1364EXAMPLE: example primecoeffs; shows an example
1365"
1366{
1367   int q,ii,n,mark;;
1368   if (size(#) == 0)
1369   {
1370      q=32003;
1371   }
1372   else
1373   {
1374     if( typeof(#[1]) != "int")
1375     {
1376       ERROR("2nd parameter must be of type int"+newline);
1377     }
1378     q=#[1];
1379   }
1380
1381   if (defined(basering) == 0)
1382   {
1383     mark=1;
1384     ring r = 0,x,dp;
1385   }
1386   def I = ideal(matrix(J));
1387   poly p = product(maxideal(1));
1388   matrix Coef=coef(I[1],p);
1389   ideal id, jd, rest;
1390   intvec v,re;
1391   list result,l;
1392   for(ii=2; ii<=ncols(I); ii++)
1393   {
1394     Coef=concat(Coef,coef(I[ii],p));
1395   }
1396   id = Coef[2,1..ncols(Coef)];
1397   id = simplify(id,6);
1398   for (ii=1; ii<=size(id); ii++)
1399   {
1400     l = primefactors(number(id[ii]),q);
1401     jd = jd,l[1];
1402     rest = rest,l[3];
1403   }
1404   jd = simplify(jd,6);
1405   for (ii=1; ii<=size(jd); ii++)
1406   {
1407     v[ii]=int(jd[ii]);
1408   }
1409   v = sort(v)[1];
1410   rest = simplify(rest,6);
1411   id = sort(id)[1];
1412   if (mark)
1413   {
1414     for (ii=1; ii<=size(rest); ii++)
1415     {
1416       re[ii] = int(rest[ii]);
1417     }
1418     result = v,re;
1419   }
1420   else
1421   {
1422      result = v,rest;
1423   }
1424   return(result);
1425}
1426example
1427{ "EXAMPLE:"; echo = 2;
1428   primecoeffs(intvec(7*8*121,7*8));"";
1429   ring r = 0,(b,c,t),dp;
1430   ideal I = -13b6c3t+4b5c4t,-10b4c2t-5b4ct2;
1431   primecoeffs(I);
1432}
1433///////////////////////////////////////////////////////////////////////////////
1434proc timeFactorize(poly i,list #)
1435"USAGE:  timeFactorize(p,d); poly p , integer d
1436RETURN:  factorize(p) if the factorization finished after d-1
1437         seconds otherwhise f is considered to be irreducible
1438EXAMPLE: example timeFactorize; shows an example
1439"
1440{
1441  def P=basering;
1442  if (size(#) > 0)
1443  {
1444    if (system("with", "MP"))
1445    {
1446      if ((typeof(#[1]) == "int")&&(#[1]))
1447      {
1448        int wait = #[1];
1449        int j = 10;
1450
1451        string bs = nameof(basering);
1452        link l_fork = "MPtcp:fork";
1453        open(l_fork);
1454        write(l_fork, quote(system("pid")));
1455        int pid = read(l_fork);
1456        write(l_fork, quote(timeFactorize(eval(i))));
1457
1458        // sleep in small intervalls for appr. one second
1459        if (wait > 0)
1460        {
1461          while(j < 1000000)
1462          {
1463            if (status(l_fork, "read", "ready", j)) {break;}
1464            j = j + j;
1465          }
1466        }
1467
1468        // sleep in intervalls of one second from now on
1469        j = 1;
1470        while (j < wait)
1471        {
1472          if (status(l_fork, "read", "ready", 1000000)) {break;}
1473          j = j + 1;
1474        }
1475
1476        if (status(l_fork, "read", "ready"))
1477        {
1478          def result = read(l_fork);
1479          if (bs != nameof(basering))
1480          {
1481            def PP = basering;
1482            setring P;
1483            def result = imap(PP, result);
1484            kill PP;
1485          }
1486          kill l_fork;
1487        }
1488        else
1489        {
1490          list result;
1491          intvec v=1,1;
1492          result[1]=list(1,i);
1493          result[2]=v;
1494          j = system("sh", "kill " + string(pid));
1495        }
1496        return (result);
1497      }
1498    }
1499  }
1500  return(factorH(i));
1501}
1502example
1503{ "EXAMPLE:"; echo = 2;
1504   ring r=0,(x,y),dp;
1505   poly p=((x2+y3)^2+xy6)*((x3+y2)^2+x10y);
1506   p=p^2;
1507   //timeFactorize(p,2);
1508   //timeFactorize(p,20);
1509}
1510
1511proc timeStd(ideal i,list #)
1512"USAGE:  timeStd(i,d), i ideal, d integer
1513RETURN:  std(i) if the standard basis computation finished after
1514         d-1 seconds and i otherwhise
1515EXAMPLE: example timeStd; shows an example
1516"
1517{
1518  def P=basering;
1519  if (size(#) > 0)
1520  {
1521    if (system("with", "MP"))
1522    {
1523      if ((typeof(#[1]) == "int")&&(#[1]))
1524      {
1525        int wait = #[1];
1526        int j = 10;
1527
1528        string bs = nameof(basering);
1529        link l_fork = "MPtcp:fork";
1530        open(l_fork);
1531        write(l_fork, quote(system("pid")));
1532        int pid = read(l_fork);
1533        write(l_fork, quote(timeStd(eval(i))));
1534
1535        // sleep in small intervalls for appr. one second
1536        if (wait > 0)
1537        {
1538          while(j < 1000000)
1539          {
1540            if (status(l_fork, "read", "ready", j)) {break;}
1541            j = j + j;
1542          }
1543        }
1544        j = 1;
1545        while (j < wait)
1546        {
1547          if (status(l_fork, "read", "ready", 1000000)) {break;}
1548          j = j + 1;
1549        }
1550        if (status(l_fork, "read", "ready"))
1551        {
1552          def result = read(l_fork);
1553          if (bs != nameof(basering))
1554          {
1555            def PP = basering;
1556            setring P;
1557            def result = imap(PP, result);
1558            kill PP;
1559          }
1560          kill l_fork;
1561        }
1562        else
1563        {
1564          ideal result=i;
1565          j = system("sh", "kill " + string(pid));
1566        }
1567        return (result);
1568      }
1569    }
1570  }
1571  return(std(i));
1572}
1573example
1574{ "EXAMPLE:"; echo = 2;
1575   ring r=32003,(a,b,c,d,e),dp;
1576   int n=6;
1577   ideal i=
1578   a^n-b^n,
1579   b^n-c^n,
1580   c^n-d^n,
1581   d^n-e^n,
1582   a^(n-1)*b+b^(n-1)*c+c^(n-1)*d+d^(n-1)*e+e^(n-1)*a;
1583   timeStd(i,2);
1584   timeStd(i,20);
1585}
1586
1587proc factorH(poly p)
1588"USAGE:  factorH(p)  p poly
1589RETURN:  factorize(p)
1590NOTE:    changes variables to make the last variable the principal
1591         one in the multivariate factorization and factorizes then
1592         the polynomial
1593EXAMPLE: example factorH; shows an example
1594"
1595{
1596   def R=basering;
1597   int i,j;
1598   int n=1;
1599   int d=nrows(coeffs(p,var(1)));
1600   for(i=1;i<=nvars(R);i++)
1601   {
1602      j=nrows(coeffs(p,var(i)));
1603      if(d>j)
1604      {
1605         n=i;
1606         d=j;
1607      }
1608   }
1609   ideal ma=maxideal(1);    //die letzte Variable ist die Hauptvariable
1610   ma[nvars(R)]=var(n);
1611   ma[n]=var(nvars(R));
1612   map phi=R,ma;
1613   list fac=factorize(phi(p));
1614   list re=phi(fac);
1615   return(re);
1616}
1617example
1618{ "EXAMPLE:"; echo = 2;
1619  system("random",992851144);
1620  ring r=32003,(x,y,z,w,t),lp;
1621  poly p=y2w9+yz7t-yz5w4-z2w4t4-w8t3;
1622  factorize(p);  //fast
1623  system("random",992851262);
1624  //factorize(p);  //slow
1625  system("random",992851262);
1626  factorH(p);
1627}
Note: See TracBrowser for help on using the repository browser.