# source:git/Singular/LIB/general.lib@3fb0d8

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