source:git/Singular/LIB/general.lib

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