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

spielwiese
Last change on this file since 7708934 was 7708934, checked in by Gert-Martin Greuel <greuel@…>, 23 years ago
* GMG Prozeduren sum und product neu ueberarbeitet, leere Summe und leeres * Produkt mit 0 bzw 1 zurueckgegeben. Funktionalitaet erweitert. git-svn-id: file:///usr/local/Singular/svn/trunk@4994 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 37.7 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.33 2000-12-30 03:00:57 greuel 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 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]
25 watchdog(i,cmd);       only wait for result of command cmd for i seconds
26 which(command);        search for command and return absolute path, if found
27           (parameters in square brackets [] are optional)
28";
29
30LIB "inout.lib";
31///////////////////////////////////////////////////////////////////////////////
32
33proc A_Z (string s,int n)
34"USAGE:   A_Z(\"a\",n);  a any letter, n integer (-26<= n <=26, !=0)
35RETURN:  string of n small (if a is small) or capital (if a is capital)
36         letters, comma separated, beginning with a, in alphabetical
37         order (or revers alphabetical order if n<0)
38EXAMPLE: example A_Z; shows an example
39"
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;
84   execute(sR);
85   R;
86}
87///////////////////////////////////////////////////////////////////////////////
88proc ASCII (list #)
89"USAGE:   ASCII([n,m]); n,m= integers (32 <= n <= m <= 126)
90RETURN:   string of printable ASCII characters (no native language support)
91          ASCII():    string of  all ASCII characters with its numbers,
92          ASCII(n):   n-th ASCII character
93          ASCII(n,m): n-th up to m-th ASCII character (inclusive)
94EXAMPLE: example ASCII; shows an example
95"
96{
97  string s1 =
98 "     !    \"    #    $    %    &    '    (    )    *    +    ,    -    .
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///////////////////////////////////////////////////////////////////////////////
136
137proc binomial (int n, int k, list #)
138"USAGE:   binomial(n,k[,p]); n,k,p integers
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
145NOTE:    In any characteristic, binomial(n,k) = coefficient of x^k in (1+x)^n
146SEE ALSO: prime
147EXAMPLE: example binomial; shows an example
148"
149{
150   int str,p;
151//---------------------------- initialization -------------------------------
152   if ( size(#) == 0 )
153   {  str = 1;
154      ring bin = 0,x,dp;
155      number r=1;
156   }
157   if ( size(#) > 0 )
158   {
159      p = (#[1]!=0)*prime(#[1]);
160      if ( defined(basering) )
161      {
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;
175      }
176   }
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 }
191example
192{ "EXAMPLE:"; echo = 2;
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);
260}
261///////////////////////////////////////////////////////////////////////////////
262
263proc factorial (int n, list #)
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
272"
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++)
302   {
303      r=r*l;
304   }
305   if ( str==1 ) { return(string(r)); }
306   else { return(r); }
307}
308example
309{ "EXAMPLE:"; echo = 2;
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;
314}
315///////////////////////////////////////////////////////////////////////////////
316
317proc fibonacci (int n, list #)
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
326EXAMPLE: example fibonacci; shows an example
327"
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 --------------------------------
356   for (ii=3; ii<=n; ii=ii+1)
357   {
358      h = f+g; f = g; g = h;
359    }
360   if ( str==1 ) { return(string(h)); }
361   else { return(h); }
362}
363example
364{ "EXAMPLE:"; echo = 2;
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
368   b;
369}
370///////////////////////////////////////////////////////////////////////////////
371
372proc kmemory (list #)
373"USAGE:   kmemory([n,[v]]); n,v integers
374RETURN:  memory in kilobyte of type int
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
379DISPLAY: detailed information about allocated and used memory if v!=0
380NOTE:    kmemory uses internal function 'memory' to compute kilobyte, and
381         is the same as 'memory' for n!=0,1,2
382EXAMPLE: example kmemory; shows an example
383"
384{
385   int n;
386   int verb;
387   if (size(#) != 0)
388   {
389     n=#[1];
390     if (size(#) >1)
391     { verb=#[2]; }
392   }
393
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   }
403   return ((memory(n)+1023)/1024);
404}
405example
406{ "EXAMPLE:"; echo = 2;
407   kmemory();
408   kmemory(1,1);
409}
410///////////////////////////////////////////////////////////////////////////////
411
412proc killall
413"USAGE:   killall(); (no parameter)
414         killall(\"type_name\");
415         killall(\"not\", \"type_name\");
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
425NOTE:    killall should never be used inside a procedure
426EXAMPLE: example killall; shows an example AND KILLS ALL YOUR VARIABLES
427"
428{
429  list @marie=names();
430  int no_kill, @joni;
431  for ( @joni=1; @joni<=size(#);  @joni++)
432  {
433    if (typeof(#[@joni]) != "string")
434    {
435      ERROR("Need string as " + string(@joni) + "th argument");
436    }
437  }
438 
439  // kills all user-defined variables but not loaded procedures
440  if( size(#)==0 )
441  {
442    for ( @joni=size(@marie); @joni>0; @joni-- )
443    {
444      if( @marie[@joni]!="LIB" and typeof(`@marie[@joni]`)!="proc" )
445      { kill `@marie[@joni]`; }
446    }
447  }
448  else
449  {
450    // kills all user-defined variables
451    if( size(#)==1 )
452    {
453      // of type proc
454      if( #[1] == "proc" )
455      {
456        for ( @joni=size(@marie); @joni>0; @joni-- )
457        {
458          if((@marie[@joni]!="killall") and (@marie[@joni]=="LIB" or
459                                       typeof(`@marie[@joni]`)=="proc"))
460          { kill `@marie[@joni]`; }
461        }
462      }
463      else
464      { 
465        // other types
466        for ( @joni=size(@marie); @joni>2; @joni-- )
467        {
468          if(typeof(`@marie[@joni]`)==#[1] and @marie[@joni]!="LIB"
469             and typeof(`@marie[@joni]`)!="proc")
470          { kill `@marie[@joni]`; }
471        }
472      }
473    }
474    else
475    {
476      // kills all user-defined variables whose name or type is not #i
477      for ( @joni=size(@marie); @joni>2; @joni-- )
478      {
479        if ( @marie[@joni] != "LIB" && typeof(`@marie[@joni]`) != "proc")
480        {
481          no_kill = 0;
482          for (j=2; j<= size(#); j++)
483          {
484            if (typeof(`@marie[@joni]`)==#[j] or @marie[@joni] == #[j])
485            {
486              no_kill = 1;
487              break;
488            }
489          }
490          if (! no_kill)
491          {
492            kill `@marie[@joni]`;
493          }
494        }
495      }
496    }
497  }
498}
499example
500{ "EXAMPLE:"; echo = 2;
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();
510}
511///////////////////////////////////////////////////////////////////////////////
512
513proc number_e (int n)
514"USAGE:   number_e(n);  n integer
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
520EXAMPLE: example number_e; shows an example
521"
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];
537         u[m] = s div (m+1);
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   }
543   if( t==1 )
544   { dbprint(printlevel-voice+2,"// "+result[1,n+1]);
545     return(r);
546   }
547   return(result[1,n+1]);
548}
549example
550{ "EXAMPLE:"; echo = 2;
551   number_e(30);"";
552   ring R = 0,t,lp;
553   number e = number_e(30);
554   e;
555}
556///////////////////////////////////////////////////////////////////////////////
557
558proc number_pi (int n)
559"USAGE:   number_pi(n);  n positive integer
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
565EXAMPLE: example number_pi; shows an example
566"
567{
568   int i,m,t,e,q,N;
569   intvec r,p,B,Prelim;
570   string result,prelim;
571   N = (10*n) div 3 + 2;
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];
585         q    = e div B[m];
586         e    = q*(m-1)+p[m-1];
587      }
588      r[1] = e%10;
589      q    = e div 10;
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      {
603         Prelim = (Prelim+1)-((Prelim+1) div 10)*10;
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];
617   if( t==1 )
618   { dbprint(printlevel-voice+2,"// "+result);
619     return(pi);
620   }
621   return(result);
622}
623example
624{ "EXAMPLE:"; echo = 2;
625   number_pi(11);"";
626   ring r = (real,10),t,dp;
627   number pi = number_pi(11); pi;
628}
629///////////////////////////////////////////////////////////////////////////////
630
631proc primes (int n, int m)
632"USAGE:   primes(n,m);  n,m integers
633RETURN:  intvec, consisting of all primes p, prime(n)<=p<=m, in increasing
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
637EXAMPLE: example primes; shows an example
638"
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;
648    primes(50,100);"";
649    intvec v = primes(37,1); v;
650}
651///////////////////////////////////////////////////////////////////////////////
652
653proc product (id, list #)
654"USAGE:    product(id[,v]); id ideal/vector/module/matrix/intvec/intmat/list,
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).
662@*        If v is outside the range of id, we have the empty product and the
663          result will be 1 (of type int).
664EXAMPLE:  example product; shows an example
665"
666
667//-------------------- initialization and special feature ---------------------
668   int n,j,tt;
669   string ty;                                         //will become type of id
670   list l;
671
672// We wish to allow something like product(x(1..10)) if x(1),...,x(10) are
673// variables. x(1..10) is a list of polys and enters the procedure with
674// id=x(1) and # a list with 9 polys, #[1]= x(2),...,#[9]= x(10). Hence, in
675// this case # is never empty. If an additional intvec v is given,
676// it is added to #, so we have to separate it first and make
677// the rest a list which has to be multiplied.
678
679   int s = size(#);
680   if( s!=0 )
681   {  if ( typeof(#[s])=="intvec" or typeof(#[s])=="int")   
682      { 
683         intvec v = #[s];
684         tt=1;
685         s=s-1;
686         if ( s>0 ) { # = #[1..s]; }
687      }
688   }
689   if ( s>0 )
690   {
691      l = list(id)+#;
692      kill id;
693      list id = l;                                    //case: id = list
694      ty = "list";
695      n = size(id);
696   }
697   else
698   {
699      ty = typeof(id);
700      if( ty == "list" )
701      { n = size(id); }
702   }
703//------------------------------ reduce to 3 cases ---------------------------
704   if( ty=="poly" or ty=="ideal" or ty=="vector"
705       or ty=="module" or ty=="matrix" )
706   {
707      ideal i = ideal(matrix(id));
708      kill id;
709      ideal id = i;                                   //case: id = ideal
710      n = ncols(id);
711   }
712   if( ty=="int" or ty=="intvec" or ty=="intmat" )
713   {
714      if ( ty == "int" ) { intmat S =id; }
715      else { intmat S = intmat(id); }
716      intvec i = S[1..nrows(S),1..ncols(S)];
717      kill id;
718      intvec id = i;                                  //case: id = intvec
719      n = size(id);
720   }
721//--------------- consider intvec v and empty product  -----------------------
722   if( tt!=0 )
723   {
724      for (j=1; j<=size(v); j++)
725      {
726         if ( v[j] <= 0 or v[j] > n )                 //v outside range of id
727         {
728            return(1);                                //empty product is 1
729         }                               
730      }
731      id = id[v];                                     //consider part of id
732   }                                                  //corresponding to v
733//--------------------- special case: one factor is zero  ---------------------
734   if ( typeof(id) == "ideal")
735   {
736      if( size(id) < ncols(id) )
737      {
738          poly f; return(f);
739      }
740   }
741//-------------------------- finally, multiply objects  -----------------------
742   n = size(id);
743   def f(1) = id[1];
744   for( j=2; j<=n; j=j+1 ) { def f(j)=f(j-1)*id[j]; }
745   return(f(n));
746}
747example
748{  "EXAMPLE:"; echo = 2;
749   ring r= 0,(x,y,z),dp;
750   ideal m = maxideal(1);
751   product(m);
752   product(m[2..3]);
753   matrix M[2][3] = 1,x,2,y,3,z;
754   product(M);
755   intvec v=2,4,6;
756   product(M,v);
757   intvec iv = 1,2,3,4,5,6,7,8,9;
758   v=1..5,7,9;
759   product(iv,v);
760   intmat A[2][3] = 1,1,1,2,2,2;
761   product(A,3..5);
762}
763///////////////////////////////////////////////////////////////////////////////
764proc ringweights (list # )
765"USAGE:   ringweights(P); P=name of an existing ring (true name, not a string)
766RETURN:  intvec consisting of the weights of the variables of P, as they
767         appear when typing P;.
768NOTE:    This is useful when enlarging P but keeping the weights of the old
769         variables.
770EXAMPLE: example ringweights; shows an example
771"
772{
773   int ii,q,fi,fo,fia;
774   intvec rw,nw;
775   string os;
776   def P = #[1];
777   string osP = ordstr(P);
778   fo  = 1;
779//------------------------- find weights in ordstr(P) -------------------------
780   fi  = find(osP,"(",fo);
781   fia = find(osP,"a",fo)+find(osP,"w",fo)+find(osP,"W",fo);
782   while ( fia )
783   {
784      os = osP[fi+1,find(osP,")",fi)-fi-1];
785      if( find(os,",") )
786      {
787         execute("nw = "+os+";");
788         if( size(nw) > ii )
789         {
790             rw = rw,nw[ii+1..size(nw)];
791         }
792         else  {  ii = ii - size(nw); }
793
794         if( find(osP[1,fi],"a",fo) ) { ii = size(nw); }
795      }
796      else
797      {
798         execute("q = "+os+";");
799         if( q > ii )
800         {
801            nw = 0; nw[q-ii] = 0;
802            nw = nw + 1;          //creates an intvec 1,...,1 of length q-ii
803            rw = rw,nw;
804         }
805         else { ii = ii - q; }
806      }
807      fo  = fi+1;
808      fi  = find(osP,"(",fo);
809      fia = find(osP,"a",fo)+find(osP,"w",fo)+find(osP,"W",fo);
810   }
811//-------------- adjust weight vector to length = nvars(P)  -------------------
812   if( fo > 1 )
813   {                                            // case when weights were found
814      rw = rw[2..size(rw)];
815      if( size(rw) > nvars(P) )
816      {
817         rw = rw[1..nvars(P)];
818      }
819      if( size(rw) < nvars(P) )
820      {
821         nw=0; nw[nvars(P)-size(rw)]=0; nw=nw+1; rw=rw,nw;
822      }
823   }
824   else
825   {                                         // case when no weights were found
826      rw[nvars(P)]= 0; rw=rw+1;
827   }
828   return(rw);
829}
830example
831{"EXAMPLE:";  echo = 2;
832  ring r0 = 0,(x,y,z),dp;
833  ringweights(r0);
834  ring r1 = 0,x(1..5),(ds(3),wp(2,3));
835  ringweights(r1);"";
836  // an example for enlarging the ring, keeping the first weights:
837  intvec v = ringweights(r1),6,2,3,4,5;
838  ring R = 0,x(1..10),(a(v),dp);
839  ordstr(R);
840}
841
842///////////////////////////////////////////////////////////////////////////////
843
844proc sort (id, list #)
845"USAGE:   sort(id[v,o,n]); id=ideal/module/intvec/list(of intvec's or int's)
846@*       sort may be called with 1, 2 or 3 arguments in the following way:
847@*       sort(id[v,n]); v=intvec of positive integers, n=integer,
848@*       sort(id[o,n]); o=string (any allowed ordstr of a ring), n=integer
849RETURN:  a list l of two elements:
850@format
851        l[1]: object of same type as input but sorted in the following way:
852           - if id=ideal/module: generators of id are sorted w.r.t. intvec v
853             (id[v[1]] becomes 1-st, id[v[2]]  2-nd element, etc.). If no v is
854             present, id is sorted w.r.t. ordering o (if o is given) or w.r.t.
855             actual monomial ordering (if no o is given):
856             NOTE: generators with SMALLER(!) leading term come FIRST
857             (e.g. sort(id); sorts backwards to actual monomial ordering)
858           - if id=list of intvec's or int's: consider a list element, say
859             id[1]=3,2,5, as exponent vector of the monomial x^3*y^2*z^5;
860             the corresponding monomials are ordered w.r.t. intvec v (s.a.).
861             If no v is present, the monomials are sorted w.r.t. ordering o
862             (if o is given) or w.r.t. lexicographical ordering (if no o is
863             given). The corresponding ordered list of exponent vectors is
864             returned.
865             (e.g. sort(id); sorts lexicographically, smaller int's come first)
866             WARNING: Since negative exponents create the 0 polynomial in
867             Singular, id should not contain negative integers: the result
868             might not be as expected
869           - if id=intvec: id is treated as list of integers
870           - if n!=0 the ordering is inverse, i.e. w.r.t. v(size(v)..1)
871             default: n=0
872         l[2]: intvec, describing the permutation of the input (hence l[2]=v
873             if v is given (with positive integers))
874@end format
875NOTE:    If v is given id may be any simply indexed object (e.g. any list or
876         string); if v[i]<0 and i<=size(id) v[i] is set internally to i;
877         entries of v must be pairwise distinct to get a permutation if id.
878         Zero generators of ideal/module are deleted
879EXAMPLE: example sort; shows an example
880"
881{  int ii,jj,s,n = 0,0,1,0;
882   intvec v;
883   if ( defined(basering) ) { def P = basering; }
884   if ( size(#)==0 and (typeof(id)=="ideal" or typeof(id)=="module") )
885   {
886      id = simplify(id,2);
887      for ( ii=1; ii<size(id); ii++ )
888      {
889         if ( id[ii]!=id[ii+1] ) { break;}
890      }
891      if ( ii != size(id) ) { v = sortvec(id); }
892      else  { v = size(id)..1; }
893   }
894   if ( size(#)>=1 and (typeof(id)=="ideal" or typeof(id)=="module") )
895   {
896      if ( typeof(#[1])=="string" )
897      {
898         execute("ring r1 =("+charstr(P)+"),("+varstr(P)+"),("+#[1]+");");
899         def i = imap(P,id);
900         v = sortvec(i);
901         setring P;
902         n=2;
903      }
904   }
905   if ( typeof(id)=="intvec" or typeof(id)=="list" and n==0 )
906   {
907      string o;
908      if ( size(#)==0 ) { o = "lp"; n=1; }
909      if ( size(#)>=1 )
910      {
911         if ( typeof(#[1])=="string" ) { o = #[1]; n=1; }
912      }
913   }
914   if ( typeof(id)=="intvec" or typeof(id)=="list" and n==1 )
915   {
916      if ( typeof(id)=="list" )
917      {
918         for (ii=1; ii<=size(id); ii++)
919         {
920            if (typeof(id[ii]) != "intvec" and typeof(id[ii]) != "int")
921               { "// list elements must be intvec/int"; return(); }
922            else
923               { s=size(id[ii])*(s < size(id[ii])) + s*(s >= size(id[ii])); }
924         }
925      }
926      execute("ring r=0,x(1..s),("+o+");");
927      ideal i;
928      poly f;
929      for (ii=1; ii<=size(id); ii++)
930      {
931         f=1;
932         for (jj=1; jj<=size(id[ii]); jj++)
933         {
934            f=f*x(jj)^(id[ii])[jj];
935         }
936         i[ii]=f;
937      }
938      v = sort(i)[2];
939   }
940   if ( size(#)!=0 and n==0 ) { v = #[1]; }
941   if( size(#)==2 )
942   {
943      if ( #[2] != 0 ) { v = v[size(v)..1]; }
944   }
945   s = size(v);
946   if( size(id) < s ) { s = size(id); }
947   def m = id;
948   if ( size(m) != 0 )
949   {
950      for ( jj=1; jj<=s; jj=jj+1)
951      {
952         if ( v[jj]<=0 ) { v[jj]=jj; }
953         m[jj] = id[v[jj]];
954      }
955   }
956   if ( v == 0 ) { v = 1; }
957   list L=m,v;
958   return(L);
959}
960example
961{  "EXAMPLE:"; echo = 2;
962   ring r0 = 0,(x,y,z,t),lp;
963   ideal i = x3,z3,xyz;
964   sort(i);            //sorts using lex ordering, smaller polys come first
965 
966   sort(i,3..1);
967
968   sort(i,"ls")[1];     //sort w.r.t. negative lex ordering
969
970   intvec v =1,10..5,2..4;v;
971   sort(v)[1];          // sort v lexicographically
972
973   sort(v,"Dp",1)[1];   // sort v w.r.t (total sum, reverse lex)
974}
975///////////////////////////////////////////////////////////////////////////////
976proc sum (id, list #)
977"USAGE:    sum(id[,v]); id ideal/vector/module/matrix/intvec/intmat/list,
978          v intvec  (default: v=1..number of entries of id)
979ASSUME:   list members can be added.
980RETURN:   The sum of all entries of id [with index given by v] of type
981          depending on the entries of id.
982NOTE:     If id is not a list, id is treated as a list of polys resp. integers.
983          A module m is identified with the corresponding matrix M (columns
984          of M generate m).
985@*        If v is outside the range of id, we have the empty sum and the
986          result will be 0 (of type int).
987EXAMPLE:  example sum; shows an example
988"
989{
990//-------------------- initialization and special feature ---------------------
991   int n,j,tt;
992   string ty;                                  // will become type of id
993   list l;
994
995// We wish to allow something like sum(x(1..10)) if x(1),...,x(10) are
996// variables. x(1..10) is a list of polys and enters the procedure with
997// id=x(1) and # a list with 9 polys, #[1]= x(2),...,#[9]= x(10). Hence, in
998// this case # is never empty. If an additional intvec v is given,
999// it is added to #, so we have to separate it first and make
1000// the rest a list which has to be added.
1001
1002   int s = size(#);
1003   if( s!=0 )
1004   {  if ( typeof(#[s])=="intvec" or typeof(#[s])=="int")
1005      {  intvec v = #[s];
1006         tt=1;
1007         s=s-1;
1008         if ( s>0 ) { # = #[1..s]; }
1009      }
1010   }
1011   if ( s>0 )
1012   {
1013      l = list(id)+#;
1014      kill id;
1015      list id = l;                                    //case: id = list
1016      ty = "list";
1017   }
1018   else
1019   {
1020      ty = typeof(id);
1021   }
1022//------------------------------ reduce to 3 cases ---------------------------
1023   if( ty=="poly" or ty=="ideal" or ty=="vector"
1024       or ty=="module" or ty=="matrix" )
1025   {                                                 //case: id = ideal
1026      ideal i = ideal(matrix(id));
1027      kill id;
1028      ideal id = simplify(i,2);                      //delete 0 entries
1029   }
1030   if( ty=="int" or ty=="intvec" or ty=="intmat" )
1031   {
1032      if ( ty == "int" ) { intmat S =id; }
1033      else { intmat S = intmat(id); }
1034      intvec i = S[1..nrows(S),1..ncols(S)];
1035      kill id;
1036      intvec id = i;                                 //case: id = intvec
1037   }
1038//------------------- consider intvec v and empty sum  -----------------------
1039   if( tt!=0 )
1040   {
1041      for (j=1; j<=size(v); j++)
1042      {
1043         if ( v[j] <= 0 or v[j] > size(id) )         //v outside range of id
1044         {
1045            return(0);                               //empty sum is 0
1046         }                               
1047      }
1048      id = id[v];                                    //consider part of id
1049   }                                                 //corresponding to v
1050
1051//-------------------------- finally, add objects  ---------------------------
1052   n = size(id);
1053   def f(1) = id[1];
1054   for( j=2; j<=n; j=j+1 ) { def f(j)=f(j-1)+id[j]; }
1055   return(f(n));   int n,j,tt;
1056}
1057example
1058{  "EXAMPLE:"; echo = 2;
1059   ring r= 0,(x,y,z),dp;
1060   vector pv = [xy,xz,yz,x2,y2,z2];
1061   sum(pv);
1062   sum(pv,2..5);
1063   matrix M[2][3] = 1,x,2,y,3,z;
1064   intvec w=2,4,6;
1065   sum(M,w);
1066   intvec iv = 1,2,3,4,5,6,7,8,9;
1067   sum(iv,2..4);
1068}
1069///////////////////////////////////////////////////////////////////////////////
1070
1071proc which (command)
1072"USAGE:    which(command); command = string expression
1073RETURN:   Absolute pathname of command, if found in search path.
1074          Empty string, otherwise.
1075NOTE:     Based on the Unix command 'which'.
1076EXAMPLE:  example which; shows an example
1077"
1078{
1079   int rs;
1080   int i;
1081   string fn = "which_" + string(system("pid"));
1082   string pn;
1083   string cmd;
1084   if( typeof(command) != "string")
1085   {
1086     return (pn);
1087   }
1088   if (system("uname") != "ix86-Win")
1089   {
1090     cmd = "which ";
1091   }
1092   else
1093   {
1094     // unfortunately, it does not take -path
1095     cmd = "type ";
1096   }
1097   i = system("sh", cmd + command + " > " + fn);
1098   pn = read(fn);
1099   if (system("uname") != "ix86-Win")
1100   {
1101     // TBC: Hmm... should parse output to get rid of 'command is '
1102     pn[size(pn)] = "";
1103     i = 1;
1104     while ((pn[i] != " ") and (pn[i] != ""))
1105     {
1106       i = i+1;
1107     }
1108     if (pn[i] == " ") {pn[i] = "";}
1109     rs = system("sh", "ls " + pn + " > " + fn + " 2>&1 ");
1110   }
1111   else
1112   {
1113     rs = 0;
1114   }
1115   i = system("sh", "rm " + fn);
1116   if (rs == 0) {return (pn);}
1117   else
1118   {
1119     print (command + " not found ");
1120     return ("");
1121   }
1122}
1123example
1124{  "EXAMPLE:"; echo = 2;
1125    which("sh");
1126}
1127///////////////////////////////////////////////////////////////////////////////
1128
1129proc watchdog(int i, string cmd)
1130"USAGE:   watchdog(i,cmd); i integer; cmd string
1131RETURN:  Result of cmd, if the result can be computed in i seconds.
1132         Otherwise the computation is interrupted after i seconds,
1133         the string "Killed" is returned and the global variable
1134         'watchdog_interrupt' is defined.
1135NOTE:    * the MP package must be enabled
1136         * the current basering should not be watchdog_rneu, since
1137           watchdog_rneu will be killed
1138         * if there are variable names of the structure x(i) all
1139           polynomials have to be put into eval(...) in order to be
1140           interpreted correctly
1141         * a second Singular process is started by this procedure
1142EXAMPLE: example watchdog; shows an example
1143"
1144{
1145  string rname=nameof(basering);
1146  if (defined(watchdog_rneu))
1147  {
1148    kill watchdog_rneu;
1149  }
1150// If we do not have MP-links, watchdog cannot be used
1151  if (system("with","MP"))
1152  {
1153    if ( i > 0 )
1154    {
1155      int j=10;
1156      int k=999999;
1157// fork, get the pid of the child and send it the command   
1158      link l_fork="MPtcp:fork";
1159      open(l_fork);
1160      write(l_fork,quote(system("pid")));
1161      int pid=read(l_fork);
1162      execute("write(l_fork,quote(" + cmd + "));");
1163
1164
1165// sleep in small, but growing intervals for appr. 1 second
1166      while(j < k)
1167      {
1168        if (status(l_fork, "read", "ready", j)) {break;}
1169        j = j + j;
1170      }
1171
1172// sleep in intervals of one second
1173      j = 1;
1174      if (!status(l_fork,"read","ready"))
1175      {
1176        while (j < i)
1177        {
1178          if (status(l_fork, "read", "ready", k)) {break;}
1179          j = j + 1;
1180        }
1181      }
1182// check, whether we have a result, and return it
1183      if (status(l_fork, "read", "ready"))
1184      {
1185        def result = read(l_fork);
1186        if (nameof(basering)!=rname)
1187        {
1188          def watchdog_rneu=basering;
1189        }
1190        if(defined(watchdog_interrupt))
1191        {
1192          kill (watchdog_interrupt);
1193        }
1194        close(l_fork);
1195      }
1196      else
1197      {
1198        string result="Killed";
1199        if(!defined(watchdog_interrupt))
1200        {
1201          int watchdog_interrupt=1;
1202          export watchdog_interrupt;
1203        }
1204        close(l_fork);
1205        j = system("sh","kill " + string(pid));
1206      }
1207      if (defined(watchdog_rneu))
1208      {
1209        keepring watchdog_rneu;
1210      }
1211      return(result);
1212    }
1213    else
1214    {
1215      ERROR("First argument of watchdog has to be a positive integer.");
1216    }
1217    ERROR("MP-support is not enabled in this version of Singular.");
1218  }
1219}
1220example
1221{ "EXAMPLE:"; echo=2;
1222  ring r=0,(x,y,z),dp;
1223  poly f=x^30+y^30;
1224  watchdog(1,"factorize(eval("+string(f)+"))");
1225  watchdog(100,"factorize(eval("+string(f)+"))");
1226}
1227///////////////////////////////////////////////////////////////////////////////
1228
1229proc deleteSublist(intvec v,list l)
1230"USAGE:   deleteSublist(v,l); intvec v; list l
1231         where the entries of the integer vector v correspond to the
1232         positions of the elements to be deleted
1233RETURN:  list without the deleted elements
1234EXAMPLE: example deleteSublist; shows an example"
1235{
1236  list k;
1237  int i,j,skip;
1238  j=1;
1239  skip=0;
1240  intvec vs=sort(v)[1];
1241  for ( i=1 ; i <=size(vs) ; i++)
1242  {
1243    while ((j+skip) < vs[i])
1244    {
1245      k[j] = l[j+skip];
1246      j++;
1247    }
1248    skip++;
1249  }
1250  if(vs[size(vs)]<size(l))
1251  {
1252    k=k+list(l[(vs[size(vs)]+1)..size(l)]);
1253  }
1254  return(k);
1255}
1256example
1257{ "EXAMPLE:"; echo=2;
1258   list l=1,2,3,4,5;
1259   intvec v=1,3,4;
1260   l=deleteSublist(v,l);
1261   l;
1262}
1263///////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.