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

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