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

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