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

spielwiese
Last change on this file since d4f44f was d4f44f, checked in by Motsak Oleksandr <motsak@…>, 15 years ago
the procedure sort is back to the original, description improvement & new example git-svn-id: file:///usr/local/Singular/svn/trunk@11980 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 45.9 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: general.lib,v 1.64 2009-07-15 14:43:52 motsak Exp $";
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  if (system("with","Namespaces"))
348  {
349    list @marie=names(Top);
350  }
351  else
352  {
353    list @marie=names();
354  }
355  int j, no_kill, @joni;
356  for ( @joni=1; @joni<=size(#);  @joni++)
357  {
358    if (typeof(#[@joni]) != "string")
359    {
360      ERROR("Need string as " + string(@joni) + "th argument");
361    }
362  }
363
364  // kills all user-defined variables but not loaded procedures
365  if( size(#)==0 )
366  {
367    for ( @joni=size(@marie); @joni>0; @joni-- )
368    {
369      if( @marie[@joni]!="LIB" and typeof(`@marie[@joni]`)!="proc"
370      and typeof(`@marie[@joni]`)!="package")
371      { kill `@marie[@joni]`; }
372    }
373  }
374  else
375  {
376    // kills all user-defined variables
377    if( size(#)==1 )
378    {
379      // of type proc
380      if( #[1] == "proc" )
381      {
382        for ( @joni=size(@marie); @joni>0; @joni-- )
383        {
384          if( (@marie[@joni]!="General")
385          and (@marie[@joni]!="Top")
386          and (@marie[@joni]!="killall")
387          and (@marie[@joni]!="LIB") and
388              ((typeof(`@marie[@joni]`)=="package") or
389              (typeof(`@marie[@joni]`)=="proc")))
390          {
391            if (defined(`@marie[@joni]`)) {kill `@marie[@joni]`;}
392          }
393          if (!defined(@joni)) break;
394        }
395        if ((system("with","Namespaces")) && defined(General))
396        {
397          @marie=names(General);
398          for ( @joni=size(@marie); @joni>0; @joni-- )
399          {
400            if(@marie[@joni]!="killall"
401            and typeof(`@marie[@joni]`)=="proc")
402            { kill General::`@marie[@joni]`; }
403          }
404          kill General::killall;
405        }
406      }
407      else
408      {
409        // other types
410        for ( @joni=size(@marie); @joni>2; @joni-- )
411        {
412          if(typeof(`@marie[@joni]`)==#[1] and @marie[@joni]!="LIB"
413             and typeof(`@marie[@joni]`)!="proc")
414          { kill `@marie[@joni]`; }
415        }
416      }
417    }
418    else
419    {
420      // kills all user-defined variables whose name or type is not #i
421      for ( @joni=size(@marie); @joni>2; @joni-- )
422      {
423        if ( @marie[@joni] != "LIB" && @marie[@joni] != "Top"
424             && typeof(`@marie[@joni]`) != "proc")
425        {
426          no_kill = 0;
427          for (j=2; j<= size(#); j++)
428          {
429            if (typeof(`@marie[@joni]`)==#[j] or @marie[@joni] == #[j])
430            {
431              no_kill = 1;
432              break;
433            }
434          }
435          if (! no_kill)
436          {
437            kill `@marie[@joni]`;
438          }
439        }
440        if (!defined(@joni)) break;
441      }
442    }
443  }
444}
445example
446{ "EXAMPLE:"; echo = 2;
447   ring rtest; ideal i=x,y,z; string str="hi"; int j = 3;
448   export rtest,i,str,j;       //this makes the local variables global
449   listvar();
450   killall("ring");            // kills all rings
451   listvar();
452   killall("not", "int");      // kills all variables except int's (and procs)
453   listvar();
454   killall();                  // kills all vars except loaded procs
455   listvar();
456}
457///////////////////////////////////////////////////////////////////////////////
458
459proc number_e (int n)
460"USAGE:   number_e(n);  n integer
461RETURN:  Euler number e=exp(1) up to n decimal digits (no rounding)
462@*       - of type string if no basering of char 0 is defined
463@*       - of type number if a basering of char 0 is defined
464DISPLAY: decimal format of e if printlevel > 0 (default:printlevel=0 )
465NOTE:    procedure uses algorithm of A.H.J. Sale
466EXAMPLE: example number_e; shows an example
467"
468{
469   int i,m,s,t;
470   intvec u,e;
471   u[n+2]=0; e[n+1]=0; e=e+1;
472   if( defined(basering) )
473   {
474      if( char(basering)==0 ) { number r=2; t=1; }
475   }
476   string result = "2.";
477   for( i=1; i<=n+1; i=i+1 )
478   {
479      e = e*10;
480      for( m=n+1; m>=1; m=m-1 )
481      {
482         s    = e[m]+u[m+1];
483         u[m] = s div (m+1);
484         e[m] = s%(m+1);
485      }
486      result = result+string(u[1]);
487      if( t==1 ) { r = r+number(u[1])/number(10)^i; }
488   }
489   if( t==1 )
490   { dbprint(printlevel-voice+2,"// "+result[1,n+1]);
491     return(r);
492   }
493   return(result[1,n+1]);
494}
495example
496{ "EXAMPLE:"; echo = 2;
497   number_e(30);"";
498   ring R = 0,t,lp;
499   number e = number_e(30);
500   e;
501}
502///////////////////////////////////////////////////////////////////////////////
503
504proc number_pi (int n)
505"USAGE:   number_pi(n);  n positive integer
506RETURN:  pi (area of unit circle) up to n decimal digits (no rounding)
507@*       - of type string if no basering of char 0 is defined,
508@*       - of type number, if a basering of char 0 is defined
509DISPLAY: decimal format of pi if printlevel > 0 (default:printlevel=0 )
510NOTE:    procedure uses algorithm of S. Rabinowitz
511EXAMPLE: example number_pi; shows an example
512"
513{
514   int i,m,t,e,q,N;
515   intvec r,p,B,Prelim;
516   string result,prelim;
517   N = (10*n) div 3 + 2;
518   p[N+1]=0; p=p+2; r=p;
519   for( i=1; i<=N+1; i=i+1 ) { B[i]=2*i-1; }
520   if( defined(basering) )
521   {
522      if( char(basering)==0 ) { number pi; number pri; t=1; }
523   }
524   for( i=0; i<=n; i=i+1 )
525   {
526      p = r*10;
527      e = p[N+1];
528      for( m=N+1; m>=2; m=m-1 )
529      {
530         r[m] = e%B[m];
531         q    = e div B[m];
532         e    = q*(m-1)+p[m-1];
533      }
534      r[1] = e%10;
535      q    = e div 10;
536      if( q!=10 and q!=9 )
537      {
538         result = result+prelim;
539         Prelim = q;
540         prelim = string(q);
541      }
542      if( q==9 )
543      {
544         Prelim = Prelim,9;
545         prelim = prelim+"9";
546      }
547      if( q==10 )
548      {
549         Prelim = (Prelim+1)-((Prelim+1) div 10)*10;
550         for( m=size(Prelim); m>0; m=m-1)
551         {
552            prelim[m] = string(Prelim[m]);
553         }
554         result = result+prelim;
555         if( t==1 ) { pi=pi+pri; }
556         Prelim = 0;
557         prelim = "0";
558      }
559      if( t==1 ) { pi=pi+number(q)/number(10)^i; }
560   }
561   result = result,prelim[1];
562   result = "3."+result[2,n-1];
563   if( t==1 )
564   { dbprint(printlevel-voice+2,"// "+result);
565     return(pi);
566   }
567   return(result);
568}
569example
570{ "EXAMPLE:"; echo = 2;
571   number_pi(11);"";
572   ring r = (real,10),t,dp;
573   number pi = number_pi(11); pi;
574}
575///////////////////////////////////////////////////////////////////////////////
576
577proc primes (int n, int m)
578"USAGE:   primes(n,m);  n,m integers
579RETURN:  intvec, consisting of all primes p, prime(n)<=p<=m, in increasing
580         order if n<=m, resp. prime(m)<=p<=n, in decreasing order if m<n.
581NOTE:    prime(n); returns the biggest prime number <= min(n,32003)
582         if n>=2, else 2
583EXAMPLE: example primes; shows an example
584"
585{  int change;
586   if ( n>m ) { change=n; n=m ; m=change; change=1; }
587   int q,p = prime(m),prime(n); intvec v = q; q = q-1;
588   while ( q>=p ) { q = prime(q); v = q,v; q = q-1; }
589   if ( change==1 ) { v = v[size(v)..1]; }
590   return(v);
591}
592example
593{  "EXAMPLE:"; echo = 2;
594    primes(50,100);"";
595    intvec v = primes(37,1); v;
596}
597///////////////////////////////////////////////////////////////////////////////
598
599proc product (id, list #)
600"USAGE:    product(id[,v]); id ideal/vector/module/matrix/intvec/intmat/list,
601          v intvec  (default: v=1..number of entries of id)
602ASSUME:   list members can be multiplied.
603RETURN:   The product of all entries of id [with index given by v] of type
604          depending on the entries of id.
605NOTE:     If id is not a list, id is treated as a list of polys resp. integers.
606          A module m is identified with the corresponding matrix M (columns
607          of M generate m).
608@*        If v is outside the range of id, we have the empty product and the
609          result will be 1 (of type int).
610EXAMPLE:  example product; shows an example
611"
612{
613//-------------------- initialization and special feature ---------------------
614   int n,j,tt;
615   string ty;                                //will become type of id
616   list l;
617
618// We wish to allow something like product(x(1..10)) if x(1),...,x(10) are
619// variables. x(1..10) is a list of polys and enters the procedure with
620// id=x(1) and # a list with 9 polys, #[1]= x(2),...,#[9]= x(10). Hence, in
621// this case # is never empty. If an additional intvec v is given,
622// it is added to #, so we have to separate it first and make
623// the rest a list which has to be multiplied.
624
625   int s = size(#);
626   if( s!=0 )
627   {  if ( typeof(#[s])=="intvec" or typeof(#[s])=="int")
628      {
629         intvec v = #[s];
630         tt=1;
631         s=s-1;
632         if ( s>0 ) { # = #[1..s]; }
633      }
634   }
635   if ( s>0 )
636   {
637      l = list(id)+#;
638      kill id;
639      list id = l;                                    //case: id = list
640      ty = "list";
641      n = size(id);
642   }
643   else
644   {
645      ty = typeof(id);
646      if( ty == "list" )
647      { n = size(id); }
648   }
649//------------------------------ reduce to 3 cases ---------------------------
650   if( ty=="poly" or ty=="ideal" or ty=="vector"
651       or ty=="module" or ty=="matrix" )
652   {
653      ideal i = ideal(matrix(id));
654      kill id;
655      ideal id = i;                                   //case: id = ideal
656      n = ncols(id);
657   }
658   if( ty=="int" or ty=="intvec" or ty=="intmat" )
659   {
660      if ( ty == "int" ) { intmat S =id; }
661      else { intmat S = intmat(id); }
662      intvec i = S[1..nrows(S),1..ncols(S)];
663      kill id;
664      intvec id = i;                                  //case: id = intvec
665      n = size(id);
666   }
667//--------------- consider intvec v and empty product  -----------------------
668   if( tt!=0 )
669   {
670      for (j=1; j<=size(v); j++)
671      {
672         if ( v[j] <= 0 or v[j] > n )                 //v outside range of id
673         {
674            return(1);                                //empty product is 1
675         }
676      }
677      id = id[v];                                     //consider part of id
678   }                                                  //corresponding to v
679//--------------------- special case: one factor is zero  ---------------------
680   if ( typeof(id) == "ideal")
681   {
682      if( size(id) < ncols(id) )
683      {
684          poly f; return(f);
685      }
686   }
687//-------------------------- finally, multiply objects  -----------------------
688   n = size(id);
689   def f(1) = id[1];
690   for( j=2; j<=n; j=j+1 ) { def f(j)=f(j-1)*id[j]; }
691   return(f(n));
692}
693example
694{  "EXAMPLE:"; echo = 2;
695   ring r= 0,(x,y,z),dp;
696   ideal m = maxideal(1);
697   product(m);
698   product(m[2..3]);
699   matrix M[2][3] = 1,x,2,y,3,z;
700   product(M);
701   intvec v=2,4,6;
702   product(M,v);
703   intvec iv = 1,2,3,4,5,6,7,8,9;
704   v=1..5,7,9;
705   product(iv,v);
706   intmat A[2][3] = 1,1,1,2,2,2;
707   product(A,3..5);
708}
709///////////////////////////////////////////////////////////////////////////////
710
711proc sort (id, list #)
712"USAGE:   sort(id[,v,o,n]); id = ideal/module/intvec/list(of intvec's or int's)
713@*       sort may be called with 1, 2 or 3 arguments in the following way:
714@*       sort(id[,v,n]); v=intvec of positive integers, n=integer,
715@*       sort(id[,o,n]); o=string (any allowed ordstr of a ring), n=integer
716RETURN:  a list l of two elements:
717@format
718        l[1]: object of same type as input but sorted in the following way:
719           - if id=ideal/module: generators of id are sorted w.r.t. intvec v
720             (id[v[1]] becomes 1-st, id[v[2]]  2-nd element, etc.). If no v is
721             present, id is sorted w.r.t. ordering o (if o is given) or w.r.t.
722             actual monomial ordering (if no o is given):
723             NOTE: generators with SMALLER(!) leading term come FIRST
724             (e.g. sort(id); sorts backwards to actual monomial ordering)
725           - if id=list of intvec's or int's: consider a list element, say
726             id[1]=3,2,5, as exponent vector of the monomial x^3*y^2*z^5;
727             the corresponding monomials are ordered w.r.t. intvec v (s.a.).
728             If no v is present, the monomials are sorted w.r.t. ordering o
729             (if o is given) or w.r.t. lexicographical ordering (if no o is
730             given). The corresponding ordered list of exponent vectors is
731             returned.
732             (e.g. sort(id); sorts lexicographically, smaller int's come first)
733             WARNING: Since negative exponents create the 0 polynomial in
734             Singular, id should not contain negative integers: the result
735             might not be as expected
736           - if id=intvec: id is treated as list of integers
737           - if n!=0 the ordering is inverse, i.e. w.r.t. v(size(v)..1)
738             default: n=0
739         l[2]: intvec, describing the permutation of the input (hence l[2]=v
740             if v is given (with positive integers))
741@end format
742NOTE:    If v is given id may be any simply indexed object (e.g. any list or
743         string); if v[i]<0 and i<=size(id) v[i] is set internally to i;
744         entries of v must be pairwise distinct to get a permutation id.
745         Zero generators of ideal/module are deleted
746         If 'o' is given, the input is sorted by considering leading terms
747         w.r.t. the new ring ordering given by 'o'
748EXAMPLE: example sort; shows an example
749"
750{  int ii,jj,s,n = 0,0,1,0;
751   intvec v;
752   if ( defined(basering) ) { def P = basering; }
753   if ( size(#)==0 and (typeof(id)=="ideal" or typeof(id)=="module"
754                        or typeof(id)=="matrix"))
755   {
756      id = simplify(id,2);
757      for ( ii=1; ii<ncols(id); ii++ )
758      {
759         if ( id[ii]!=id[ii+1] ) { break;}
760      }
761      if ( ii != ncols(id) ) { v = sortvec(id); }
762      else  { v = ncols(id)..1; }
763   }
764   if ( size(#)>=1 and (typeof(id)=="ideal" or typeof(id)=="module"
765                        or typeof(id)=="matrix") )
766   {
767      if ( typeof(#[1])=="string" )
768      {
769         execute("ring r1 =("+charstr(P)+"),("+varstr(P)+"),("+#[1]+");");
770         def i = imap(P,id);
771         v = sortvec(i);
772         setring P;
773         n=2;
774      }
775   }
776   if ( typeof(id)=="intvec" or typeof(id)=="list" and n==0 )
777   {
778      string o;
779      if ( size(#)==0 ) { o = "lp"; n=1; }
780      if ( size(#)>=1 )
781      {
782         if ( typeof(#[1])=="string" ) { o = #[1]; n=1; }
783      }
784   }
785   if ( typeof(id)=="intvec" or typeof(id)=="list" and n==1 )
786   {
787      if ( typeof(id)=="list" )
788      {
789         for (ii=1; ii<=size(id); ii++)
790         {
791            if (typeof(id[ii]) != "intvec" and typeof(id[ii]) != "int")
792               { ERROR("// list elements must be intvec/int"); }
793            else
794               { s=size(id[ii])*(s < size(id[ii])) + s*(s >= size(id[ii])); }
795         }
796      }
797      execute("ring r=0,x(1..s),("+o+");");
798      ideal i;
799      poly f;
800      for (ii=1; ii<=size(id); ii++)
801      {
802         f=1;
803         for (jj=1; jj<=size(id[ii]); jj++)
804         {
805            f=f*x(jj)^(id[ii])[jj];
806         }
807         i[ii]=f;
808      }
809      v = sort(i)[2];
810   }
811   if ( size(#)!=0 and n==0 ) { v = #[1]; }
812   if( size(#)==2 )
813   {
814      if ( #[2] != 0 ) { v = v[size(v)..1]; }
815   }
816   s = size(v);
817   if( size(id) < s ) { s = size(id); }
818   def m = id;
819   if ( size(m) != 0 )
820   {
821      for ( jj=1; jj<=s; jj++)
822      {
823         if ( v[jj]<=0 ) { v[jj]=jj; }
824         m[jj] = id[v[jj]];
825      }
826   }
827   if ( v == 0 ) { v = 1; }
828   list L=m,v;
829   return(L);
830}
831example
832{  "EXAMPLE:"; echo = 2;
833   ring r0 = 0,(x,y,z,t),lp;
834   ideal i = x3,z3,xyz;
835   sort(i);            //sorts using lex ordering, smaller polys come first
836
837   sort(i,3..1);
838
839   sort(i,"ls")[1];     //sort w.r.t. negative lex ordering
840
841   intvec v =1,10..5,2..4;v;
842   sort(v)[1];          // sort v lexicographically
843
844   sort(v,"Dp",1)[1];   // sort v w.r.t (total sum, reverse lex)
845
846   // Note that in general: lead(sort(M)) != sort(lead(M)), e.g:
847   module M = [0, 1, 1, 0], [1, 0, 0, 1]; M;
848   sort(lead(M), "c, dp")[1];
849   lead(sort(M, "c, dp")[1]); 
850
851   // In order to sort M wrt a NEW ordering by considering OLD leading
852   // terms use one of the following equivalent commands:
853   module( M[ sort(lead(M), "c,dp")[2] ] );
854   sort( M, sort(lead(M), "c,dp")[2] )[1];
855}
856///////////////////////////////////////////////////////////////////////////////
857
858static proc lsum (int n, list l)
859{ if (n>10)
860  { return( lsum(n/2,list(l[1..(n/2)])) + lsum(n-n/2, list(l[(n/2+1)..n])) );
861  }
862  else
863  { def Summe=l[1];
864    for (int i=2;i<=n;i++)
865    { Summe=Summe+l[i];
866    }
867    return(Summe);
868  }
869}
870
871///////////////////////////////////////////////////////////////////////////////
872
873proc sum (id, list #)
874"USAGE:    sum(id[,v]); id ideal/vector/module/matrix/intvec/intmat/list,
875          v intvec  (default: v=1..number of entries of id)
876ASSUME:   list members can be added.
877RETURN:   The sum of all entries of id [with index given by v] of type
878          depending on the entries of id.
879NOTE:     If id is not a list, id is treated as a list of polys resp. integers.
880          A module m is identified with the corresponding matrix M (columns
881          of M generate m).
882@*        If v is outside the range of id, we have the empty sum and the
883          result will be 0 (of type int).
884EXAMPLE:  example sum; shows an example
885"
886{
887//-------------------- initialization and special feature ---------------------
888   int n,j,tt;
889   string ty;                                  // will become type of id
890   list l;
891
892// We wish to allow something like sum(x(1..10)) if x(1),...,x(10) are
893// variables. x(1..10) is a list of polys and enters the procedure with
894// id=x(1) and # a list with 9 polys, #[1]= x(2),...,#[9]= x(10). Hence, in
895// this case # is never empty. If an additional intvec v is given,
896// it is added to #, so we have to separate it first and make
897// the rest a list which has to be added.
898
899   int s = size(#);
900   if( s!=0 )
901   {  if ( typeof(#[s])=="intvec" or typeof(#[s])=="int")
902      {  intvec v = #[s];
903         tt=1;
904         s=s-1;
905         if ( s>0 ) { # = #[1..s]; }
906      }
907   }
908   if ( s>0 )
909   {
910      l = list(id)+#;
911      kill id;
912      list id = l;                                    //case: id = list
913      ty = "list";
914   }
915   else
916   {
917      ty = typeof(id);
918   }
919//------------------------------ reduce to 3 cases ---------------------------
920   if( ty=="poly" or ty=="ideal" or ty=="vector"
921       or ty=="module" or ty=="matrix" )
922   {                                                 //case: id = ideal
923      ideal i = ideal(matrix(id));
924      kill id;
925      ideal id = simplify(i,2);                      //delete 0 entries
926   }
927   if( ty=="int" or ty=="intvec" or ty=="intmat" )
928   {                                                 //case: id = intvec
929      if ( ty == "int" ) { intmat S =id; }
930      else { intmat S = intmat(id); }
931      intvec i = S[1..nrows(S),1..ncols(S)];
932      kill id;
933      intvec id = i;
934   }
935//------------------- consider intvec v and empty sum  -----------------------
936   if( tt!=0 )
937   {
938      for (j=1; j<=size(v); j++)
939      {
940         if ( v[j] <= 0 or v[j] > size(id) )         //v outside range of id
941         {
942            return(0);                               //empty sum is 0
943         }
944      }
945      id = id[v];                                    //consider part of id
946   }                                                 //corresponding to v
947
948//-------------------------- finally, add objects  ---------------------------
949   n = size(id);
950   if (n>10)
951   { return( lsum(n/2,list(id[1..(n/2)])) + lsum(n-n/2, list(id[(n/2+1)..n])) );
952   }
953   else
954   { def Summe=id[1];
955     for (int lauf=2;lauf<=n;lauf++)
956     { Summe=Summe+id[lauf];
957     }
958     return(Summe);
959   }
960}
961example
962{  "EXAMPLE:"; echo = 2;
963   ring r1 = 0,(x,y,z),dp;
964   vector pv = [xy,xz,yz,x2,y2,z2];
965   sum(pv);
966   sum(pv,2..5);
967   matrix M[2][3] = 1,x,2,y,3,z;
968   intvec w=2,4,6;
969   sum(M,w);
970   intvec iv = 1,2,3,4,5,6,7,8,9;
971   sum(iv,2..4);
972   iv = intvec(1..100);
973   sum(iv);
974   ring r2 = 0,(x(1..10)),dp;
975   sum(x(3..7),intvec(1,3,5));
976}
977///////////////////////////////////////////////////////////////////////////////
978
979
980///////////////////////////////////////////////////////////////////////////////
981
982proc which (command)
983"USAGE:    which(command); command = string expression
984RETURN:   Absolute pathname of command, if found in search path.
985          Empty string, otherwise.
986NOTE:     Based on the Unix command 'which'.
987EXAMPLE:  example which; shows an example
988"
989{
990   int rs;
991   int i;
992   string fn = "which_" + string(system("pid"));
993   string pn;
994   string cmd;
995   if( typeof(command) != "string")
996   {
997     return (pn);
998   }
999   if (system("uname") != "ix86-Win")
1000   {
1001     cmd = "which ";
1002   }
1003   else
1004   {
1005     // unfortunately, it does not take -path
1006     cmd = "type ";
1007   }
1008   i = system("sh", cmd + command + " > " + fn);
1009   pn = read(fn);
1010   if (system("uname") != "ix86-Win")
1011   {
1012     // TBC: Hmm... should parse output to get rid of 'command is '
1013     pn[size(pn)] = "";
1014     i = 1;
1015     while ((pn[i] != " ") and (pn[i] != ""))
1016     {
1017       i = i+1;
1018     }
1019     if (pn[i] == " ") {pn[i] = "";}
1020     rs = system("sh", "ls " + pn + " > " + fn + " 2>&1 ");
1021   }
1022   else
1023   {
1024     rs = 0;
1025   }
1026   i = system("sh", "rm " + fn);
1027   if (rs == 0) {return (pn);}
1028   else
1029   {
1030     print (command + " not found ");
1031     return ("");
1032   }
1033}
1034example
1035{  "EXAMPLE:"; echo = 2;
1036    which("sh");
1037}
1038///////////////////////////////////////////////////////////////////////////////
1039
1040proc watchdog(int i, string cmd)
1041"USAGE:   watchdog(i,cmd); i integer, cmd string
1042RETURN:  Result of cmd, if the result can be computed in i seconds.
1043         Otherwise the computation is interrupted after i seconds,
1044         the string "Killed" is returned and the global variable
1045         'watchdog_interrupt' is defined.
1046NOTE:    * the MP package must be enabled
1047         * the current basering should not be watchdog_rneu, since
1048           watchdog_rneu will be killed
1049         * if there are variable names of the structure x(i) all
1050           polynomials have to be put into eval(...) in order to be
1051           interpreted correctly
1052         * a second Singular process is started by this procedure
1053EXAMPLE: example watchdog; shows an example
1054"
1055{
1056  string rname=nameof(basering);
1057  def rsave=basering;
1058  if (defined(watchdog_rneu))
1059  {
1060    kill watchdog_rneu;
1061  }
1062// If we do not have MP-links, watchdog cannot be used
1063  if (system("with","MP"))
1064  {
1065    if ( i > 0 )
1066    {
1067      int j=10;
1068      int k=999999;
1069// fork, get the pid of the child and send it the command
1070      link l_fork="MPtcp:fork";
1071      open(l_fork);
1072      write(l_fork,quote(system("pid")));
1073      int pid=read(l_fork);
1074      execute("write(l_fork,quote(" + cmd + "));");
1075
1076
1077// sleep in small, but growing intervals for appr. 1 second
1078      while(j < k)
1079      {
1080        if (status(l_fork, "read", "ready", j)) {break;}
1081        j = j + j;
1082      }
1083
1084// sleep in intervals of one second
1085      j = 1;
1086      if (!status(l_fork,"read","ready"))
1087      {
1088        while (j < i)
1089        {
1090          if (status(l_fork, "read", "ready", k)) {break;}
1091          j = j + 1;
1092        }
1093      }
1094// check, whether we have a result, and return it
1095      if (status(l_fork, "read", "ready"))
1096      {
1097        def result = read(l_fork);
1098        if (nameof(basering)!=rname)
1099        {
1100          def watchdog_rneu=basering;
1101          setring rsave;
1102          if (!defined(result))
1103          {
1104            def result=fetch(watchdog_rneu,result);
1105          }
1106        }
1107        if(defined(watchdog_interrupt))
1108        {
1109          kill watchdog_interrupt;
1110        }
1111        close(l_fork);
1112      }
1113      else
1114      {
1115        string result="Killed";
1116        if(!defined(watchdog_interrupt))
1117        {
1118          int watchdog_interrupt=1;
1119          export watchdog_interrupt;
1120        }
1121        close(l_fork);
1122        j = system("sh","kill " + string(pid));
1123      }
1124      return(result);
1125    }
1126    else
1127    {
1128      ERROR("First argument of watchdog has to be a positive integer.");
1129    }
1130  }
1131  else
1132  {
1133    ERROR("MP-support is not enabled in this version of Singular.");
1134  }
1135}
1136example
1137{ "EXAMPLE:"; echo=2;
1138  ring r=0,(x,y,z),dp;
1139  poly f=x^30+y^30;
1140  watchdog(1,"factorize(eval("+string(f)+"))");
1141  watchdog(100,"factorize(eval("+string(f)+"))");
1142}
1143///////////////////////////////////////////////////////////////////////////////
1144
1145proc deleteSublist(intvec v,list l)
1146"USAGE:   deleteSublist(v,l); intvec v; list l
1147         where the entries of the integer vector v correspond to the
1148         positions of the elements to be deleted
1149RETURN:  list without the deleted elements
1150EXAMPLE: example deleteSublist; shows an example"
1151{
1152  list k;
1153  int i,j,skip;
1154  j=1;
1155  skip=0;
1156  intvec vs=sort(v)[1];
1157  for ( i=1 ; i <=size(vs) ; i++)
1158  {
1159    while ((j+skip) < vs[i])
1160    {
1161      k[j] = l[j+skip];
1162      j++;
1163    }
1164    skip++;
1165  }
1166  if(vs[size(vs)]<size(l))
1167  {
1168    k=k+list(l[(vs[size(vs)]+1)..size(l)]);
1169  }
1170  return(k);
1171}
1172example
1173{ "EXAMPLE:"; echo=2;
1174   list l=1,2,3,4,5;
1175   intvec v=1,3,4;
1176   l=deleteSublist(v,l);
1177   l;
1178}
1179///////////////////////////////////////////////////////////////////////////////
1180proc primefactors (n, list #)
1181"USAGE:   primefactors(n [,p]); n = int or number, p = integer
1182COMPUTE: primefactors <= min(p,32003) of n (default p = 32003)
1183RETURN:  a list, say l,
1184         l[1] : primefactors <= min(p,32003) of n
1185         l[2] : l[2][i] = multiplicity of l[1][i]
1186         l[3] : remaining factor ( n=product{ (l[1][i]^l[2][i])*l[3]} )
1187                type(l[3])=typeof(n)
1188NOTE:    If n is a long integer (of type number) then the procedure
1189         finds primefactors <= min(p,32003) but n may be as larger as
1190         2147483647 (max. integer representation)
1191WARNING: the procedure works for small integers only, just by testing all
1192         primes (not to be considerd as serious prime factorization!)
1193EXAMPLE: example primefactors; shows an example
1194"
1195{
1196   int ii,jj,z,p,num,w3,q;
1197   intvec w1,w2,v;
1198   list l;
1199   if (size(#) == 0)
1200   {
1201      p=32003;
1202   }
1203   else
1204   {
1205     if( typeof(#[1]) != "int")
1206     {
1207       ERROR("2nd parameter must be of type int"+newline);
1208     }
1209     p=#[1];
1210   }
1211   if( n<0) { n=-n;};
1212
1213// ----------------- case: 1st parameter is a number --------------------
1214   if (typeof(n) =="number")
1215   {
1216     kill w3;
1217     number w3;
1218     if( n > 2147483647 )           //2147483647 max. integer representation
1219     {
1220        v = primes(2,p);
1221        number m;
1222        for( ii=1; ii<=size(v); ii++)
1223        {
1224          jj=0;
1225          while(1)
1226          {
1227            q  = v[ii];
1228            jj = jj+1;
1229            m  = n/q;                  //divide n as often as possible
1230            if (denominator(m)!=1) { break;  }
1231            n=m;
1232          }
1233          if( jj>1 )
1234          {
1235             w1 = w1,v[ii];          //primes
1236             w2 = w2,jj-1;           //powers
1237          }
1238          if( n <= 2147483647 ) { break;  }
1239        }
1240     }
1241
1242     if( n >  2147483647 )            //n is still too big
1243     {
1244       if( size(w1) >1 )               //at least 1 primefactor was found
1245       {
1246         w1 = w1[2..size(w1)];
1247         w2 = w2[2..size(w2)];
1248       }
1249       else                           //no primefactor was found
1250       {
1251         w1 = 1; w2 = 1;
1252       }
1253       l  = w1,w2,n;
1254       return(l);
1255     }
1256
1257     if( n <= 2147483647 )          //n is in inter range
1258     {
1259       num = int(n);
1260       kill n;
1261       int n = num;
1262     }
1263   }
1264
1265// --------------------------- trivial cases --------------------
1266   if( n==0 )
1267   {
1268     w1=1; w2=1; w3=0; l=w1,w2,w3;
1269     return(l);
1270   }
1271
1272   if( n==1 )
1273   {
1274       w3=1;
1275       if( size(w1) >1 )               //at least 1 primefactor was found
1276       {
1277         w1 = w1[2..size(w1)];
1278         w2 = w2[2..size(w2)];
1279       }
1280       else                           //no primefactor was found
1281       {
1282         w1 = 1; w2 = 1;
1283       }
1284      l=w1,w2,w3;
1285      return(l);
1286   }
1287   if ( prime(n)==n )         //note: prime(n) <= 32003 in Singular
1288   {                          //case n is a prime
1289      if (p > n)
1290      {
1291        w1=w1,n; w2=w2,1; w3=1;
1292        w1 = w1[2..size(w1)];
1293        w2 = w2[2..size(w2)];
1294        l=w1,w2,w3;
1295        return(l);
1296      }
1297      else
1298      {
1299        w3=n;
1300        if( size(w1) >1 )               //at least 1 primefactor was found
1301        {
1302           w1 = w1[2..size(w1)];
1303           w2 = w2[2..size(w2)];
1304         }
1305         else                           //no primefactor was found
1306         {
1307           w1 = 1; w2 = 1;
1308         }
1309         l=w1,w2,w3;
1310        return(l);
1311      }
1312   }
1313   else
1314   {
1315      if ( p >= n)
1316      {
1317         v = primes(q,n div 2 + 1);
1318      }
1319      else
1320      {
1321         v = primes(q,p);
1322      }
1323//------------- search for primfactors <= last entry of v ------------
1324      for(ii=1; ii<=size(v); ii++)
1325      {
1326         z=0;
1327         while( (n mod v[ii]) == 0 )
1328         {
1329            z=z+1;
1330            n = n div v[ii];
1331         }
1332         if (z!=0)
1333         {
1334            w1 = w1,v[ii];        //primes
1335            w2 = w2,z;            //multiplicities
1336         }
1337      }
1338   }
1339//--------------- case:at least 1 primefactor was found ---------------
1340   if( size(w1) >1 )               //at least 1 primefactor was found
1341   {
1342      w1 = w1[2..size(w1)];
1343      w2 = w2[2..size(w2)];
1344   }
1345   else                           //no primefactor was found
1346   {
1347     w1 = 1; w2 = 1;
1348   }
1349   w3 = n;
1350   l  = w1,w2,w3;
1351   return(l);
1352}
1353example
1354{ "EXAMPLE:"; echo = 2;
1355   primefactors(7*8*121);
1356   ring r = 0,x,dp;
1357   primefactors(123456789100);
1358}
1359
1360///////////////////////////////////////////////////////////////////////////////
1361proc primecoeffs(J, list #)
1362"USAGE:   primecoeffs(J[,p]); J any type which can be converted to a matrix
1363         e.g. ideal, matrix, vector, module, int, intvec
1364         p = integer
1365COMPUTE: primefactors <= min(p,32003) of coeffs of J (default p = 32003)
1366RETURN:  a list, say l, of two intvectors:@*
1367         l[1] : the different primefactors of all coefficients of J@*
1368         l[2] : the different remaining factors
1369NOTE:    the procedure works for small integers only, just by testing all
1370         primes (not to be considered as serious prime factorization!)
1371EXAMPLE: example primecoeffs; shows an example
1372"
1373{
1374   int q,ii,n,mark;;
1375   if (size(#) == 0)
1376   {
1377      q=32003;
1378   }
1379   else
1380   {
1381     if( typeof(#[1]) != "int")
1382     {
1383       ERROR("2nd parameter must be of type int"+newline);
1384     }
1385     q=#[1];
1386   }
1387
1388   if (defined(basering) == 0)
1389   {
1390     mark=1;
1391     ring r = 0,x,dp;
1392   }
1393   def I = ideal(matrix(J));
1394   poly p = product(maxideal(1));
1395   matrix Coef=coef(I[1],p);
1396   ideal id, jd, rest;
1397   intvec v,re;
1398   list result,l;
1399   for(ii=2; ii<=ncols(I); ii++)
1400   {
1401     Coef=concat(Coef,coef(I[ii],p));
1402   }
1403   id = Coef[2,1..ncols(Coef)];
1404   id = simplify(id,6);
1405   for (ii=1; ii<=size(id); ii++)
1406   {
1407     l = primefactors(number(id[ii]),q);
1408     jd = jd,l[1];
1409     rest = rest,l[3];
1410   }
1411   jd = simplify(jd,6);
1412   for (ii=1; ii<=size(jd); ii++)
1413   {
1414     v[ii]=int(jd[ii]);
1415   }
1416   v = sort(v)[1];
1417   rest = simplify(rest,6);
1418   id = sort(id)[1];
1419   if (mark)
1420   {
1421     for (ii=1; ii<=size(rest); ii++)
1422     {
1423       re[ii] = int(rest[ii]);
1424     }
1425     result = v,re;
1426   }
1427   else
1428   {
1429      result = v,rest;
1430   }
1431   return(result);
1432}
1433example
1434{ "EXAMPLE:"; echo = 2;
1435   primecoeffs(intvec(7*8*121,7*8));"";
1436   ring r = 0,(b,c,t),dp;
1437   ideal I = -13b6c3t+4b5c4t,-10b4c2t-5b4ct2;
1438   primecoeffs(I);
1439}
1440///////////////////////////////////////////////////////////////////////////////
1441proc timeFactorize(poly i,list #)
1442"USAGE:  timeFactorize(p,d); poly p , integer d
1443RETURN:  factorize(p) if the factorization finished after d-1
1444         seconds otherwhise f is considered to be irreducible
1445EXAMPLE: example timeFactorize; shows an example
1446"
1447{
1448  def P=basering;
1449  if (size(#) > 0)
1450  {
1451    if (system("with", "MP"))
1452    {
1453      if ((typeof(#[1]) == "int")&&(#[1]))
1454      {
1455        int wait = #[1];
1456        int j = 10;
1457
1458        string bs = nameof(basering);
1459        link l_fork = "MPtcp:fork";
1460        open(l_fork);
1461        write(l_fork, quote(system("pid")));
1462        int pid = read(l_fork);
1463        write(l_fork, quote(timeFactorize(eval(i))));
1464
1465        // sleep in small intervalls for appr. one second
1466        if (wait > 0)
1467        {
1468          while(j < 1000000)
1469          {
1470            if (status(l_fork, "read", "ready", j)) {break;}
1471            j = j + j;
1472          }
1473        }
1474
1475        // sleep in intervalls of one second from now on
1476        j = 1;
1477        while (j < wait)
1478        {
1479          if (status(l_fork, "read", "ready", 1000000)) {break;}
1480          j = j + 1;
1481        }
1482
1483        if (status(l_fork, "read", "ready"))
1484        {
1485          def result = read(l_fork);
1486          if (bs != nameof(basering))
1487          {
1488            def PP = basering;
1489            setring P;
1490            def result = imap(PP, result);
1491            kill PP;
1492          }
1493          kill l_fork;
1494        }
1495        else
1496        {
1497          list result;
1498          intvec v=1,1;
1499          result[1]=list(1,i);
1500          result[2]=v;
1501          j = system("sh", "kill " + string(pid));
1502        }
1503        return (result);
1504      }
1505    }
1506  }
1507  return(factorH(i));
1508}
1509example
1510{ "EXAMPLE:"; echo = 2;
1511   ring r=0,(x,y),dp;
1512   poly p=((x2+y3)^2+xy6)*((x3+y2)^2+x10y);
1513   p=p^2;
1514   //timeFactorize(p,2);
1515   //timeFactorize(p,20);
1516}
1517
1518proc timeStd(ideal i,list #)
1519"USAGE:  timeStd(i,d), i ideal, d integer
1520RETURN:  std(i) if the standard basis computation finished after
1521         d-1 seconds and i otherwhise
1522EXAMPLE: example timeStd; shows an example
1523"
1524{
1525  def P=basering;
1526  if (size(#) > 0)
1527  {
1528    if (system("with", "MP"))
1529    {
1530      if ((typeof(#[1]) == "int")&&(#[1]))
1531      {
1532        int wait = #[1];
1533        int j = 10;
1534
1535        string bs = nameof(basering);
1536        link l_fork = "MPtcp:fork";
1537        open(l_fork);
1538        write(l_fork, quote(system("pid")));
1539        int pid = read(l_fork);
1540        write(l_fork, quote(timeStd(eval(i))));
1541
1542        // sleep in small intervalls for appr. one second
1543        if (wait > 0)
1544        {
1545          while(j < 1000000)
1546          {
1547            if (status(l_fork, "read", "ready", j)) {break;}
1548            j = j + j;
1549          }
1550        }
1551        j = 1;
1552        while (j < wait)
1553        {
1554          if (status(l_fork, "read", "ready", 1000000)) {break;}
1555          j = j + 1;
1556        }
1557        if (status(l_fork, "read", "ready"))
1558        {
1559          def result = read(l_fork);
1560          if (bs != nameof(basering))
1561          {
1562            def PP = basering;
1563            setring P;
1564            def result = imap(PP, result);
1565            kill PP;
1566          }
1567          kill l_fork;
1568        }
1569        else
1570        {
1571          ideal result=i;
1572          j = system("sh", "kill " + string(pid));
1573        }
1574        return (result);
1575      }
1576    }
1577  }
1578  return(std(i));
1579}
1580example
1581{ "EXAMPLE:"; echo = 2;
1582   ring r=32003,(a,b,c,d,e),dp;
1583   int n=6;
1584   ideal i=
1585   a^n-b^n,
1586   b^n-c^n,
1587   c^n-d^n,
1588   d^n-e^n,
1589   a^(n-1)*b+b^(n-1)*c+c^(n-1)*d+d^(n-1)*e+e^(n-1)*a;
1590   timeStd(i,2);
1591   timeStd(i,20);
1592}
1593
1594proc factorH(poly p)
1595"USAGE:  factorH(p)  p poly
1596RETURN:  factorize(p)
1597NOTE:    changes variables to make the last variable the principal
1598         one in the multivariate factorization and factorizes then
1599         the polynomial
1600EXAMPLE: example factorH; shows an example
1601"
1602{
1603   def R=basering;
1604   int i,j;
1605   int n=1;
1606   int d=nrows(coeffs(p,var(1)));
1607   for(i=1;i<=nvars(R);i++)
1608   {
1609      j=nrows(coeffs(p,var(i)));
1610      if(d>j)
1611      {
1612         n=i;
1613         d=j;
1614      }
1615   }
1616   ideal ma=maxideal(1);    //die letzte Variable ist die Hauptvariable
1617   ma[nvars(R)]=var(n);
1618   ma[n]=var(nvars(R));
1619   map phi=R,ma;
1620   list fac=factorize(phi(p));
1621   list re=phi(fac);
1622   return(re);
1623}
1624example
1625{ "EXAMPLE:"; echo = 2;
1626  system("random",992851144);
1627  ring r=32003,(x,y,z,w,t),lp;
1628  poly p=y2w9+yz7t-yz5w4-z2w4t4-w8t3;
1629  factorize(p);  //fast
1630  system("random",992851262);
1631  //factorize(p);  //slow
1632  system("random",992851262);
1633  factorH(p);
1634}
Note: See TracBrowser for help on using the repository browser.