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

spielwiese
Last change on this file since c94e60f was c94e60f, checked in by Hans Schoenemann <hannes@…>, 14 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
RevLine 
[d694de]1//GMG, last modified 18.6.99
[ebbe4a]2//anne, added deleteSublist and watchdog 12.12.2000
[ca41246]3//eric, added absValue 11.04.2002
[3d124a7]4///////////////////////////////////////////////////////////////////////////////
[341696]5version="$Id$";
[49998f]6category="General purpose";
[5480da]7info="
[803c5a1]8LIBRARY:  general.lib   Elementary Computations of General Type
[3d124a7]9
[f34c37c]10PROCEDURES:
[0b59f5]11 A_Z(\"a\",n);          string a,b,... of n comma separated letters
[63be42]12 ASCII([n,m]);          string of printable ASCII characters (number n to m)
[ca41246]13 absValue(c);           absolute value of c
[fc62b67]14 binomial(n,m[,../..]); n choose m (type int), [type bigint]
[ebbe4a]15 deleteSublist(iv,l);   delete entries given by iv from list l
[fc62b67]16 factorial(n[,../..]);  n factorial (=n!) (type int), [type bigint]
[3d124a7]17 fibonacci(n[,p]);      nth Fibonacci number [char p]
[dd2aa36]18 kmemory([n[,v]]);      active [allocated] memory in kilobyte
[3d124a7]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]
[ebbe4a]26 watchdog(i,cmd);       only wait for result of command cmd for i seconds
[63be42]27 which(command);        search for command and return absolute path, if found
[298d0a]28 primecoeffs(J[,q]);    primefactors <= min(p,32003) of coeffs of J
[90d772]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
[194f5e5]37           (parameters in square brackets [] are optional)
[5480da]38";
39
[3d124a7]40LIB "inout.lib";
[8b87364]41LIB "poly.lib";
42LIB "matrix.lib";
[3d124a7]43///////////////////////////////////////////////////////////////////////////////
44
45proc A_Z (string s,int n)
[d2b2a7]46"USAGE:   A_Z(\"a\",n);  a any letter, n integer (-26<= n <=26, !=0)
[3d124a7]47RETURN:  string of n small (if a is small) or capital (if a is capital)
[0b59f5]48         letters, comma separated, beginning with a, in alphabetical
[a7a00b]49         order (or reverse alphabetical order if n<0)
[3d124a7]50EXAMPLE: example A_Z; shows an example
[d2b2a7]51"
[3d124a7]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;
[034ce1]96   execute(sR);
[3d124a7]97   R;
98}
99///////////////////////////////////////////////////////////////////////////////
[63be42]100proc ASCII (list #)
101"USAGE:   ASCII([n,m]); n,m= integers (32 <= n <= m <= 126)
[a7a00b]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@*
[b42ab6]105          ASCII(n,m): n-th up to m-th ASCII character (inclusive)
[63be42]106EXAMPLE: example ASCII; shows an example
107"
108{
109  string s1 =
[008846]110 "     !    \"    #    $    %    &    '    (    )    *    +    ,    -    .
[63be42]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///////////////////////////////////////////////////////////////////////////////
[3d124a7]148
[ca41246]149proc absValue(def c)
[298d0a]150"USAGE:  absValue(c); c int, number or poly
[ca41246]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
[298d0a]154@*       operators are defined.
[ca41246]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
[fc62b67]171proc binomial (int n, int k)
172"USAGE:   binomial(n,k); n,k integers
[b42ab6]173RETURN:  binomial(n,k); binomial coefficient n choose k
[fc62b67]174@*       - of type bigint (computed in characteristic 0)
[c860e9]175NOTE:    In any characteristic, binomial(n,k) = coefficient of x^k in (1+x)^n
[b42ab6]176SEE ALSO: prime
[3d124a7]177EXAMPLE: example binomial; shows an example
[d2b2a7]178"
[3d124a7]179{
[fc62b67]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 )
[c860e9]189   {
[fc62b67]190      r=r*(nn+1-l)/l;
[c860e9]191   }
[fc62b67]192   return(r);
193}
[3d124a7]194example
195{ "EXAMPLE:"; echo = 2;
[883957]196   binomial(200,100);"";                 //type bigint
[c860e9]197   int n,k = 200,100;
[7f3ad4]198   bigint b1 = binomial(n,k);
[c860e9]199   ring r = 0,x,dp;
[883957]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
[c860e9]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);
[3d124a7]242}
243///////////////////////////////////////////////////////////////////////////////
244
[fc62b67]245proc factorial (int n)
246"USAGE:    factorial(n);  n integer
247RETURN:   factorial(n):   n! of type bigint.
[b42ab6]248SEE ALSO: prime
249EXAMPLE:  example factorial; shows an example
[d2b2a7]250"
[fc62b67]251{
252   bigint r=1;
[c860e9]253//------------------------------ computation --------------------------------
[a292211]254   for (int l=2; l<=n; l++)
[3d124a7]255   {
[f937e2]256      r=r*l;
[3d124a7]257   }
[fc62b67]258   return(r);
[3d124a7]259}
260example
261{ "EXAMPLE:"; echo = 2;
[fc62b67]262   factorial(37);"";                 //37! (as long integer)
[3d124a7]263}
264///////////////////////////////////////////////////////////////////////////////
265
[fc62b67]266proc fibonacci (int n)
267"USAGE:    fibonacci(n);  n,p integers
[b42ab6]268RETURN:   fibonacci(n): nth Fibonacci number, f(0)=f(1)=1, f(i+1)=f(i-1)+f(i)
[fc62b67]269@*        - computed in characteristic 0, of type bigint
[b42ab6]270SEE ALSO: prime
[3d124a7]271EXAMPLE: example fibonacci; shows an example
[d2b2a7]272"
[fc62b67]273{
274  bigint f,g,h=1,1,1;
[c860e9]275//------------------------------ computation --------------------------------
[f937e2]276   for (ii=3; ii<=n; ii=ii+1)
[3d124a7]277   {
[f937e2]278      h = f+g; f = g; g = h;
[fc62b67]279   }
280   return(h);
[3d124a7]281}
282example
283{ "EXAMPLE:"; echo = 2;
[b42ab6]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
[c860e9]287   b;
[3d124a7]288}
289///////////////////////////////////////////////////////////////////////////////
290
[d694de]291proc kmemory (list #)
[b42ab6]292"USAGE:   kmemory([n,[v]]); n,v integers
[a7a00b]293RETURN:  memory in kilobyte of type bigint
[b42ab6]294@*       n=0: memory used by active variables (same as no parameters)
295@*       n=1: total memory allocated by Singular
[dd2aa36]296DISPLAY: detailed information about allocated and used memory if v!=0
[d694de]297NOTE:    kmemory uses internal function 'memory' to compute kilobyte, and
298         is the same as 'memory' for n!=0,1,2
[3d124a7]299EXAMPLE: example kmemory; shows an example
[d2b2a7]300"
[917fb5]301{
[d694de]302   int n;
[dd2aa36]303   int verb;
[d694de]304   if (size(#) != 0)
305   {
306     n=#[1];
[dd2aa36]307     if (size(#) >1)
308     { verb=#[2]; }
[d694de]309   }
[917fb5]310
[dd2aa36]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):"); }
[c60d60]319  }
320  return ((memory(n)+1023) div 1024);
[3d124a7]321}
322example
323{ "EXAMPLE:"; echo = 2;
324   kmemory();
[dd2aa36]325   kmemory(1,1);
[3d124a7]326}
327///////////////////////////////////////////////////////////////////////////////
328
329proc killall
[d2b2a7]330"USAGE:   killall(); (no parameter)
331         killall(\"type_name\");
332         killall(\"not\", \"type_name\");
[b42ab6]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
[65546eb]339@*       - killall(\"not\", \"name_1\", \"name_2\", ...);
340           kills all user-defined variables, except those of name \"name_i\"
[b42ab6]341           and except loaded procedures
[3d124a7]342NOTE:    killall should never be used inside a procedure
343EXAMPLE: example killall; shows an example AND KILLS ALL YOUR VARIABLES
[d2b2a7]344"
[3d124a7]345{
[01cda73]346  list @marie=names(Top);
[09f420]347  int j, no_kill, @joni;
[b42ab6]348  for ( @joni=1; @joni<=size(#);  @joni++)
[5c187b]349  {
[b42ab6]350    if (typeof(#[@joni]) != "string")
[5c187b]351    {
[b42ab6]352      ERROR("Need string as " + string(@joni) + "th argument");
[5c187b]353    }
354  }
[65546eb]355
[5c187b]356  // kills all user-defined variables but not loaded procedures
357  if( size(#)==0 )
358  {
[b42ab6]359    for ( @joni=size(@marie); @joni>0; @joni-- )
[5c187b]360    {
[48c165a]361      if( @marie[@joni]!="LIB" and typeof(`@marie[@joni]`)!="proc"
362      and typeof(`@marie[@joni]`)!="package")
[b42ab6]363      { kill `@marie[@joni]`; }
[5c187b]364    }
365  }
366  else
367  {
368    // kills all user-defined variables
369    if( size(#)==1 )
370    {
371      // of type proc
372      if( #[1] == "proc" )
[3d124a7]373      {
[b42ab6]374        for ( @joni=size(@marie); @joni>0; @joni-- )
[5c187b]375        {
[298d0a]376          if( (@marie[@joni]!="General")
[48c165a]377          and (@marie[@joni]!="Top")
378          and (@marie[@joni]!="killall")
[298d0a]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;
[48c165a]386        }
[01cda73]387        if (defined(General))
[48c165a]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;
[5c187b]397        }
[3d124a7]398      }
[5c187b]399      else
[65546eb]400      {
[5c187b]401        // other types
[b42ab6]402        for ( @joni=size(@marie); @joni>2; @joni-- )
[5c187b]403        {
[65546eb]404          if(typeof(`@marie[@joni]`)==#[1] and @marie[@joni]!="LIB"
405             and typeof(`@marie[@joni]`)!="proc")
[b42ab6]406          { kill `@marie[@joni]`; }
[5c187b]407        }
408      }
409    }
410    else
411    {
[65546eb]412      // kills all user-defined variables whose name or type is not #i
[b42ab6]413      for ( @joni=size(@marie); @joni>2; @joni-- )
[5c187b]414      {
[a1b1dd]415        if ( @marie[@joni] != "LIB" && @marie[@joni] != "Top"
416             && typeof(`@marie[@joni]`) != "proc")
[5c187b]417        {
418          no_kill = 0;
[09f420]419          for (j=2; j<= size(#); j++)
[6f2edc]420          {
[b42ab6]421            if (typeof(`@marie[@joni]`)==#[j] or @marie[@joni] == #[j])
[5c187b]422            {
423              no_kill = 1;
424              break;
425            }
[6f2edc]426          }
[5c187b]427          if (! no_kill)
[6f2edc]428          {
[b42ab6]429            kill `@marie[@joni]`;
[6f2edc]430          }
431        }
[298d0a]432        if (!defined(@joni)) break;
[5c187b]433      }
434    }
[6f2edc]435  }
[3d124a7]436}
437example
438{ "EXAMPLE:"; echo = 2;
[c860e9]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();
[3d124a7]448}
449///////////////////////////////////////////////////////////////////////////////
450
451proc number_e (int n)
[d2b2a7]452"USAGE:   number_e(n);  n integer
[b42ab6]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
[3d124a7]458EXAMPLE: example number_e; shows an example
[d2b2a7]459"
[3d124a7]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];
[18dd47]475         u[m] = s div (m+1);
[3d124a7]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   }
[c860e9]481   if( t==1 )
482   { dbprint(printlevel-voice+2,"// "+result[1,n+1]);
483     return(r);
484   }
[3d124a7]485   return(result[1,n+1]);
486}
487example
488{ "EXAMPLE:"; echo = 2;
[c860e9]489   number_e(30);"";
[3d124a7]490   ring R = 0,t,lp;
[c860e9]491   number e = number_e(30);
[3d124a7]492   e;
493}
494///////////////////////////////////////////////////////////////////////////////
495
496proc number_pi (int n)
[d2b2a7]497"USAGE:   number_pi(n);  n positive integer
[b42ab6]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
[3d124a7]503EXAMPLE: example number_pi; shows an example
[d2b2a7]504"
[3d124a7]505{
506   int i,m,t,e,q,N;
507   intvec r,p,B,Prelim;
508   string result,prelim;
[18dd47]509   N = (10*n) div 3 + 2;
[3d124a7]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];
[18dd47]523         q    = e div B[m];
[3d124a7]524         e    = q*(m-1)+p[m-1];
525      }
526      r[1] = e%10;
[18dd47]527      q    = e div 10;
[3d124a7]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      {
[18dd47]541         Prelim = (Prelim+1)-((Prelim+1) div 10)*10;
[3d124a7]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];
[c860e9]555   if( t==1 )
556   { dbprint(printlevel-voice+2,"// "+result);
557     return(pi);
558   }
[3d124a7]559   return(result);
560}
561example
562{ "EXAMPLE:"; echo = 2;
[c860e9]563   number_pi(11);"";
564   ring r = (real,10),t,dp;
565   number pi = number_pi(11); pi;
[3d124a7]566}
567///////////////////////////////////////////////////////////////////////////////
568
569proc primes (int n, int m)
[d2b2a7]570"USAGE:   primes(n,m);  n,m integers
[3d124a7]571RETURN:  intvec, consisting of all primes p, prime(n)<=p<=m, in increasing
[b42ab6]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
[3d124a7]575EXAMPLE: example primes; shows an example
[d2b2a7]576"
[3d124a7]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;
[c860e9]586    primes(50,100);"";
587    intvec v = primes(37,1); v;
[3d124a7]588}
589///////////////////////////////////////////////////////////////////////////////
590
591proc product (id, list #)
[c860e9]592"USAGE:    product(id[,v]); id ideal/vector/module/matrix/intvec/intmat/list,
[b42ab6]593          v intvec  (default: v=1..number of entries of id)
594ASSUME:   list members can be multiplied.
[65546eb]595RETURN:   The product of all entries of id [with index given by v] of type
[b42ab6]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).
[7708934]600@*        If v is outside the range of id, we have the empty product and the
601          result will be 1 (of type int).
[3d124a7]602EXAMPLE:  example product; shows an example
[d2b2a7]603"
[65546eb]604{
[7708934]605//-------------------- initialization and special feature ---------------------
[c860e9]606   int n,j,tt;
[65546eb]607   string ty;                                //will become type of id
[c860e9]608   list l;
[7708934]609
610// We wish to allow something like product(x(1..10)) if x(1),...,x(10) are
[65546eb]611// variables. x(1..10) is a list of polys and enters the procedure with
[7708934]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
[c860e9]617   int s = size(#);
618   if( s!=0 )
[65546eb]619   {  if ( typeof(#[s])=="intvec" or typeof(#[s])=="int")
620      {
[7708934]621         intvec v = #[s];
[65546eb]622         tt=1;
[7708934]623         s=s-1;
[c860e9]624         if ( s>0 ) { # = #[1..s]; }
625      }
626   }
627   if ( s>0 )
628   {
[7708934]629      l = list(id)+#;
630      kill id;
631      list id = l;                                    //case: id = list
632      ty = "list";
633      n = size(id);
[c860e9]634   }
635   else
[65546eb]636   {
[7708934]637      ty = typeof(id);
[65546eb]638      if( ty == "list" )
[7708934]639      { n = size(id); }
[c860e9]640   }
[7708934]641//------------------------------ reduce to 3 cases ---------------------------
[c860e9]642   if( ty=="poly" or ty=="ideal" or ty=="vector"
643       or ty=="module" or ty=="matrix" )
[3d124a7]644   {
645      ideal i = ideal(matrix(id));
[c860e9]646      kill id;
[7708934]647      ideal id = i;                                   //case: id = ideal
648      n = ncols(id);
[3d124a7]649   }
[c860e9]650   if( ty=="int" or ty=="intvec" or ty=="intmat" )
[3d124a7]651   {
[c860e9]652      if ( ty == "int" ) { intmat S =id; }
[3d124a7]653      else { intmat S = intmat(id); }
654      intvec i = S[1..nrows(S),1..ncols(S)];
[c860e9]655      kill id;
[7708934]656      intvec id = i;                                  //case: id = intvec
[65546eb]657      n = size(id);
[7708934]658   }
659//--------------- consider intvec v and empty product  -----------------------
[65546eb]660   if( tt!=0 )
[7708934]661   {
662      for (j=1; j<=size(v); j++)
663      {
664         if ( v[j] <= 0 or v[j] > n )                 //v outside range of id
[65546eb]665         {
[7708934]666            return(1);                                //empty product is 1
[65546eb]667         }
[7708934]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      }
[3d124a7]678   }
[7708934]679//-------------------------- finally, multiply objects  -----------------------
[65546eb]680   n = size(id);
[7708934]681   def f(1) = id[1];
[c860e9]682   for( j=2; j<=n; j=j+1 ) { def f(j)=f(j-1)*id[j]; }
683   return(f(n));
[3d124a7]684}
685example
686{  "EXAMPLE:"; echo = 2;
687   ring r= 0,(x,y,z),dp;
688   ideal m = maxideal(1);
689   product(m);
[c860e9]690   product(m[2..3]);
[3d124a7]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 #)
[a7a00b]704"USAGE:   sort(id[,v,o,n]); id = ideal/module/intvec/list(of intvec's or int's)
[b42ab6]705@*       sort may be called with 1, 2 or 3 arguments in the following way:
[6e52cf]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
[b42ab6]708RETURN:  a list l of two elements:
709@format
710        l[1]: object of same type as input but sorted in the following way:
[3d124a7]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):
[b42ab6]715             NOTE: generators with SMALLER(!) leading term come FIRST
716             (e.g. sort(id); sorts backwards to actual monomial ordering)
[3d124a7]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)
[a30caa3]725             WARNING: Since negative exponents create the 0 polynomial in
[63be42]726             Singular, id should not contain negative integers: the result
[a30caa3]727             might not be as expected
[3d124a7]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
[b42ab6]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
[63be42]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;
[f32177]736         entries of v must be pairwise distinct to get a permutation id.
[3d124a7]737         Zero generators of ideal/module are deleted
[d4f44f]738         If 'o' is given, the input is sorted by considering leading terms
739         w.r.t. the new ring ordering given by 'o'
[3d124a7]740EXAMPLE: example sort; shows an example
[d2b2a7]741"
[c860e9]742{  int ii,jj,s,n = 0,0,1,0;
[3d124a7]743   intvec v;
744   if ( defined(basering) ) { def P = basering; }
[b5726c]745   if ( size(#)==0 and (typeof(id)=="ideal" or typeof(id)=="module"
746                        or typeof(id)=="matrix"))
[3d124a7]747   {
748      id = simplify(id,2);
[558209]749      for ( ii=1; ii<ncols(id); ii++ )
[3d124a7]750      {
751         if ( id[ii]!=id[ii+1] ) { break;}
752      }
[558209]753      if ( ii != ncols(id) ) { v = sortvec(id); }
754      else  { v = ncols(id)..1; }
[3d124a7]755   }
[b5726c]756   if ( size(#)>=1 and (typeof(id)=="ideal" or typeof(id)=="module"
757                        or typeof(id)=="matrix") )
[3d124a7]758   {
759      if ( typeof(#[1])=="string" )
760      {
[034ce1]761         execute("ring r1 =("+charstr(P)+"),("+varstr(P)+"),("+#[1]+");");
[d4f44f]762         def i = imap(P,id);
[3d124a7]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")
[bea07f]784               { ERROR("// list elements must be intvec/int"); }
[3d124a7]785            else
786               { s=size(id[ii])*(s < size(id[ii])) + s*(s >= size(id[ii])); }
787         }
788      }
[034ce1]789      execute("ring r=0,x(1..s),("+o+");");
[3d124a7]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);
[63be42]809   if( size(id) < s ) { s = size(id); }
[3d124a7]810   def m = id;
[63be42]811   if ( size(m) != 0 )
812   {
[558209]813      for ( jj=1; jj<=s; jj++)
[63be42]814      {
815         if ( v[jj]<=0 ) { v[jj]=jj; }
816         m[jj] = id[v[jj]];
817      }
818   }
819   if ( v == 0 ) { v = 1; }
[3d124a7]820   list L=m,v;
821   return(L);
822}
823example
824{  "EXAMPLE:"; echo = 2;
[c860e9]825   ring r0 = 0,(x,y,z,t),lp;
826   ideal i = x3,z3,xyz;
[584f84d]827   sort(i);            //sorts using lex ordering, smaller polys come first
[65546eb]828
[c860e9]829   sort(i,3..1);
[b42ab6]830
[584f84d]831   sort(i,"ls")[1];     //sort w.r.t. negative lex ordering
[b42ab6]832
833   intvec v =1,10..5,2..4;v;
[584f84d]834   sort(v)[1];          // sort v lexicographically
[b42ab6]835
[584f84d]836   sort(v,"Dp",1)[1];   // sort v w.r.t (total sum, reverse lex)
[d4f44f]837
838   // Note that in general: lead(sort(M)) != sort(lead(M)), e.g:
839   module M = [0, 1, 1, 0], [1, 0, 0, 1]; M;
[3f4e52]840   sort(lead(M), "c, dp")[1];
841   lead(sort(M, "c, dp")[1]);
[d4f44f]842
843   // In order to sort M wrt a NEW ordering by considering OLD leading
844   // terms use one of the following equivalent commands:
[3f4e52]845   module( M[ sort(lead(M), "c,dp")[2] ] );
846   sort( M, sort(lead(M), "c,dp")[2] )[1];
[3d124a7]847}
848///////////////////////////////////////////////////////////////////////////////
[15e59a]849
[84375a]850static proc lsum (int n, list l)
[15e59a]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
[3d124a7]865proc sum (id, list #)
[b42ab6]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.
[65546eb]869RETURN:   The sum of all entries of id [with index given by v] of type
[b42ab6]870          depending on the entries of id.
[7708934]871NOTE:     If id is not a list, id is treated as a list of polys resp. integers.
[b42ab6]872          A module m is identified with the corresponding matrix M (columns
873          of M generate m).
[7708934]874@*        If v is outside the range of id, we have the empty sum and the
875          result will be 0 (of type int).
[3d124a7]876EXAMPLE:  example sum; shows an example
[d2b2a7]877"
[3d124a7]878{
[7708934]879//-------------------- initialization and special feature ---------------------
[b42ab6]880   int n,j,tt;
[7708934]881   string ty;                                  // will become type of id
[b42ab6]882   list l;
[7708934]883
884// We wish to allow something like sum(x(1..10)) if x(1),...,x(10) are
[65546eb]885// variables. x(1..10) is a list of polys and enters the procedure with
[7708934]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
[b42ab6]891   int s = size(#);
892   if( s!=0 )
[7708934]893   {  if ( typeof(#[s])=="intvec" or typeof(#[s])=="int")
[b42ab6]894      {  intvec v = #[s];
[65546eb]895         tt=1;
[7708934]896         s=s-1;
[b42ab6]897         if ( s>0 ) { # = #[1..s]; }
898      }
899   }
900   if ( s>0 )
901   {
[7708934]902      l = list(id)+#;
903      kill id;
904      list id = l;                                    //case: id = list
905      ty = "list";
[b42ab6]906   }
907   else
[65546eb]908   {
[7708934]909      ty = typeof(id);
[b42ab6]910   }
[7708934]911//------------------------------ reduce to 3 cases ---------------------------
[b42ab6]912   if( ty=="poly" or ty=="ideal" or ty=="vector"
913       or ty=="module" or ty=="matrix" )
[7708934]914   {                                                 //case: id = ideal
[3d124a7]915      ideal i = ideal(matrix(id));
[b42ab6]916      kill id;
[7708934]917      ideal id = simplify(i,2);                      //delete 0 entries
[3d124a7]918   }
[b42ab6]919   if( ty=="int" or ty=="intvec" or ty=="intmat" )
[15e59a]920   {                                                 //case: id = intvec
[b42ab6]921      if ( ty == "int" ) { intmat S =id; }
[3d124a7]922      else { intmat S = intmat(id); }
923      intvec i = S[1..nrows(S),1..ncols(S)];
[b42ab6]924      kill id;
[15e59a]925      intvec id = i;
[3d124a7]926   }
[7708934]927//------------------- consider intvec v and empty sum  -----------------------
[65546eb]928   if( tt!=0 )
[7708934]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
[65546eb]933         {
[7708934]934            return(0);                               //empty sum is 0
[65546eb]935         }
[7708934]936      }
937      id = id[v];                                    //consider part of id
938   }                                                 //corresponding to v
939
940//-------------------------- finally, add objects  ---------------------------
[65546eb]941   n = size(id);
[15e59a]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   }
[b42ab6]952}
[3d124a7]953example
954{  "EXAMPLE:"; echo = 2;
[15e59a]955   ring r1 = 0,(x,y,z),dp;
[3d124a7]956   vector pv = [xy,xz,yz,x2,y2,z2];
957   sum(pv);
[c860e9]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);
[15e59a]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));
[3d124a7]968}
969///////////////////////////////////////////////////////////////////////////////
[6f2edc]970
[15e59a]971
972///////////////////////////////////////////////////////////////////////////////
973
[6f2edc]974proc which (command)
[d2b2a7]975"USAGE:    which(command); command = string expression
[6f2edc]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
[d2b2a7]980"
[6f2edc]981{
982   int rs;
983   int i;
[a70441f]984   string fn = "which_" + string(system("pid"));
[6f2edc]985   string pn;
[a70441f]986   string cmd;
[82716e]987   if( typeof(command) != "string")
[6f2edc]988   {
[82716e]989     return (pn);
[6f2edc]990   }
[a70441f]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);
[6f2edc]1001   pn = read(fn);
[a70441f]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
[6f2edc]1015   {
[a70441f]1016     rs = 0;
[6f2edc]1017   }
1018   i = system("sh", "rm " + fn);
1019   if (rs == 0) {return (pn);}
[82716e]1020   else
[6f2edc]1021   {
1022     print (command + " not found ");
1023     return ("");
1024   }
1025}
1026example
1027{  "EXAMPLE:"; echo = 2;
[a70441f]1028    which("sh");
[6f2edc]1029}
1030///////////////////////////////////////////////////////////////////////////////
[ebbe4a]1031
1032proc watchdog(int i, string cmd)
[a7a00b]1033"USAGE:   watchdog(i,cmd); i integer, cmd string
[b42ab6]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.
[65546eb]1038NOTE:    * the MP package must be enabled
1039         * the current basering should not be watchdog_rneu, since
[b42ab6]1040           watchdog_rneu will be killed
[ebbe4a]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
[a6606e6]1044         * a second Singular process is started by this procedure
[ebbe4a]1045EXAMPLE: example watchdog; shows an example
1046"
1047{
1048  string rname=nameof(basering);
[17ee4f]1049  def rsave=basering;
[ebbe4a]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;
[65546eb]1061// fork, get the pid of the child and send it the command
[ebbe4a]1062      link l_fork="MPtcp:fork";
1063      open(l_fork);
1064      write(l_fork,quote(system("pid")));
1065      int pid=read(l_fork);
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      {
1072        if (status(l_fork, "read", "ready", j)) {break;}
1073        j = j + j;
1074      }
1075
1076// sleep in intervals of one second
1077      j = 1;
1078      if (!status(l_fork,"read","ready"))
1079      {
1080        while (j < i)
1081        {
1082          if (status(l_fork, "read", "ready", k)) {break;}
1083          j = j + 1;
1084        }
1085      }
1086// check, whether we have a result, and return it
1087      if (status(l_fork, "read", "ready"))
1088      {
1089        def result = read(l_fork);
1090        if (nameof(basering)!=rname)
1091        {
1092          def watchdog_rneu=basering;
[17ee4f]1093          setring rsave;
[7f3ad4]1094          if (!defined(result))
1095          {
1096            def result=fetch(watchdog_rneu,result);
1097          }
[ebbe4a]1098        }
1099        if(defined(watchdog_interrupt))
1100        {
[3b77465]1101          kill watchdog_interrupt;
[ebbe4a]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    }
[50cbdc]1122  }
1123  else
1124  {
[ebbe4a]1125    ERROR("MP-support is not enabled in this version of Singular.");
[65546eb]1126  }
[ebbe4a]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)
[803c5a1]1138"USAGE:   deleteSublist(v,l); intvec v; list l
[ebbe4a]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}
[8b87364]1171
1172///////////////////////////////////////////////////////////////////////////////
1173proc primecoeffs(J, list #)
[a7a00b]1174"USAGE:   primecoeffs(J[,p]); J any type which can be converted to a matrix
[8b87364]1175         e.g. ideal, matrix, vector, module, int, intvec
[a7a00b]1176         p = integer
[c94e60f]1177COMPUTE: primefactors <= p of coeffs of J (default p = 32003)
[a7a00b]1178RETURN:  a list, say l, of two intvectors:@*
1179         l[1] : the different primefactors of all coefficients of J@*
[8b87364]1180         l[2] : the different remaining factors
1181EXAMPLE: example primecoeffs; shows an example
1182"
1183{
1184   int q,ii,n,mark;;
[298d0a]1185   if (size(#) == 0)
[8b87364]1186   {
[298d0a]1187      q=32003;
[8b87364]1188   }
[298d0a]1189   else
[8b87364]1190   {
1191     if( typeof(#[1]) != "int")
1192     {
1193       ERROR("2nd parameter must be of type int"+newline);
1194     }
1195     q=#[1];
[be8dbb]1196     if (q > 32003) { q = 32003; }
[8b87364]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));
[298d0a]1206   matrix Coef=coef(I[1],p);
[8b87364]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);
[298d0a]1216   for (ii=1; ii<=size(id); ii++)
1217   {
1218     l = primefactors(number(id[ii]),q);
[c94e60f]1219     list jdl=l[2];
1220     jd = jd,jdl[1..size(jdl)];
1221     kill jdl;
1222     rest=rest, l[1];
[298d0a]1223   }
[8b87364]1224   jd = simplify(jd,6);
[298d0a]1225   for (ii=1; ii<=size(jd); ii++)
1226   {
[8b87364]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   {
[298d0a]1242      result = v,rest;
[8b87364]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);
[298d0a]1252}
[8b87364]1253///////////////////////////////////////////////////////////////////////////////
[90d772]1254proc timeFactorize(poly i,list #)
[a7a00b]1255"USAGE:  timeFactorize(p,d); poly p , integer d
[3c4dcc]1256RETURN:  factorize(p) if the factorization finished after d-1
[90d772]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);
1272        link l_fork = "MPtcp:fork";
1273        open(l_fork);
1274        write(l_fork, quote(system("pid")));
1275        int pid = read(l_fork);
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          {
1283            if (status(l_fork, "read", "ready", j)) {break;}
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        {
1292          if (status(l_fork, "read", "ready", 1000000)) {break;}
1293          j = j + 1;
1294        }
1295
1296        if (status(l_fork, "read", "ready"))
1297        {
1298          def result = read(l_fork);
1299          if (bs != nameof(basering))
1300          {
1301            def PP = basering;
1302            setring P;
1303            def result = imap(PP, result);
1304            kill PP;
1305          }
[3b77465]1306          kill l_fork;
[90d772]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);
1349        link l_fork = "MPtcp:fork";
1350        open(l_fork);
1351        write(l_fork, quote(system("pid")));
1352        int pid = read(l_fork);
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          {
1360            if (status(l_fork, "read", "ready", j)) {break;}
1361            j = j + j;
1362          }
1363        }
1364        j = 1;
1365        while (j < wait)
1366        {
1367          if (status(l_fork, "read", "ready", 1000000)) {break;}
1368          j = j + 1;
1369        }
1370        if (status(l_fork, "read", "ready"))
1371        {
1372          def result = read(l_fork);
1373          if (bs != nameof(basering))
1374          {
1375            def PP = basering;
1376            setring P;
1377            def result = imap(PP, result);
1378            kill PP;
1379          }
[3b77465]1380          kill l_fork;
[90d772]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)
[a7a00b]1410NOTE:    changes variables to make the last variable the principal
[90d772]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.