jengelh-datetimespielwiese
Last change on this file since 7f3ad4 was 7f3ad4, checked in by Hans Schönemann <hannes@…>, 14 years ago
• Property mode set to `100644`
File size: 45.4 KB
Line
2//anne, added deleteSublist and watchdog 12.12.2000
4///////////////////////////////////////////////////////////////////////////////
5version="\$Id: general.lib,v 1.59 2009-01-14 16:07:04 Singular 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
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.
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
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\"
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 if id.
745         Zero generators of ideal/module are deleted
746EXAMPLE: example sort; shows an example
747"
748{  int ii,jj,s,n = 0,0,1,0;
749   intvec v;
750   if ( defined(basering) ) { def P = basering; }
751   if ( size(#)==0 and (typeof(id)=="ideal" or typeof(id)=="module"
752                        or typeof(id)=="matrix"))
753   {
754      id = simplify(id,2);
755      for ( ii=1; ii<ncols(id); ii++ )
756      {
757         if ( id[ii]!=id[ii+1] ) { break;}
758      }
759      if ( ii != ncols(id) ) { v = sortvec(id); }
760      else  { v = ncols(id)..1; }
761   }
762   if ( size(#)>=1 and (typeof(id)=="ideal" or typeof(id)=="module"
763                        or typeof(id)=="matrix") )
764   {
765      if ( typeof(#[1])=="string" )
766      {
767         execute("ring r1 =("+charstr(P)+"),("+varstr(P)+"),("+#[1]+");");
768         def i = imap(P,id);
769         v = sortvec(i);
770         setring P;
771         n=2;
772      }
773   }
774   if ( typeof(id)=="intvec" or typeof(id)=="list" and n==0 )
775   {
776      string o;
777      if ( size(#)==0 ) { o = "lp"; n=1; }
778      if ( size(#)>=1 )
779      {
780         if ( typeof(#[1])=="string" ) { o = #[1]; n=1; }
781      }
782   }
783   if ( typeof(id)=="intvec" or typeof(id)=="list" and n==1 )
784   {
785      if ( typeof(id)=="list" )
786      {
787         for (ii=1; ii<=size(id); ii++)
788         {
789            if (typeof(id[ii]) != "intvec" and typeof(id[ii]) != "int")
790               { ERROR("// list elements must be intvec/int"); }
791            else
792               { s=size(id[ii])*(s < size(id[ii])) + s*(s >= size(id[ii])); }
793         }
794      }
795      execute("ring r=0,x(1..s),("+o+");");
796      ideal i;
797      poly f;
798      for (ii=1; ii<=size(id); ii++)
799      {
800         f=1;
801         for (jj=1; jj<=size(id[ii]); jj++)
802         {
803            f=f*x(jj)^(id[ii])[jj];
804         }
805         i[ii]=f;
806      }
807      v = sort(i)[2];
808   }
809   if ( size(#)!=0 and n==0 ) { v = #[1]; }
810   if( size(#)==2 )
811   {
812      if ( #[2] != 0 ) { v = v[size(v)..1]; }
813   }
814   s = size(v);
815   if( size(id) < s ) { s = size(id); }
816   def m = id;
817   if ( size(m) != 0 )
818   {
819      for ( jj=1; jj<=s; jj++)
820      {
821         if ( v[jj]<=0 ) { v[jj]=jj; }
822         m[jj] = id[v[jj]];
823      }
824   }
825   if ( v == 0 ) { v = 1; }
826   list L=m,v;
827   return(L);
828}
829example
830{  "EXAMPLE:"; echo = 2;
831   ring r0 = 0,(x,y,z,t),lp;
832   ideal i = x3,z3,xyz;
833   sort(i);            //sorts using lex ordering, smaller polys come first
834
835   sort(i,3..1);
836
837   sort(i,"ls")[1];     //sort w.r.t. negative lex ordering
838
839   intvec v =1,10..5,2..4;v;
840   sort(v)[1];          // sort v lexicographically
841
842   sort(v,"Dp",1)[1];   // sort v w.r.t (total sum, reverse lex)
843}
844///////////////////////////////////////////////////////////////////////////////
845
846static proc lsum (int n, list l)
847{ if (n>10)
848  { return( lsum(n/2,list(l[1..(n/2)])) + lsum(n-n/2, list(l[(n/2+1)..n])) );
849  }
850  else
851  { def Summe=l[1];
852    for (int i=2;i<=n;i++)
853    { Summe=Summe+l[i];
854    }
855    return(Summe);
856  }
857}
858
859///////////////////////////////////////////////////////////////////////////////
860
861proc sum (id, list #)
862"USAGE:    sum(id[,v]); id ideal/vector/module/matrix/intvec/intmat/list,
863          v intvec  (default: v=1..number of entries of id)
864ASSUME:   list members can be added.
865RETURN:   The sum of all entries of id [with index given by v] of type
866          depending on the entries of id.
867NOTE:     If id is not a list, id is treated as a list of polys resp. integers.
868          A module m is identified with the corresponding matrix M (columns
869          of M generate m).
870@*        If v is outside the range of id, we have the empty sum and the
871          result will be 0 (of type int).
872EXAMPLE:  example sum; shows an example
873"
874{
875//-------------------- initialization and special feature ---------------------
876   int n,j,tt;
877   string ty;                                  // will become type of id
878   list l;
879
880// We wish to allow something like sum(x(1..10)) if x(1),...,x(10) are
881// variables. x(1..10) is a list of polys and enters the procedure with
882// id=x(1) and # a list with 9 polys, #[1]= x(2),...,#[9]= x(10). Hence, in
883// this case # is never empty. If an additional intvec v is given,
884// it is added to #, so we have to separate it first and make
885// the rest a list which has to be added.
886
887   int s = size(#);
888   if( s!=0 )
889   {  if ( typeof(#[s])=="intvec" or typeof(#[s])=="int")
890      {  intvec v = #[s];
891         tt=1;
892         s=s-1;
893         if ( s>0 ) { # = #[1..s]; }
894      }
895   }
896   if ( s>0 )
897   {
898      l = list(id)+#;
899      kill id;
900      list id = l;                                    //case: id = list
901      ty = "list";
902   }
903   else
904   {
905      ty = typeof(id);
906   }
907//------------------------------ reduce to 3 cases ---------------------------
908   if( ty=="poly" or ty=="ideal" or ty=="vector"
909       or ty=="module" or ty=="matrix" )
910   {                                                 //case: id = ideal
911      ideal i = ideal(matrix(id));
912      kill id;
913      ideal id = simplify(i,2);                      //delete 0 entries
914   }
915   if( ty=="int" or ty=="intvec" or ty=="intmat" )
916   {                                                 //case: id = intvec
917      if ( ty == "int" ) { intmat S =id; }
918      else { intmat S = intmat(id); }
919      intvec i = S[1..nrows(S),1..ncols(S)];
920      kill id;
921      intvec id = i;
922   }
923//------------------- consider intvec v and empty sum  -----------------------
924   if( tt!=0 )
925   {
926      for (j=1; j<=size(v); j++)
927      {
928         if ( v[j] <= 0 or v[j] > size(id) )         //v outside range of id
929         {
930            return(0);                               //empty sum is 0
931         }
932      }
933      id = id[v];                                    //consider part of id
934   }                                                 //corresponding to v
935
937   n = size(id);
938   if (n>10)
939   { return( lsum(n/2,list(id[1..(n/2)])) + lsum(n-n/2, list(id[(n/2+1)..n])) );
940   }
941   else
942   { def Summe=id[1];
943     for (int lauf=2;lauf<=n;lauf++)
944     { Summe=Summe+id[lauf];
945     }
946     return(Summe);
947   }
948}
949example
950{  "EXAMPLE:"; echo = 2;
951   ring r1 = 0,(x,y,z),dp;
952   vector pv = [xy,xz,yz,x2,y2,z2];
953   sum(pv);
954   sum(pv,2..5);
955   matrix M[2][3] = 1,x,2,y,3,z;
956   intvec w=2,4,6;
957   sum(M,w);
958   intvec iv = 1,2,3,4,5,6,7,8,9;
959   sum(iv,2..4);
960   iv = intvec(1..100);
961   sum(iv);
962   ring r2 = 0,(x(1..10)),dp;
963   sum(x(3..7),intvec(1,3,5));
964}
965///////////////////////////////////////////////////////////////////////////////
966
967
968///////////////////////////////////////////////////////////////////////////////
969
970proc which (command)
971"USAGE:    which(command); command = string expression
972RETURN:   Absolute pathname of command, if found in search path.
973          Empty string, otherwise.
974NOTE:     Based on the Unix command 'which'.
975EXAMPLE:  example which; shows an example
976"
977{
978   int rs;
979   int i;
980   string fn = "which_" + string(system("pid"));
981   string pn;
982   string cmd;
983   if( typeof(command) != "string")
984   {
985     return (pn);
986   }
987   if (system("uname") != "ix86-Win")
988   {
989     cmd = "which ";
990   }
991   else
992   {
993     // unfortunately, it does not take -path
994     cmd = "type ";
995   }
996   i = system("sh", cmd + command + " > " + fn);
998   if (system("uname") != "ix86-Win")
999   {
1000     // TBC: Hmm... should parse output to get rid of 'command is '
1001     pn[size(pn)] = "";
1002     i = 1;
1003     while ((pn[i] != " ") and (pn[i] != ""))
1004     {
1005       i = i+1;
1006     }
1007     if (pn[i] == " ") {pn[i] = "";}
1008     rs = system("sh", "ls " + pn + " > " + fn + " 2>&1 ");
1009   }
1010   else
1011   {
1012     rs = 0;
1013   }
1014   i = system("sh", "rm " + fn);
1015   if (rs == 0) {return (pn);}
1016   else
1017   {
1019     return ("");
1020   }
1021}
1022example
1023{  "EXAMPLE:"; echo = 2;
1024    which("sh");
1025}
1026///////////////////////////////////////////////////////////////////////////////
1027
1028proc watchdog(int i, string cmd)
1029"USAGE:   watchdog(i,cmd); i integer, cmd string
1030RETURN:  Result of cmd, if the result can be computed in i seconds.
1031         Otherwise the computation is interrupted after i seconds,
1032         the string "Killed" is returned and the global variable
1033         'watchdog_interrupt' is defined.
1034NOTE:    * the MP package must be enabled
1035         * the current basering should not be watchdog_rneu, since
1036           watchdog_rneu will be killed
1037         * if there are variable names of the structure x(i) all
1038           polynomials have to be put into eval(...) in order to be
1039           interpreted correctly
1040         * a second Singular process is started by this procedure
1041EXAMPLE: example watchdog; shows an example
1042"
1043{
1044  string rname=nameof(basering);
1045  def rsave=basering;
1046  if (defined(watchdog_rneu))
1047  {
1048    kill watchdog_rneu;
1049  }
1050// If we do not have MP-links, watchdog cannot be used
1051  if (system("with","MP"))
1052  {
1053    if ( i > 0 )
1054    {
1055      int j=10;
1056      int k=999999;
1057// fork, get the pid of the child and send it the command
1059      open(l_fork);
1060      write(l_fork,quote(system("pid")));
1062      execute("write(l_fork,quote(" + cmd + "));");
1063
1064
1065// sleep in small, but growing intervals for appr. 1 second
1066      while(j < k)
1067      {
1069        j = j + j;
1070      }
1071
1072// sleep in intervals of one second
1073      j = 1;
1075      {
1076        while (j < i)
1077        {
1079          j = j + 1;
1080        }
1081      }
1082// check, whether we have a result, and return it
1084      {
1086        if (nameof(basering)!=rname)
1087        {
1088          def watchdog_rneu=basering;
1089          setring rsave;
1090          if (!defined(result))
1091          {
1092            def result=fetch(watchdog_rneu,result);
1093          }
1094        }
1095        if(defined(watchdog_interrupt))
1096        {
1097          kill watchdog_interrupt;
1098        }
1099        close(l_fork);
1100      }
1101      else
1102      {
1103        string result="Killed";
1104        if(!defined(watchdog_interrupt))
1105        {
1106          int watchdog_interrupt=1;
1107          export watchdog_interrupt;
1108        }
1109        close(l_fork);
1110        j = system("sh","kill " + string(pid));
1111      }
1112      return(result);
1113    }
1114    else
1115    {
1116      ERROR("First argument of watchdog has to be a positive integer.");
1117    }
1118  }
1119  else
1120  {
1121    ERROR("MP-support is not enabled in this version of Singular.");
1122  }
1123}
1124example
1125{ "EXAMPLE:"; echo=2;
1126  ring r=0,(x,y,z),dp;
1127  poly f=x^30+y^30;
1128  watchdog(1,"factorize(eval("+string(f)+"))");
1129  watchdog(100,"factorize(eval("+string(f)+"))");
1130}
1131///////////////////////////////////////////////////////////////////////////////
1132
1133proc deleteSublist(intvec v,list l)
1134"USAGE:   deleteSublist(v,l); intvec v; list l
1135         where the entries of the integer vector v correspond to the
1136         positions of the elements to be deleted
1137RETURN:  list without the deleted elements
1138EXAMPLE: example deleteSublist; shows an example"
1139{
1140  list k;
1141  int i,j,skip;
1142  j=1;
1143  skip=0;
1144  intvec vs=sort(v)[1];
1145  for ( i=1 ; i <=size(vs) ; i++)
1146  {
1147    while ((j+skip) < vs[i])
1148    {
1149      k[j] = l[j+skip];
1150      j++;
1151    }
1152    skip++;
1153  }
1154  if(vs[size(vs)]<size(l))
1155  {
1156    k=k+list(l[(vs[size(vs)]+1)..size(l)]);
1157  }
1158  return(k);
1159}
1160example
1161{ "EXAMPLE:"; echo=2;
1162   list l=1,2,3,4,5;
1163   intvec v=1,3,4;
1164   l=deleteSublist(v,l);
1165   l;
1166}
1167///////////////////////////////////////////////////////////////////////////////
1168proc primefactors (n, list #)
1169"USAGE:   primefactors(n [,p]); n = int or number, p = integer
1170COMPUTE: primefactors <= min(p,32003) of n (default p = 32003)
1171RETURN:  a list, say l,
1172         l[1] : primefactors <= min(p,32003) of n
1173         l[2] : l[2][i] = multiplicity of l[1][i]
1174         l[3] : remaining factor ( n=product{ (l[1][i]^l[2][i])*l[3]} )
1175                type(l[3])=typeof(n)
1176NOTE:    If n is a long integer (of type number) then the procedure
1177         finds primefactors <= min(p,32003) but n may as larger as
1178         2147483647 (max. integer representation)
1179WARNING: the procedure works for small integers only, just by testing all
1180         primes (not to be considerd as serious prime factorization!)
1181EXAMPLE: example primefactors; shows an example
1182"
1183{
1184   int ii,jj,z,p,num,w3,q;
1185   intvec w1,w2,v;
1186   list l;
1187   if (size(#) == 0)
1188   {
1189      p=32003;
1190   }
1191   else
1192   {
1193     if( typeof(#[1]) != "int")
1194     {
1195       ERROR("2nd parameter must be of type int"+newline);
1196     }
1197     p=#[1];
1198   }
1199   if( n<0) { n=-n;};
1200
1201// ----------------- case: 1st parameter is a number --------------------
1202   if (typeof(n) =="number")
1203   {
1204     kill w3;
1205     number w3;
1206     if( n > 2147483647 )           //2147483647 max. integer representation
1207     {
1208        v = primes(2,p);
1209        number m;
1210        for( ii=1; ii<=size(v); ii++)
1211        {
1212          jj=0;
1213          while(1)
1214          {
1215            q  = v[ii];
1216            jj = jj+1;
1217            m  = n/q;                  //divide n as often as possible
1218            if (denominator(m)!=1) { break;  }
1219            n=m;
1220          }
1221          if( jj>1 )
1222          {
1223             w1 = w1,v[ii];          //primes
1224             w2 = w2,jj-1;           //powers
1225          }
1226          if( n <= 2147483647 ) { break;  }
1227        }
1228     }
1229
1230     if( n >  2147483647 )            //n is still too big
1231     {
1232       if( size(w1) >1 )               //at least 1 primefactor was found
1233       {
1234         w1 = w1[2..size(w1)];
1235         w2 = w2[2..size(w2)];
1236       }
1237       else                           //no primefactor was found
1238       {
1239         w1 = 1; w2 = 1;
1240       }
1241       l  = w1,w2,n;
1242       return(l);
1243     }
1244
1245     if( n <= 2147483647 )          //n is in inter range
1246     {
1247       num = int(n);
1248       kill n;
1249       int n = num;
1250     }
1251   }
1252
1253// --------------------------- trivial cases --------------------
1254   if( n==0 )
1255   {
1256     w1=1; w2=1; w3=0; l=w1,w2,w3;
1257     return(l);
1258   }
1259
1260   if( n==1 )
1261   {
1262       w3=1;
1263       if( size(w1) >1 )               //at least 1 primefactor was found
1264       {
1265         w1 = w1[2..size(w1)];
1266         w2 = w2[2..size(w2)];
1267       }
1268       else                           //no primefactor was found
1269       {
1270         w1 = 1; w2 = 1;
1271       }
1272      l=w1,w2,w3;
1273      return(l);
1274   }
1275   if ( prime(n)==n )         //note: prime(n) <= 32003 in Singular
1276   {                          //case n is a prime
1277      if (p > n)
1278      {
1279        w1=w1,n; w2=w2,1; w3=1;
1280        w1 = w1[2..size(w1)];
1281        w2 = w2[2..size(w2)];
1282        l=w1,w2,w3;
1283        return(l);
1284      }
1285      else
1286      {
1287        w3=n;
1288        if( size(w1) >1 )               //at least 1 primefactor was found
1289        {
1290           w1 = w1[2..size(w1)];
1291           w2 = w2[2..size(w2)];
1292         }
1293         else                           //no primefactor was found
1294         {
1295           w1 = 1; w2 = 1;
1296         }
1297         l=w1,w2,w3;
1298        return(l);
1299      }
1300   }
1301   else
1302   {
1303      if ( p >= n)
1304      {
1305         v = primes(q,n div 2 + 1);
1306      }
1307      else
1308      {
1309         v = primes(q,p);
1310      }
1311//------------- search for primfactors <= last entry of v ------------
1312      for(ii=1; ii<=size(v); ii++)
1313      {
1314         z=0;
1315         while( (n mod v[ii]) == 0 )
1316         {
1317            z=z+1;
1318            n = n div v[ii];
1319         }
1320         if (z!=0)
1321         {
1322            w1 = w1,v[ii];        //primes
1323            w2 = w2,z;            //multiplicities
1324         }
1325      }
1326   }
1327//--------------- case:at least 1 primefactor was found ---------------
1328   if( size(w1) >1 )               //at least 1 primefactor was found
1329   {
1330      w1 = w1[2..size(w1)];
1331      w2 = w2[2..size(w2)];
1332   }
1333   else                           //no primefactor was found
1334   {
1335     w1 = 1; w2 = 1;
1336   }
1337   w3 = n;
1338   l  = w1,w2,w3;
1339   return(l);
1340}
1341example
1342{ "EXAMPLE:"; echo = 2;
1343   primefactors(7*8*121);
1344   ring r = 0,x,dp;
1345   primefactors(123456789100);
1346}
1347
1348///////////////////////////////////////////////////////////////////////////////
1349proc primecoeffs(J, list #)
1350"USAGE:   primecoeffs(J[,p]); J any type which can be converted to a matrix
1351         e.g. ideal, matrix, vector, module, int, intvec
1352         p = integer
1353COMPUTE: primefactors <= min(p,32003) of coeffs of J (default p = 32003)
1354RETURN:  a list, say l, of two intvectors:@*
1355         l[1] : the different primefactors of all coefficients of J@*
1356         l[2] : the different remaining factors
1357NOTE:    the procedure works for small integers only, just by testing all
1358         primes (not to be considered as serious prime factorization!)
1359EXAMPLE: example primecoeffs; shows an example
1360"
1361{
1362   int q,ii,n,mark;;
1363   if (size(#) == 0)
1364   {
1365      q=32003;
1366   }
1367   else
1368   {
1369     if( typeof(#[1]) != "int")
1370     {
1371       ERROR("2nd parameter must be of type int"+newline);
1372     }
1373     q=#[1];
1374   }
1375
1376   if (defined(basering) == 0)
1377   {
1378     mark=1;
1379     ring r = 0,x,dp;
1380   }
1381   def I = ideal(matrix(J));
1382   poly p = product(maxideal(1));
1383   matrix Coef=coef(I[1],p);
1384   ideal id, jd, rest;
1385   intvec v,re;
1386   list result,l;
1387   for(ii=2; ii<=ncols(I); ii++)
1388   {
1389     Coef=concat(Coef,coef(I[ii],p));
1390   }
1391   id = Coef[2,1..ncols(Coef)];
1392   id = simplify(id,6);
1393   for (ii=1; ii<=size(id); ii++)
1394   {
1395     l = primefactors(number(id[ii]),q);
1396     jd = jd,l[1];
1397     rest = rest,l[3];
1398   }
1399   jd = simplify(jd,6);
1400   for (ii=1; ii<=size(jd); ii++)
1401   {
1402     v[ii]=int(jd[ii]);
1403   }
1404   v = sort(v)[1];
1405   rest = simplify(rest,6);
1406   id = sort(id)[1];
1407   if (mark)
1408   {
1409     for (ii=1; ii<=size(rest); ii++)
1410     {
1411       re[ii] = int(rest[ii]);
1412     }
1413     result = v,re;
1414   }
1415   else
1416   {
1417      result = v,rest;
1418   }
1419   return(result);
1420}
1421example
1422{ "EXAMPLE:"; echo = 2;
1423   primecoeffs(intvec(7*8*121,7*8));"";
1424   ring r = 0,(b,c,t),dp;
1425   ideal I = -13b6c3t+4b5c4t,-10b4c2t-5b4ct2;
1426   primecoeffs(I);
1427}
1428///////////////////////////////////////////////////////////////////////////////
1429proc timeFactorize(poly i,list #)
1430"USAGE:  timeFactorize(p,d); poly p , integer d
1431RETURN:  factorize(p) if the factorization finished after d-1
1432         seconds otherwhise f is considered to be irreducible
1433EXAMPLE: example timeFactorize; shows an example
1434"
1435{
1436  def P=basering;
1437  if (size(#) > 0)
1438  {
1439    if (system("with", "MP"))
1440    {
1441      if ((typeof(#[1]) == "int")&&(#[1]))
1442      {
1443        int wait = #[1];
1444        int j = 10;
1445
1446        string bs = nameof(basering);
1448        open(l_fork);
1449        write(l_fork, quote(system("pid")));
1451        write(l_fork, quote(timeFactorize(eval(i))));
1452
1453        // sleep in small intervalls for appr. one second
1454        if (wait > 0)
1455        {
1456          while(j < 1000000)
1457          {
1459            j = j + j;
1460          }
1461        }
1462
1463        // sleep in intervalls of one second from now on
1464        j = 1;
1465        while (j < wait)
1466        {
1468          j = j + 1;
1469        }
1470
1472        {
1474          if (bs != nameof(basering))
1475          {
1476            def PP = basering;
1477            setring P;
1478            def result = imap(PP, result);
1479            kill PP;
1480          }
1481          kill l_fork;
1482        }
1483        else
1484        {
1485          list result;
1486          intvec v=1,1;
1487          result[1]=list(1,i);
1488          result[2]=v;
1489          j = system("sh", "kill " + string(pid));
1490        }
1491        return (result);
1492      }
1493    }
1494  }
1495  return(factorH(i));
1496}
1497example
1498{ "EXAMPLE:"; echo = 2;
1499   ring r=0,(x,y),dp;
1500   poly p=((x2+y3)^2+xy6)*((x3+y2)^2+x10y);
1501   p=p^2;
1502   //timeFactorize(p,2);
1503   //timeFactorize(p,20);
1504}
1505
1506proc timeStd(ideal i,list #)
1507"USAGE:  timeStd(i,d), i ideal, d integer
1508RETURN:  std(i) if the standard basis computation finished after
1509         d-1 seconds and i otherwhise
1510EXAMPLE: example timeStd; shows an example
1511"
1512{
1513  def P=basering;
1514  if (size(#) > 0)
1515  {
1516    if (system("with", "MP"))
1517    {
1518      if ((typeof(#[1]) == "int")&&(#[1]))
1519      {
1520        int wait = #[1];
1521        int j = 10;
1522
1523        string bs = nameof(basering);
1525        open(l_fork);
1526        write(l_fork, quote(system("pid")));
1528        write(l_fork, quote(timeStd(eval(i))));
1529
1530        // sleep in small intervalls for appr. one second
1531        if (wait > 0)
1532        {
1533          while(j < 1000000)
1534          {
1536            j = j + j;
1537          }
1538        }
1539        j = 1;
1540        while (j < wait)
1541        {
1543          j = j + 1;
1544        }
1546        {
1548          if (bs != nameof(basering))
1549          {
1550            def PP = basering;
1551            setring P;
1552            def result = imap(PP, result);
1553            kill PP;
1554          }
1555          kill l_fork;
1556        }
1557        else
1558        {
1559          ideal result=i;
1560          j = system("sh", "kill " + string(pid));
1561        }
1562        return (result);
1563      }
1564    }
1565  }
1566  return(std(i));
1567}
1568example
1569{ "EXAMPLE:"; echo = 2;
1570   ring r=32003,(a,b,c,d,e),dp;
1571   int n=6;
1572   ideal i=
1573   a^n-b^n,
1574   b^n-c^n,
1575   c^n-d^n,
1576   d^n-e^n,
1577   a^(n-1)*b+b^(n-1)*c+c^(n-1)*d+d^(n-1)*e+e^(n-1)*a;
1578   timeStd(i,2);
1579   timeStd(i,20);
1580}
1581
1582proc factorH(poly p)
1583"USAGE:  factorH(p)  p poly
1584RETURN:  factorize(p)
1585NOTE:    changes variables to make the last variable the principal
1586         one in the multivariate factorization and factorizes then
1587         the polynomial
1588EXAMPLE: example factorH; shows an example
1589"
1590{
1591   def R=basering;
1592   int i,j;
1593   int n=1;
1594   int d=nrows(coeffs(p,var(1)));
1595   for(i=1;i<=nvars(R);i++)
1596   {
1597      j=nrows(coeffs(p,var(i)));
1598      if(d>j)
1599      {
1600         n=i;
1601         d=j;
1602      }
1603   }
1604   ideal ma=maxideal(1);    //die letzte Variable ist die Hauptvariable
1605   ma[nvars(R)]=var(n);
1606   ma[n]=var(nvars(R));
1607   map phi=R,ma;
1608   list fac=factorize(phi(p));
1609   list re=phi(fac);
1610   return(re);
1611}
1612example
1613{ "EXAMPLE:"; echo = 2;
1614  system("random",992851144);
1615  ring r=32003,(x,y,z,w,t),lp;
1616  poly p=y2w9+yz7t-yz5w4-z2w4t4-w8t3;
1617  factorize(p);  //fast
1618  system("random",992851262);
1619  //factorize(p);  //slow
1620  system("random",992851262);
1621  factorH(p);
1622}
Note: See TracBrowser for help on using the repository browser.