# source:git/Singular/LIB/general.lib@1b2216

jengelh-datetimespielwiese
Last change on this file since 1b2216 was 1b2216, checked in by Hans Schoenemann <hannes@…>, 10 years ago
• Property mode set to `100644`
File size: 36.6 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
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 or intvec: sorted w.r.t. < (indep. of other arguments)
670           - if n!=0 the ordering is inverse, i.e. w.r.t. v(size(v)..1)
671             default: n=0
672         l[2]: intvec, describing the permutation of the input (hence l[2]=v
673             if v is given (with positive integers))
674@end format
675NOTE:    If v is given id may be any simply indexed object (e.g. any list or
676         string); if v[i]<0 and i<=size(id) v[i] is set internally to i;
677         entries of v must be pairwise distinct to get a permutation id.
678         Zero generators of ideal/module are deleted
679         If 'o' is given, the input is sorted by considering leading terms
680         w.r.t. the new ring ordering given by 'o'
681EXAMPLE: example sort; shows an example
682"
683{  int ii,jj,s,n = 0,0,1,0;
684   intvec v;
685   if ( defined(basering) ) { def P = basering; }
686   if ( size(#)==0 and (typeof(id)=="ideal" or typeof(id)=="module"
687                        or typeof(id)=="matrix"))
688   {
689      id = simplify(id,2);
690      for ( ii=1; ii<ncols(id); ii++ )
691      {
692         if ( id[ii]!=id[ii+1] ) { break;}
693      }
694      if ( ii != ncols(id) ) { v = sortvec(id); }
695      else  { v = ncols(id)..1; }
696   }
697   if ( size(#)>=1 and (typeof(id)=="ideal" or typeof(id)=="module"
698                        or typeof(id)=="matrix") )
699   {
700      if ( typeof(#[1])=="string" )
701      {
702         execute("ring r1 =("+charstr(P)+"),("+varstr(P)+"),("+#[1]+");");
703         def i = imap(P,id);
704         v = sortvec(i);
705         setring P;
706         n=2;
707      }
708   }
709   if ( typeof(id)=="intvec" or typeof(id)=="list" )
710   {
711      int Bn,Bi,Bj,Bb;
712      intvec pivot;
713      for(Bi=size(id);Bi>0;Bi--) { pivot[Bi]=Bi; }
714      while(Bj==0)
715      {
716        Bi++;
717        Bj=1;
718        for(Bn=1;Bn<=size(id)-Bi;Bn++)
719        {
720          if(id[pivot[Bn]]>id[pivot[Bn+1]])
721          {
722            Bb=pivot[Bn];
723            pivot[Bn]=pivot[Bn+1];
724            pivot[Bn+1]=Bb;
725            Bj=0;
726          }
727        }
728      }
729      def Br=id;
730      for(Bi=size(id);Bi>0;Bi--) { Br[Bi]=id[pivot[Bi]]; }
731      return(list(Br,pivot));
732   }
733   if ( size(#)!=0 and n==0 ) { v = #[1]; }
734   if( size(#)==2 )
735   {
736      if ( #[2] != 0 ) { v = v[size(v)..1]; }
737   }
738   s = size(v);
739   if( size(id) < s ) { s = size(id); }
740   def m = id;
741   if ( size(m) != 0 )
742   {
743      for ( jj=1; jj<=s; jj++)
744      {
745         if ( v[jj]<=0 ) { v[jj]=jj; }
746         m[jj] = id[v[jj]];
747      }
748   }
749   if ( v == 0 ) { v = 1; }
750   list L=m,v;
751   return(L);
752}
753example
754{  "EXAMPLE:"; echo = 2;
755   ring r0 = 0,(x,y,z,t),lp;
756   ideal i = x3,z3,xyz;
757   sort(i);            //sorts using lex ordering, smaller polys come first
758
759   sort(i,3..1);
760
761   sort(i,"ls")[1];     //sort w.r.t. negative lex ordering
762
763   intvec v =1,10..5,2..4;v;
764   sort(v)[1];          // sort v lexicographically
765
766   sort(v,"Dp",1)[1];   // sort v w.r.t (total sum, reverse lex)
767
769   module M = [0, 1, 1, 0], [1, 0, 0, 1]; M;
772
773   // In order to sort M wrt a NEW ordering by considering OLD leading
774   // terms use one of the following equivalent commands:
775   module( M[ sort(lead(M), "c,dp")[2] ] );
776   sort( M, sort(lead(M), "c,dp")[2] )[1];
777
778   // BUG: Please, don't use this sort for integer vectors or lists
779   // with them if there can be negative integers!
780   // TODO: for some HiWi
781   sort(3..-3)[1];
782   sort(list(-v, v))[1];
783
784}
785///////////////////////////////////////////////////////////////////////////////
786
787static proc lsum (int n, list l)
788{ if (n>10)
789  { return( lsum(n div 2,list(l[1..(n div 2)])) + lsum(n-n div 2, list(l[(n div 2+1)..n])) );
790  }
791  else
792  { def Summe=l[1];
793    for (int i=2;i<=n;i++)
794    { Summe=Summe+l[i];
795    }
796    return(Summe);
797  }
798}
799
800///////////////////////////////////////////////////////////////////////////////
801
802proc sum (id, list #)
803"USAGE:    sum(id[,v]); id ideal/vector/module/matrix/intvec/intmat/list,
804          v intvec  (default: v=1..number of entries of id)
805ASSUME:   list members can be added.
806RETURN:   The sum of all entries of id [with index given by v] of type
807          depending on the entries of id.
808NOTE:     If id is not a list, id is treated as a list of polys resp. integers.
809          A module m is identified with the corresponding matrix M (columns
810          of M generate m).
811@*        If v is outside the range of id, we have the empty sum and the
812          result will be 0 (of type int).
813EXAMPLE:  example sum; shows an example
814"
815{
816//-------------------- initialization and special feature ---------------------
817   int n,j,tt;
818   string ty;                                  // will become type of id
819   list l;
820
821// We wish to allow something like sum(x(1..10)) if x(1),...,x(10) are
822// variables. x(1..10) is a list of polys and enters the procedure with
823// id=x(1) and # a list with 9 polys, #[1]= x(2),...,#[9]= x(10). Hence, in
824// this case # is never empty. If an additional intvec v is given,
825// it is added to #, so we have to separate it first and make
826// the rest a list which has to be added.
827
828   int s = size(#);
829   if( s!=0 )
830   {  if ( typeof(#[s])=="intvec" or typeof(#[s])=="int")
831      {  intvec v = #[s];
832         tt=1;
833         s=s-1;
834         if ( s>0 ) { # = #[1..s]; }
835      }
836   }
837   if ( s>0 )
838   {
839      l = list(id)+#;
840      kill id;
841      list id = l;                                    //case: id = list
842      ty = "list";
843   }
844   else
845   {
846      ty = typeof(id);
847   }
848//------------------------------ reduce to 3 cases ---------------------------
849   if( ty=="poly" or ty=="ideal" or ty=="vector"
850       or ty=="module" or ty=="matrix" )
851   {                                                 //case: id = ideal
852      ideal i = ideal(matrix(id));
853      kill id;
854      ideal id = simplify(i,2);                      //delete 0 entries
855   }
856   if( ty=="int" or ty=="intvec" or ty=="intmat" )
857   {                                                 //case: id = intvec
858      if ( ty == "int" ) { intmat S =id; }
859      else { intmat S = intmat(id); }
860      intvec i = S[1..nrows(S),1..ncols(S)];
861      kill id;
862      intvec id = i;
863   }
864//------------------- consider intvec v and empty sum  -----------------------
865   if( tt!=0 )
866   {
867      for (j=1; j<=size(v); j++)
868      {
869         if ( v[j] <= 0 or v[j] > size(id) )         //v outside range of id
870         {
871            return(0);                               //empty sum is 0
872         }
873      }
874      id = id[v];                                    //consider part of id
875   }                                                 //corresponding to v
876
878   n = size(id);
879   if (n>10)
880   { return( lsum(n div 2,list(id[1..(n div 2)])) + lsum(n-n div 2, list(id[(n div 2+1)..n])) );
881   }
882   else
883   { def Summe=id[1];
884     for (int lauf=2;lauf<=n;lauf++)
885     { Summe=Summe+id[lauf];
886     }
887     return(Summe);
888   }
889}
890example
891{  "EXAMPLE:"; echo = 2;
892   ring r1 = 0,(x,y,z),dp;
893   vector pv = [xy,xz,yz,x2,y2,z2];
894   sum(pv);
895   sum(pv,2..5);
896   matrix M[2][3] = 1,x,2,y,3,z;
897   intvec w=2,4,6;
898   sum(M,w);
899   intvec iv = 1,2,3,4,5,6,7,8,9;
900   sum(iv,2..4);
901   iv = intvec(1..100);
902   sum(iv);
903   ring r2 = 0,(x(1..10)),dp;
904   sum(x(3..7),intvec(1,3,5));
905}
906///////////////////////////////////////////////////////////////////////////////
907
908
909///////////////////////////////////////////////////////////////////////////////
910
911proc watchdog(int i, string cmd)
912"USAGE:   watchdog(i,cmd); i integer, cmd string
913RETURN:  Result of cmd, if the result can be computed in i seconds.
914         Otherwise the computation is interrupted after i seconds,
915         the string "Killed" is returned and the global variable
916         'watchdog_interrupt' is defined.
917NOTE:    * the MP package must be enabled
918         * the current basering should not be watchdog_rneu, since
919           watchdog_rneu will be killed
920         * if there are variable names of the structure x(i) all
921           polynomials have to be put into eval(...) in order to be
922           interpreted correctly
923         * a second Singular process is started by this procedure
924EXAMPLE: example watchdog; shows an example
925"
926{
927  string rname=nameof(basering);
928  def rsave=basering;
929  if (defined(watchdog_rneu))
930  {
931    kill watchdog_rneu;
932  }
933  if ( i > 0 )
934  {
935    int j=10;
936    int k=999999;
937// fork, get the pid of the child and send it the command
939    open(l_fork);
940    execute("write(l_fork,quote(" + cmd + "));");
941
942
943// sleep in small, but growing intervals for appr. 1 second
944    while(j < k)
945    {
947      j = j + j;
948    }
949
950// sleep in intervals of one second
951    j = 1;
953    {
954      while (j < i)
955      {
957        j = j + 1;
958      }
959    }
960// check, whether we have a result, and return it
962    {
964      if (nameof(basering)!=rname)
965      {
966        def watchdog_rneu=basering;
967        setring rsave;
968        if (!defined(result))
969        {
970          def result=fetch(watchdog_rneu,result);
971        }
972      }
973      if(defined(watchdog_interrupt))
974      {
975        kill watchdog_interrupt;
976      }
977      close(l_fork);
978    }
979    else
980    {
981      string result="Killed";
982      if(!defined(watchdog_interrupt))
983      {
984        int watchdog_interrupt=1;
985        export watchdog_interrupt;
986      }
987      close(l_fork);
988    }
989    return(result);
990  }
991  else
992  {
993    ERROR("First argument of watchdog has to be a positive integer.");
994  }
995}
996example
997{ "EXAMPLE:"; echo=2;
998  ring r=0,(x,y,z),dp;
999  poly f=x^30+y^30;
1000  watchdog(1,"factorize(eval("+string(f)+"))");
1001  watchdog(100,"factorize(eval("+string(f)+"))");
1002}
1003///////////////////////////////////////////////////////////////////////////////
1004
1005proc deleteSublist(intvec v,list l)
1006"USAGE:   deleteSublist(v,l); intvec v; list l
1007         where the entries of the integer vector v correspond to the
1008         positions of the elements to be deleted
1009RETURN:  list without the deleted elements
1010EXAMPLE: example deleteSublist; shows an example"
1011{
1012  list k;
1013  int i,j,skip;
1014  j=1;
1015  skip=0;
1016  intvec vs=sort(v)[1];
1017  for ( i=1 ; i <=size(vs) ; i++)
1018  {
1019    while ((j+skip) < vs[i])
1020    {
1021      k[j] = l[j+skip];
1022      j++;
1023    }
1024    skip++;
1025  }
1026  if(vs[size(vs)]<size(l))
1027  {
1028    k=k+list(l[(vs[size(vs)]+1)..size(l)]);
1029  }
1030  return(k);
1031}
1032example
1033{ "EXAMPLE:"; echo=2;
1034   list l=1,2,3,4,5;
1035   intvec v=1,3,4;
1036   l=deleteSublist(v,l);
1037   l;
1038}
1039
1040///////////////////////////////////////////////////////////////////////////////
1041proc primecoeffs(J, list #)
1042"USAGE:   primecoeffs(J[,p]); J any type which can be converted to a matrix
1043         e.g. ideal, matrix, vector, module, int, intvec
1044         p = integer
1045COMPUTE: primefactors <= p of coeffs of J (default p = 32003)
1046RETURN:  a list, say l, of two intvectors:@*
1047         l[1] : the different primefactors of all coefficients of J@*
1048         l[2] : the different remaining factors
1049EXAMPLE: example primecoeffs; shows an example
1050"
1051{
1052   int q,ii,n,mark;;
1053   if (size(#) == 0)
1054   {
1055      q=32003;
1056   }
1057   else
1058   {
1059     if( typeof(#[1]) != "int")
1060     {
1061       ERROR("2nd parameter must be of type int"+newline);
1062     }
1063     q=#[1];
1064     if (q > 32003) { q = 32003; }
1065   }
1066
1067   if (defined(basering) == 0)
1068   {
1069     mark=1;
1070     ring r = 0,x,dp;
1071   }
1072   def I = ideal(matrix(J));
1073   poly p = product(maxideal(1));
1074   matrix Coef=coef(I[1],p);
1075   ideal id, jd, rest;
1076   intvec v,re;
1077   list result,l;
1078   for(ii=2; ii<=ncols(I); ii++)
1079   {
1080     Coef=concat(Coef,coef(I[ii],p));
1081   }
1082   id = Coef[2,1..ncols(Coef)];
1083   id = simplify(id,6);
1084   for (ii=1; ii<=size(id); ii++)
1085   {
1086     l = primefactors(number(id[ii]),q);
1087     list jdl=l[1];
1088     jd = jd,jdl[1..size(jdl)];
1089     kill jdl;
1090     rest=rest, l[3];
1091   }
1092   jd = simplify(jd,6);
1093   for (ii=1; ii<=size(jd); ii++)
1094   {
1095     v[ii]=int(jd[ii]);
1096   }
1097   v = sort(v)[1];
1098   rest = simplify(rest,6);
1099   id = sort(id)[1];
1100   if (mark)
1101   {
1102     for (ii=1; ii<=size(rest); ii++)
1103     {
1104       re[ii] = int(rest[ii]);
1105     }
1106     result = v,re;
1107   }
1108   else
1109   {
1110      result = v,rest;
1111   }
1112   return(result);
1113}
1114example
1115{ "EXAMPLE:"; echo = 2;
1116   primecoeffs(intvec(7*8*121,7*8));"";
1117   ring r = 0,(b,c,t),dp;
1118   ideal I = -13b6c3t+4b5c4t,-10b4c2t-5b4ct2;
1119   primecoeffs(I);
1120}
1121///////////////////////////////////////////////////////////////////////////////
1122proc timeFactorize(poly i,list #)
1123"USAGE:  timeFactorize(p,d); poly p , integer d
1124RETURN:  factorize(p) if the factorization finished after d-1
1125         seconds otherwhise f is considered to be irreducible
1126EXAMPLE: example timeFactorize; shows an example
1127"
1128{
1129  def P=basering;
1130  if (size(#) > 0)
1131  {
1132    if ((typeof(#[1]) == "int")&&(#[1]))
1133    {
1134      int wait = #[1];
1135      int j = 10;
1136
1137      string bs = nameof(basering);
1139      open(l_fork);
1140      write(l_fork, quote(timeFactorize(eval(i))));
1141
1142      // sleep in small intervalls for appr. one second
1143      if (wait > 0)
1144      {
1145        while(j < 1000000)
1146        {
1148          j = j + j;
1149        }
1150      }
1151
1152      // sleep in intervalls of one second from now on
1153      j = 1;
1154      while (j < wait)
1155      {
1157        j = j + 1;
1158      }
1159
1161      {
1163        if (bs != nameof(basering))
1164        {
1165          def PP = basering;
1166          setring P;
1167          def result = imap(PP, result);
1168          kill PP;
1169        }
1170        close(l_fork);
1171      }
1172      else
1173      {
1174        list result;
1175        intvec v=1,1;
1176        result[1]=list(1,i);
1177        result[2]=v;
1178        close(l_fork);
1179      }
1180      return (result);
1181    }
1182  }
1183  return(factorH(i));
1184}
1185example
1186{ "EXAMPLE:"; echo = 2;
1187   ring r=0,(x,y),dp;
1188   poly p=((x2+y3)^2+xy6)*((x3+y2)^2+x10y);
1189   p=p^2;
1190   //timeFactorize(p,2);
1191   //timeFactorize(p,20);
1192}
1193
1194proc timeStd(ideal i,list #)
1195"USAGE:  timeStd(i,d), i ideal, d integer
1196RETURN:  std(i) if the standard basis computation finished after
1197         d-1 seconds and i otherwhise
1198EXAMPLE: example timeStd; shows an example
1199"
1200{
1201  def P=basering;
1202  if (size(#) > 0)
1203  {
1204    if ((typeof(#[1]) == "int")&&(#[1]))
1205    {
1206      int wait = #[1];
1207      int j = 10;
1208
1209      string bs = nameof(basering);
1211      open(l_fork);
1212      write(l_fork, quote(system("pid")));
1214      write(l_fork, quote(timeStd(eval(i))));
1215
1216      // sleep in small intervalls for appr. one second
1217      if (wait > 0)
1218      {
1219        while(j < 1000000)
1220        {
1222          j = j + j;
1223        }
1224      }
1225      j = 1;
1226      while (j < wait)
1227      {
1229        j = j + 1;
1230      }
1232      {
1234        if (bs != nameof(basering))
1235        {
1236          def PP = basering;
1237          setring P;
1238          def result = imap(PP, result);
1239          kill PP;
1240        }
1241        close(l_fork);
1242      }
1243      else
1244      {
1245        ideal result=i;
1246        close(l_fork);
1247      }
1248      return (result);
1249    }
1250  }
1251  return(std(i));
1252}
1253example
1254{ "EXAMPLE:"; echo = 2;
1255   ring r=32003,(a,b,c,d,e),dp;
1256   int n=7;
1257   ideal i=
1258   a^n-b^n,
1259   b^n-c^n,
1260   c^n-d^n,
1261   d^n-e^n,
1262   a^(n-1)*b+b^(n-1)*c+c^(n-1)*d+d^(n-1)*e+e^(n-1)*a;
1263   def i1=timeStd(i,1);
1264   def i2=timeStd(i,20);
1265   listvar();
1266}
1267
1268proc factorH(poly p)
1269"USAGE:  factorH(p)  p poly
1270RETURN:  factorize(p)
1271NOTE:    changes variables to make the last variable the principal
1272         one in the multivariate factorization and factorizes then
1273         the polynomial
1274EXAMPLE: example factorH; shows an example
1275"
1276{
1277   def R=basering;
1278   int i,j;
1279   int n=1;
1280   int d=nrows(coeffs(p,var(1)));
1281   for(i=1;i<=nvars(R);i++)
1282   {
1283      j=nrows(coeffs(p,var(i)));
1284      if(d>j)
1285      {
1286         n=i;
1287         d=j;
1288      }
1289   }
1290   ideal ma=maxideal(1);    //die letzte Variable ist die Hauptvariable
1291   ma[nvars(R)]=var(n);
1292   ma[n]=var(nvars(R));
1293   map phi=R,ma;
1294   list fac=factorize(phi(p));
1295   list re=phi(fac);
1296   return(re);
1297}
1298example
1299{ "EXAMPLE:"; echo = 2;
1300  system("random",992851144);
1301  ring r=32003,(x,y,z,w,t),lp;
1302  poly p=y2w9+yz7t-yz5w4-z2w4t4-w8t3;
1303  factorize(p);  //fast
1304  system("random",992851262);
1305  //factorize(p);  //slow
1306  system("random",992851262);
1307  factorH(p);
1308}
Note: See TracBrowser for help on using the repository browser.