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

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