source: git/Singular/LIB/general.lib @ 49998f

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