source: git/Singular/LIB/general.lib @ 2ab830

spielwiese
Last change on this file since 2ab830 was 29c4ee, checked in by Hans Schoenemann <hannes@…>, 13 years ago
fix tr. 309 git-svn-id: file:///usr/local/Singular/svn/trunk@13814 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 37.4 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.
203SEE ALSO: prime
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
225SEE ALSO: prime
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\"
293           and except loaded procedures
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
768   // Note that in general: lead(sort(M)) != sort(lead(M)), e.g:
769   module M = [0, 1, 1, 0], [1, 0, 0, 1]; M;
770   sort(lead(M), "c, dp")[1];
771   lead(sort(M, "c, dp")[1]);
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/2,list(l[1..(n/2)])) + lsum(n-n/2, list(l[(n/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
877//-------------------------- finally, add objects  ---------------------------
878   n = size(id);
879   if (n>10)
880   { return( lsum(n/2,list(id[1..(n/2)])) + lsum(n-n/2, list(id[(n/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 we do not have MP-links, watchdog cannot be used
934  if (system("with","MP"))
935  {
936    if ( i > 0 )
937    {
938      int j=10;
939      int k=999999;
940// fork, get the pid of the child and send it the command
941      link l_fork="MPtcp:fork";
942      open(l_fork);
943      write(l_fork,quote(system("pid")));
944      int pid=read(l_fork);
945      execute("write(l_fork,quote(" + cmd + "));");
946
947
948// sleep in small, but growing intervals for appr. 1 second
949      while(j < k)
950      {
951        if (status(l_fork, "read", "ready", j)) {break;}
952        j = j + j;
953      }
954
955// sleep in intervals of one second
956      j = 1;
957      if (!status(l_fork,"read","ready"))
958      {
959        while (j < i)
960        {
961          if (status(l_fork, "read", "ready", k)) {break;}
962          j = j + 1;
963        }
964      }
965// check, whether we have a result, and return it
966      if (status(l_fork, "read", "ready"))
967      {
968        def result = read(l_fork);
969        if (nameof(basering)!=rname)
970        {
971          def watchdog_rneu=basering;
972          setring rsave;
973          if (!defined(result))
974          {
975            def result=fetch(watchdog_rneu,result);
976          }
977        }
978        if(defined(watchdog_interrupt))
979        {
980          kill watchdog_interrupt;
981        }
982        close(l_fork);
983      }
984      else
985      {
986        string result="Killed";
987        if(!defined(watchdog_interrupt))
988        {
989          int watchdog_interrupt=1;
990          export watchdog_interrupt;
991        }
992        close(l_fork);
993        j = system("sh","kill " + string(pid));
994      }
995      return(result);
996    }
997    else
998    {
999      ERROR("First argument of watchdog has to be a positive integer.");
1000    }
1001  }
1002  else
1003  {
1004    ERROR("MP-support is not enabled in this version of Singular.");
1005  }
1006}
1007example
1008{ "EXAMPLE:"; echo=2;
1009  ring r=0,(x,y,z),dp;
1010  poly f=x^30+y^30;
1011  watchdog(1,"factorize(eval("+string(f)+"))");
1012  watchdog(100,"factorize(eval("+string(f)+"))");
1013}
1014///////////////////////////////////////////////////////////////////////////////
1015
1016proc deleteSublist(intvec v,list l)
1017"USAGE:   deleteSublist(v,l); intvec v; list l
1018         where the entries of the integer vector v correspond to the
1019         positions of the elements to be deleted
1020RETURN:  list without the deleted elements
1021EXAMPLE: example deleteSublist; shows an example"
1022{
1023  list k;
1024  int i,j,skip;
1025  j=1;
1026  skip=0;
1027  intvec vs=sort(v)[1];
1028  for ( i=1 ; i <=size(vs) ; i++)
1029  {
1030    while ((j+skip) < vs[i])
1031    {
1032      k[j] = l[j+skip];
1033      j++;
1034    }
1035    skip++;
1036  }
1037  if(vs[size(vs)]<size(l))
1038  {
1039    k=k+list(l[(vs[size(vs)]+1)..size(l)]);
1040  }
1041  return(k);
1042}
1043example
1044{ "EXAMPLE:"; echo=2;
1045   list l=1,2,3,4,5;
1046   intvec v=1,3,4;
1047   l=deleteSublist(v,l);
1048   l;
1049}
1050
1051///////////////////////////////////////////////////////////////////////////////
1052proc primecoeffs(J, list #)
1053"USAGE:   primecoeffs(J[,p]); J any type which can be converted to a matrix
1054         e.g. ideal, matrix, vector, module, int, intvec
1055         p = integer
1056COMPUTE: primefactors <= p of coeffs of J (default p = 32003)
1057RETURN:  a list, say l, of two intvectors:@*
1058         l[1] : the different primefactors of all coefficients of J@*
1059         l[2] : the different remaining factors
1060EXAMPLE: example primecoeffs; shows an example
1061"
1062{
1063   int q,ii,n,mark;;
1064   if (size(#) == 0)
1065   {
1066      q=32003;
1067   }
1068   else
1069   {
1070     if( typeof(#[1]) != "int")
1071     {
1072       ERROR("2nd parameter must be of type int"+newline);
1073     }
1074     q=#[1];
1075     if (q > 32003) { q = 32003; }
1076   }
1077
1078   if (defined(basering) == 0)
1079   {
1080     mark=1;
1081     ring r = 0,x,dp;
1082   }
1083   def I = ideal(matrix(J));
1084   poly p = product(maxideal(1));
1085   matrix Coef=coef(I[1],p);
1086   ideal id, jd, rest;
1087   intvec v,re;
1088   list result,l;
1089   for(ii=2; ii<=ncols(I); ii++)
1090   {
1091     Coef=concat(Coef,coef(I[ii],p));
1092   }
1093   id = Coef[2,1..ncols(Coef)];
1094   id = simplify(id,6);
1095   for (ii=1; ii<=size(id); ii++)
1096   {
1097     l = primefactors(number(id[ii]),q);
1098     list jdl=l[1];
1099     jd = jd,jdl[1..size(jdl)];
1100     kill jdl;
1101     rest=rest, l[3];
1102   }
1103   jd = simplify(jd,6);
1104   for (ii=1; ii<=size(jd); ii++)
1105   {
1106     v[ii]=int(jd[ii]);
1107   }
1108   v = sort(v)[1];
1109   rest = simplify(rest,6);
1110   id = sort(id)[1];
1111   if (mark)
1112   {
1113     for (ii=1; ii<=size(rest); ii++)
1114     {
1115       re[ii] = int(rest[ii]);
1116     }
1117     result = v,re;
1118   }
1119   else
1120   {
1121      result = v,rest;
1122   }
1123   return(result);
1124}
1125example
1126{ "EXAMPLE:"; echo = 2;
1127   primecoeffs(intvec(7*8*121,7*8));"";
1128   ring r = 0,(b,c,t),dp;
1129   ideal I = -13b6c3t+4b5c4t,-10b4c2t-5b4ct2;
1130   primecoeffs(I);
1131}
1132///////////////////////////////////////////////////////////////////////////////
1133proc timeFactorize(poly i,list #)
1134"USAGE:  timeFactorize(p,d); poly p , integer d
1135RETURN:  factorize(p) if the factorization finished after d-1
1136         seconds otherwhise f is considered to be irreducible
1137EXAMPLE: example timeFactorize; shows an example
1138"
1139{
1140  def P=basering;
1141  if (size(#) > 0)
1142  {
1143    if (system("with", "MP"))
1144    {
1145      if ((typeof(#[1]) == "int")&&(#[1]))
1146      {
1147        int wait = #[1];
1148        int j = 10;
1149
1150        string bs = nameof(basering);
1151        link l_fork = "MPtcp:fork";
1152        open(l_fork);
1153        write(l_fork, quote(system("pid")));
1154        int pid = read(l_fork);
1155        write(l_fork, quote(timeFactorize(eval(i))));
1156
1157        // sleep in small intervalls for appr. one second
1158        if (wait > 0)
1159        {
1160          while(j < 1000000)
1161          {
1162            if (status(l_fork, "read", "ready", j)) {break;}
1163            j = j + j;
1164          }
1165        }
1166
1167        // sleep in intervalls of one second from now on
1168        j = 1;
1169        while (j < wait)
1170        {
1171          if (status(l_fork, "read", "ready", 1000000)) {break;}
1172          j = j + 1;
1173        }
1174
1175        if (status(l_fork, "read", "ready"))
1176        {
1177          def result = read(l_fork);
1178          if (bs != nameof(basering))
1179          {
1180            def PP = basering;
1181            setring P;
1182            def result = imap(PP, result);
1183            kill PP;
1184          }
1185          kill l_fork;
1186        }
1187        else
1188        {
1189          list result;
1190          intvec v=1,1;
1191          result[1]=list(1,i);
1192          result[2]=v;
1193          j = system("sh", "kill " + string(pid));
1194        }
1195        return (result);
1196      }
1197    }
1198  }
1199  return(factorH(i));
1200}
1201example
1202{ "EXAMPLE:"; echo = 2;
1203   ring r=0,(x,y),dp;
1204   poly p=((x2+y3)^2+xy6)*((x3+y2)^2+x10y);
1205   p=p^2;
1206   //timeFactorize(p,2);
1207   //timeFactorize(p,20);
1208}
1209
1210proc timeStd(ideal i,list #)
1211"USAGE:  timeStd(i,d), i ideal, d integer
1212RETURN:  std(i) if the standard basis computation finished after
1213         d-1 seconds and i otherwhise
1214EXAMPLE: example timeStd; shows an example
1215"
1216{
1217  def P=basering;
1218  if (size(#) > 0)
1219  {
1220    if (system("with", "MP"))
1221    {
1222      if ((typeof(#[1]) == "int")&&(#[1]))
1223      {
1224        int wait = #[1];
1225        int j = 10;
1226
1227        string bs = nameof(basering);
1228        link l_fork = "MPtcp:fork";
1229        open(l_fork);
1230        write(l_fork, quote(system("pid")));
1231        int pid = read(l_fork);
1232        write(l_fork, quote(timeStd(eval(i))));
1233
1234        // sleep in small intervalls for appr. one second
1235        if (wait > 0)
1236        {
1237          while(j < 1000000)
1238          {
1239            if (status(l_fork, "read", "ready", j)) {break;}
1240            j = j + j;
1241          }
1242        }
1243        j = 1;
1244        while (j < wait)
1245        {
1246          if (status(l_fork, "read", "ready", 1000000)) {break;}
1247          j = j + 1;
1248        }
1249        if (status(l_fork, "read", "ready"))
1250        {
1251          def result = read(l_fork);
1252          if (bs != nameof(basering))
1253          {
1254            def PP = basering;
1255            setring P;
1256            def result = imap(PP, result);
1257            kill PP;
1258          }
1259          kill l_fork;
1260        }
1261        else
1262        {
1263          ideal result=i;
1264          j = system("sh", "kill " + string(pid));
1265        }
1266        return (result);
1267      }
1268    }
1269  }
1270  return(std(i));
1271}
1272example
1273{ "EXAMPLE:"; echo = 2;
1274   ring r=32003,(a,b,c,d,e),dp;
1275   int n=6;
1276   ideal i=
1277   a^n-b^n,
1278   b^n-c^n,
1279   c^n-d^n,
1280   d^n-e^n,
1281   a^(n-1)*b+b^(n-1)*c+c^(n-1)*d+d^(n-1)*e+e^(n-1)*a;
1282   timeStd(i,2);
1283   timeStd(i,20);
1284}
1285
1286proc factorH(poly p)
1287"USAGE:  factorH(p)  p poly
1288RETURN:  factorize(p)
1289NOTE:    changes variables to make the last variable the principal
1290         one in the multivariate factorization and factorizes then
1291         the polynomial
1292EXAMPLE: example factorH; shows an example
1293"
1294{
1295   def R=basering;
1296   int i,j;
1297   int n=1;
1298   int d=nrows(coeffs(p,var(1)));
1299   for(i=1;i<=nvars(R);i++)
1300   {
1301      j=nrows(coeffs(p,var(i)));
1302      if(d>j)
1303      {
1304         n=i;
1305         d=j;
1306      }
1307   }
1308   ideal ma=maxideal(1);    //die letzte Variable ist die Hauptvariable
1309   ma[nvars(R)]=var(n);
1310   ma[n]=var(nvars(R));
1311   map phi=R,ma;
1312   list fac=factorize(phi(p));
1313   list re=phi(fac);
1314   return(re);
1315}
1316example
1317{ "EXAMPLE:"; echo = 2;
1318  system("random",992851144);
1319  ring r=32003,(x,y,z,w,t),lp;
1320  poly p=y2w9+yz7t-yz5w4-z2w4t4-w8t3;
1321  factorize(p);  //fast
1322  system("random",992851262);
1323  //factorize(p);  //slow
1324  system("random",992851262);
1325  factorH(p);
1326}
Note: See TracBrowser for help on using the repository browser.