source: git/Singular/LIB/general.lib @ 1b2216

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