source: git/Singular/LIB/general.lib @ 48c165a

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