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

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