# source:git/Singular/LIB/general.lib@e67c6f

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