//GMG, last modified 18.6.99 //anne, added deleteSublist and watchdog 12.12.2000 //eric, added absValue 11.04.2002 /////////////////////////////////////////////////////////////////////////////// version="$Id: general.lib,v 1.66 2009-07-28 14:17:00 Singular Exp $"; category="General purpose"; info=" LIBRARY: general.lib Elementary Computations of General Type PROCEDURES: A_Z(\"a\",n); string a,b,... of n comma separated letters ASCII([n,m]); string of printable ASCII characters (number n to m) absValue(c); absolute value of c binomial(n,m[,../..]); n choose m (type int), [type bigint] deleteSublist(iv,l); delete entries given by iv from list l factorial(n[,../..]); n factorial (=n!) (type int), [type bigint] fibonacci(n[,p]); nth Fibonacci number [char p] kmemory([n[,v]]); active [allocated] memory in kilobyte killall(); kill all user-defined variables number_e(n); compute exp(1) up to n decimal digits number_pi(n); compute pi (area of unit circle) up to n digits primes(n,m); intvec of primes p, n<=p<=m product(../..[,v]); multiply components of vector/ideal/...[indices v] sort(ideal/module); sort generators according to monomial ordering sum(vector/id/..[,v]); add components of vector/ideal/...[with indices v] watchdog(i,cmd); only wait for result of command cmd for i seconds which(command); search for command and return absolute path, if found primecoeffs(J[,q]); primefactors <= min(p,32003) of coeffs of J primefactors(n[,p]); primefactors <= min(p,32003) of n timeStd(i,d) std(i) if the standard basis computation finished after d-1 seconds and i otherwhise timeFactorize(p,d) factorize(p) if the factorization finished after d-1 seconds otherwhise f is considered to be irreducible factorH(p) changes variables to become the last variable the principal one in the multivariate factorization and factorizes then the polynomial (parameters in square brackets [] are optional) "; LIB "inout.lib"; LIB "poly.lib"; LIB "matrix.lib"; /////////////////////////////////////////////////////////////////////////////// proc A_Z (string s,int n) "USAGE: A_Z(\"a\",n); a any letter, n integer (-26<= n <=26, !=0) RETURN: string of n small (if a is small) or capital (if a is capital) letters, comma separated, beginning with a, in alphabetical order (or reverse alphabetical order if n<0) EXAMPLE: example A_Z; shows an example " { if ( n>=-26 and n<=26 and n!=0 ) { string alpha = "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,"+ "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,"+ "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,"+ "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"; int ii; int aa; for(ii=1; ii<=51; ii=ii+2) { if( alpha[ii]==s ) { aa=ii; } } if ( aa==0) { for(ii=105; ii<=155; ii=ii+2) { if( alpha[ii]==s ) { aa=ii; } } } if( aa!=0 ) { string out; if (n > 0) { out = alpha[aa,2*(n)-1]; return (out); } if (n < 0) { string beta = "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,"+ "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,"+ "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,"+ "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"; if ( aa < 52 ) { aa=52-aa; } if ( aa > 104 ) { aa=260-aa; } out = beta[aa,2*(-n)-1]; return (out); } } } } example { "EXAMPLE:"; echo = 2; A_Z("c",5); A_Z("Z",-5); string sR = "ring R = (0,"+A_Z("A",6)+"),("+A_Z("a",10)+"),dp;"; sR; execute(sR); R; } /////////////////////////////////////////////////////////////////////////////// proc ASCII (list #) "USAGE: ASCII([n,m]); n,m= integers (32 <= n <= m <= 126) RETURN: string of printable ASCII characters (no native language support)@* ASCII(): string of all ASCII characters with its numbers,@* ASCII(n): n-th ASCII character@* ASCII(n,m): n-th up to m-th ASCII character (inclusive) EXAMPLE: example ASCII; shows an example " { string s1 = " ! \" # $ % & ' ( ) * + , - . 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 / 0 1 2 3 4 5 6 7 8 9 : ; < = 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 > ? @ A B C D E F G H I J K L 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 M N O P Q R S T U V W X Y Z [ 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 \\ ] ^ _ ` a b c d e f g h i j 92 93 94 95 96 97 98 99 100 101 102 103 104 105 10 k l m n o p q r s t u v w x y 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 z { | } ~ 122 123 124 125 126 "; string s2 = " !\"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\\]^_`abcdefghijklmnopqrstuvwxyz{|}~"; if ( size(#) == 0 ) { return(s1); } if ( size(#) == 1 ) { return( s2[#[1]-31] ); } if ( size(#) == 2 ) { return( s2[#[1]-31,#[2]-#[1]+1] ); } } example { "EXAMPLE:"; echo = 2; ASCII();""; ASCII(42); ASCII(32,126); } /////////////////////////////////////////////////////////////////////////////// proc absValue(def c) "USAGE: absValue(c); c int, number or poly RETURN: absValue(c); the absolute value of c NOTE: absValue(c)=c if c>=0; absValue=-c if c<=0. @* So the function can be applied to any type, for which comparison @* operators are defined. EXAMPLE: example absValue; shows an example " { if (c>=0) { return(c); } else { return(-c); } } example { "EXAMPLE:"; echo = 2; ring r1 = 0,x,dp; absValue(-2002); poly f=-4; absValue(f); } /////////////////////////////////////////////////////////////////////////////// proc binomial (int n, int k) "USAGE: binomial(n,k); n,k integers RETURN: binomial(n,k); binomial coefficient n choose k @* - of type bigint (computed in characteristic 0) NOTE: In any characteristic, binomial(n,k) = coefficient of x^k in (1+x)^n SEE ALSO: prime EXAMPLE: example binomial; shows an example " { bigint l; bigint r=1; bigint kk=k; bigint nn=n; if ( k > n-k ) { k = n-k; } if ( k<=0 or k>n ) //trivial cases { r = (k==0)*r; } for (l=1; l<=kk; l=l+1 ) { r=r*(nn+1-l)/l; } return(r); } example { "EXAMPLE:"; echo = 2; binomial(200,100);""; //type bigint int n,k = 200,100; bigint b1 = binomial(n,k); ring r = 0,x,dp; poly b2 = coeffs((x+1)^n,x)[k+1,1]; //coefficient of x^k in (x+1)^n b1-b2; //b1 and b2 should coincide } /////////////////////////////////////////////////////////////////////////////// static proc binomp (int n, int k, int p) //computes binomial coefficient n choose k in basering of char p > 0 //binomial(n,k) = coefficient of x^k in (1+x)^n. //Let n=q*p^j, gcd(q,p)=1, then (1+x)^n = (1 + x^(p^j))^q. We have //binomial(n,k)=0 if k!=l*p^j and binomial(n,l*p^j) = binomial(q,l). //Do this reduction first. Then, in denominator and numerator //of defining formula for binomial coefficient, reduce those factors //mod p which are not divisible by p and cancel common factors p. Hence, //if n = h*p+r, k=l*p+s, r,s n-k ) { k = n-k; } if ( k<=0 or k>n) //trivial cases { r = (k==0)*r; } else { while ( q mod p == 0 ) { l = l*p; q = q div p; } //we have now n=q*l, l=p^j, gcd(q,p)=1; if (k mod l != 0 ) { r = 0; } else { l = k div l; n = q mod p; k = l mod p; //now 0<= k,n 1) { verb=#[2]; } } if ( verb != 0) { if ( n==0) { dbprint(printlevel-voice+3, "// memory used, at the moment, by active variables (kilobyte):"); } if ( n==1 ) { dbprint(printlevel-voice+3, "// total memory allocated, at the moment, by SINGULAR (kilobyte):"); } } return ((memory(n)+1023) div 1024); } example { "EXAMPLE:"; echo = 2; kmemory(); kmemory(1,1); } /////////////////////////////////////////////////////////////////////////////// proc killall "USAGE: killall(); (no parameter) killall(\"type_name\"); killall(\"not\", \"type_name\"); RETURN: killall(); kills all user-defined variables except loaded procedures, no return value. @* - killall(\"type_name\"); kills all user-defined variables, of type \"type_name\" @* - killall(\"not\", \"type_name\"); kills all user-defined variables, except those of type \"type_name\" and except loaded procedures @* - killall(\"not\", \"name_1\", \"name_2\", ...); kills all user-defined variables, except those of name \"name_i\" and except loaded procedures NOTE: killall should never be used inside a procedure EXAMPLE: example killall; shows an example AND KILLS ALL YOUR VARIABLES " { list @marie=names(Top); int j, no_kill, @joni; for ( @joni=1; @joni<=size(#); @joni++) { if (typeof(#[@joni]) != "string") { ERROR("Need string as " + string(@joni) + "th argument"); } } // kills all user-defined variables but not loaded procedures if( size(#)==0 ) { for ( @joni=size(@marie); @joni>0; @joni-- ) { if( @marie[@joni]!="LIB" and typeof(`@marie[@joni]`)!="proc" and typeof(`@marie[@joni]`)!="package") { kill `@marie[@joni]`; } } } else { // kills all user-defined variables if( size(#)==1 ) { // of type proc if( #[1] == "proc" ) { for ( @joni=size(@marie); @joni>0; @joni-- ) { if( (@marie[@joni]!="General") and (@marie[@joni]!="Top") and (@marie[@joni]!="killall") and (@marie[@joni]!="LIB") and ((typeof(`@marie[@joni]`)=="package") or (typeof(`@marie[@joni]`)=="proc"))) { if (defined(`@marie[@joni]`)) {kill `@marie[@joni]`;} } if (!defined(@joni)) break; } if (defined(General)) { @marie=names(General); for ( @joni=size(@marie); @joni>0; @joni-- ) { if(@marie[@joni]!="killall" and typeof(`@marie[@joni]`)=="proc") { kill General::`@marie[@joni]`; } } kill General::killall; } } else { // other types for ( @joni=size(@marie); @joni>2; @joni-- ) { if(typeof(`@marie[@joni]`)==#[1] and @marie[@joni]!="LIB" and typeof(`@marie[@joni]`)!="proc") { kill `@marie[@joni]`; } } } } else { // kills all user-defined variables whose name or type is not #i for ( @joni=size(@marie); @joni>2; @joni-- ) { if ( @marie[@joni] != "LIB" && @marie[@joni] != "Top" && typeof(`@marie[@joni]`) != "proc") { no_kill = 0; for (j=2; j<= size(#); j++) { if (typeof(`@marie[@joni]`)==#[j] or @marie[@joni] == #[j]) { no_kill = 1; break; } } if (! no_kill) { kill `@marie[@joni]`; } } if (!defined(@joni)) break; } } } } example { "EXAMPLE:"; echo = 2; ring rtest; ideal i=x,y,z; string str="hi"; int j = 3; export rtest,i,str,j; //this makes the local variables global listvar(); killall("ring"); // kills all rings listvar(); killall("not", "int"); // kills all variables except int's (and procs) listvar(); killall(); // kills all vars except loaded procs listvar(); } /////////////////////////////////////////////////////////////////////////////// proc number_e (int n) "USAGE: number_e(n); n integer RETURN: Euler number e=exp(1) up to n decimal digits (no rounding) @* - of type string if no basering of char 0 is defined @* - of type number if a basering of char 0 is defined DISPLAY: decimal format of e if printlevel > 0 (default:printlevel=0 ) NOTE: procedure uses algorithm of A.H.J. Sale EXAMPLE: example number_e; shows an example " { int i,m,s,t; intvec u,e; u[n+2]=0; e[n+1]=0; e=e+1; if( defined(basering) ) { if( char(basering)==0 ) { number r=2; t=1; } } string result = "2."; for( i=1; i<=n+1; i=i+1 ) { e = e*10; for( m=n+1; m>=1; m=m-1 ) { s = e[m]+u[m+1]; u[m] = s div (m+1); e[m] = s%(m+1); } result = result+string(u[1]); if( t==1 ) { r = r+number(u[1])/number(10)^i; } } if( t==1 ) { dbprint(printlevel-voice+2,"// "+result[1,n+1]); return(r); } return(result[1,n+1]); } example { "EXAMPLE:"; echo = 2; number_e(30);""; ring R = 0,t,lp; number e = number_e(30); e; } /////////////////////////////////////////////////////////////////////////////// proc number_pi (int n) "USAGE: number_pi(n); n positive integer RETURN: pi (area of unit circle) up to n decimal digits (no rounding) @* - of type string if no basering of char 0 is defined, @* - of type number, if a basering of char 0 is defined DISPLAY: decimal format of pi if printlevel > 0 (default:printlevel=0 ) NOTE: procedure uses algorithm of S. Rabinowitz EXAMPLE: example number_pi; shows an example " { int i,m,t,e,q,N; intvec r,p,B,Prelim; string result,prelim; N = (10*n) div 3 + 2; p[N+1]=0; p=p+2; r=p; for( i=1; i<=N+1; i=i+1 ) { B[i]=2*i-1; } if( defined(basering) ) { if( char(basering)==0 ) { number pi; number pri; t=1; } } for( i=0; i<=n; i=i+1 ) { p = r*10; e = p[N+1]; for( m=N+1; m>=2; m=m-1 ) { r[m] = e%B[m]; q = e div B[m]; e = q*(m-1)+p[m-1]; } r[1] = e%10; q = e div 10; if( q!=10 and q!=9 ) { result = result+prelim; Prelim = q; prelim = string(q); } if( q==9 ) { Prelim = Prelim,9; prelim = prelim+"9"; } if( q==10 ) { Prelim = (Prelim+1)-((Prelim+1) div 10)*10; for( m=size(Prelim); m>0; m=m-1) { prelim[m] = string(Prelim[m]); } result = result+prelim; if( t==1 ) { pi=pi+pri; } Prelim = 0; prelim = "0"; } if( t==1 ) { pi=pi+number(q)/number(10)^i; } } result = result,prelim[1]; result = "3."+result[2,n-1]; if( t==1 ) { dbprint(printlevel-voice+2,"// "+result); return(pi); } return(result); } example { "EXAMPLE:"; echo = 2; number_pi(11);""; ring r = (real,10),t,dp; number pi = number_pi(11); pi; } /////////////////////////////////////////////////////////////////////////////// proc primes (int n, int m) "USAGE: primes(n,m); n,m integers RETURN: intvec, consisting of all primes p, prime(n)<=p<=m, in increasing order if n<=m, resp. prime(m)<=p<=n, in decreasing order if m=2, else 2 EXAMPLE: example primes; shows an example " { int change; if ( n>m ) { change=n; n=m ; m=change; change=1; } int q,p = prime(m),prime(n); intvec v = q; q = q-1; while ( q>=p ) { q = prime(q); v = q,v; q = q-1; } if ( change==1 ) { v = v[size(v)..1]; } return(v); } example { "EXAMPLE:"; echo = 2; primes(50,100);""; intvec v = primes(37,1); v; } /////////////////////////////////////////////////////////////////////////////// proc product (id, list #) "USAGE: product(id[,v]); id ideal/vector/module/matrix/intvec/intmat/list, v intvec (default: v=1..number of entries of id) ASSUME: list members can be multiplied. RETURN: The product of all entries of id [with index given by v] of type depending on the entries of id. NOTE: If id is not a list, id is treated as a list of polys resp. integers. A module m is identified with the corresponding matrix M (columns of M generate m). @* If v is outside the range of id, we have the empty product and the result will be 1 (of type int). EXAMPLE: example product; shows an example " { //-------------------- initialization and special feature --------------------- int n,j,tt; string ty; //will become type of id list l; // We wish to allow something like product(x(1..10)) if x(1),...,x(10) are // variables. x(1..10) is a list of polys and enters the procedure with // id=x(1) and # a list with 9 polys, #[1]= x(2),...,#[9]= x(10). Hence, in // this case # is never empty. If an additional intvec v is given, // it is added to #, so we have to separate it first and make // the rest a list which has to be multiplied. int s = size(#); if( s!=0 ) { if ( typeof(#[s])=="intvec" or typeof(#[s])=="int") { intvec v = #[s]; tt=1; s=s-1; if ( s>0 ) { # = #[1..s]; } } } if ( s>0 ) { l = list(id)+#; kill id; list id = l; //case: id = list ty = "list"; n = size(id); } else { ty = typeof(id); if( ty == "list" ) { n = size(id); } } //------------------------------ reduce to 3 cases --------------------------- if( ty=="poly" or ty=="ideal" or ty=="vector" or ty=="module" or ty=="matrix" ) { ideal i = ideal(matrix(id)); kill id; ideal id = i; //case: id = ideal n = ncols(id); } if( ty=="int" or ty=="intvec" or ty=="intmat" ) { if ( ty == "int" ) { intmat S =id; } else { intmat S = intmat(id); } intvec i = S[1..nrows(S),1..ncols(S)]; kill id; intvec id = i; //case: id = intvec n = size(id); } //--------------- consider intvec v and empty product ----------------------- if( tt!=0 ) { for (j=1; j<=size(v); j++) { if ( v[j] <= 0 or v[j] > n ) //v outside range of id { return(1); //empty product is 1 } } id = id[v]; //consider part of id } //corresponding to v //--------------------- special case: one factor is zero --------------------- if ( typeof(id) == "ideal") { if( size(id) < ncols(id) ) { poly f; return(f); } } //-------------------------- finally, multiply objects ----------------------- n = size(id); def f(1) = id[1]; for( j=2; j<=n; j=j+1 ) { def f(j)=f(j-1)*id[j]; } return(f(n)); } example { "EXAMPLE:"; echo = 2; ring r= 0,(x,y,z),dp; ideal m = maxideal(1); product(m); product(m[2..3]); matrix M[2][3] = 1,x,2,y,3,z; product(M); intvec v=2,4,6; product(M,v); intvec iv = 1,2,3,4,5,6,7,8,9; v=1..5,7,9; product(iv,v); intmat A[2][3] = 1,1,1,2,2,2; product(A,3..5); } /////////////////////////////////////////////////////////////////////////////// proc sort (id, list #) "USAGE: sort(id[,v,o,n]); id = ideal/module/intvec/list(of intvec's or int's) @* sort may be called with 1, 2 or 3 arguments in the following way: @* sort(id[,v,n]); v=intvec of positive integers, n=integer, @* sort(id[,o,n]); o=string (any allowed ordstr of a ring), n=integer RETURN: a list l of two elements: @format l[1]: object of same type as input but sorted in the following way: - if id=ideal/module: generators of id are sorted w.r.t. intvec v (id[v[1]] becomes 1-st, id[v[2]] 2-nd element, etc.). If no v is present, id is sorted w.r.t. ordering o (if o is given) or w.r.t. actual monomial ordering (if no o is given): NOTE: generators with SMALLER(!) leading term come FIRST (e.g. sort(id); sorts backwards to actual monomial ordering) - if id=list of intvec's or int's: consider a list element, say id[1]=3,2,5, as exponent vector of the monomial x^3*y^2*z^5; the corresponding monomials are ordered w.r.t. intvec v (s.a.). If no v is present, the monomials are sorted w.r.t. ordering o (if o is given) or w.r.t. lexicographical ordering (if no o is given). The corresponding ordered list of exponent vectors is returned. (e.g. sort(id); sorts lexicographically, smaller int's come first) WARNING: Since negative exponents create the 0 polynomial in Singular, id should not contain negative integers: the result might not be as expected - if id=intvec: id is treated as list of integers - if n!=0 the ordering is inverse, i.e. w.r.t. v(size(v)..1) default: n=0 l[2]: intvec, describing the permutation of the input (hence l[2]=v if v is given (with positive integers)) @end format NOTE: If v is given id may be any simply indexed object (e.g. any list or string); if v[i]<0 and i<=size(id) v[i] is set internally to i; entries of v must be pairwise distinct to get a permutation id. Zero generators of ideal/module are deleted If 'o' is given, the input is sorted by considering leading terms w.r.t. the new ring ordering given by 'o' EXAMPLE: example sort; shows an example " { int ii,jj,s,n = 0,0,1,0; intvec v; if ( defined(basering) ) { def P = basering; } if ( size(#)==0 and (typeof(id)=="ideal" or typeof(id)=="module" or typeof(id)=="matrix")) { id = simplify(id,2); for ( ii=1; ii=1 and (typeof(id)=="ideal" or typeof(id)=="module" or typeof(id)=="matrix") ) { if ( typeof(#[1])=="string" ) { execute("ring r1 =("+charstr(P)+"),("+varstr(P)+"),("+#[1]+");"); def i = imap(P,id); v = sortvec(i); setring P; n=2; } } if ( typeof(id)=="intvec" or typeof(id)=="list" and n==0 ) { string o; if ( size(#)==0 ) { o = "lp"; n=1; } if ( size(#)>=1 ) { if ( typeof(#[1])=="string" ) { o = #[1]; n=1; } } } if ( typeof(id)=="intvec" or typeof(id)=="list" and n==1 ) { if ( typeof(id)=="list" ) { for (ii=1; ii<=size(id); ii++) { if (typeof(id[ii]) != "intvec" and typeof(id[ii]) != "int") { ERROR("// list elements must be intvec/int"); } else { s=size(id[ii])*(s < size(id[ii])) + s*(s >= size(id[ii])); } } } execute("ring r=0,x(1..s),("+o+");"); ideal i; poly f; for (ii=1; ii<=size(id); ii++) { f=1; for (jj=1; jj<=size(id[ii]); jj++) { f=f*x(jj)^(id[ii])[jj]; } i[ii]=f; } v = sort(i)[2]; } if ( size(#)!=0 and n==0 ) { v = #[1]; } if( size(#)==2 ) { if ( #[2] != 0 ) { v = v[size(v)..1]; } } s = size(v); if( size(id) < s ) { s = size(id); } def m = id; if ( size(m) != 0 ) { for ( jj=1; jj<=s; jj++) { if ( v[jj]<=0 ) { v[jj]=jj; } m[jj] = id[v[jj]]; } } if ( v == 0 ) { v = 1; } list L=m,v; return(L); } example { "EXAMPLE:"; echo = 2; ring r0 = 0,(x,y,z,t),lp; ideal i = x3,z3,xyz; sort(i); //sorts using lex ordering, smaller polys come first sort(i,3..1); sort(i,"ls")[1]; //sort w.r.t. negative lex ordering intvec v =1,10..5,2..4;v; sort(v)[1]; // sort v lexicographically sort(v,"Dp",1)[1]; // sort v w.r.t (total sum, reverse lex) // Note that in general: lead(sort(M)) != sort(lead(M)), e.g: module M = [0, 1, 1, 0], [1, 0, 0, 1]; M; sort(lead(M), "c, dp")[1]; lead(sort(M, "c, dp")[1]); // In order to sort M wrt a NEW ordering by considering OLD leading // terms use one of the following equivalent commands: module( M[ sort(lead(M), "c,dp")[2] ] ); sort( M, sort(lead(M), "c,dp")[2] )[1]; } /////////////////////////////////////////////////////////////////////////////// static proc lsum (int n, list l) { if (n>10) { return( lsum(n/2,list(l[1..(n/2)])) + lsum(n-n/2, list(l[(n/2+1)..n])) ); } else { def Summe=l[1]; for (int i=2;i<=n;i++) { Summe=Summe+l[i]; } return(Summe); } } /////////////////////////////////////////////////////////////////////////////// proc sum (id, list #) "USAGE: sum(id[,v]); id ideal/vector/module/matrix/intvec/intmat/list, v intvec (default: v=1..number of entries of id) ASSUME: list members can be added. RETURN: The sum of all entries of id [with index given by v] of type depending on the entries of id. NOTE: If id is not a list, id is treated as a list of polys resp. integers. A module m is identified with the corresponding matrix M (columns of M generate m). @* If v is outside the range of id, we have the empty sum and the result will be 0 (of type int). EXAMPLE: example sum; shows an example " { //-------------------- initialization and special feature --------------------- int n,j,tt; string ty; // will become type of id list l; // We wish to allow something like sum(x(1..10)) if x(1),...,x(10) are // variables. x(1..10) is a list of polys and enters the procedure with // id=x(1) and # a list with 9 polys, #[1]= x(2),...,#[9]= x(10). Hence, in // this case # is never empty. If an additional intvec v is given, // it is added to #, so we have to separate it first and make // the rest a list which has to be added. int s = size(#); if( s!=0 ) { if ( typeof(#[s])=="intvec" or typeof(#[s])=="int") { intvec v = #[s]; tt=1; s=s-1; if ( s>0 ) { # = #[1..s]; } } } if ( s>0 ) { l = list(id)+#; kill id; list id = l; //case: id = list ty = "list"; } else { ty = typeof(id); } //------------------------------ reduce to 3 cases --------------------------- if( ty=="poly" or ty=="ideal" or ty=="vector" or ty=="module" or ty=="matrix" ) { //case: id = ideal ideal i = ideal(matrix(id)); kill id; ideal id = simplify(i,2); //delete 0 entries } if( ty=="int" or ty=="intvec" or ty=="intmat" ) { //case: id = intvec if ( ty == "int" ) { intmat S =id; } else { intmat S = intmat(id); } intvec i = S[1..nrows(S),1..ncols(S)]; kill id; intvec id = i; } //------------------- consider intvec v and empty sum ----------------------- if( tt!=0 ) { for (j=1; j<=size(v); j++) { if ( v[j] <= 0 or v[j] > size(id) ) //v outside range of id { return(0); //empty sum is 0 } } id = id[v]; //consider part of id } //corresponding to v //-------------------------- finally, add objects --------------------------- n = size(id); if (n>10) { return( lsum(n/2,list(id[1..(n/2)])) + lsum(n-n/2, list(id[(n/2+1)..n])) ); } else { def Summe=id[1]; for (int lauf=2;lauf<=n;lauf++) { Summe=Summe+id[lauf]; } return(Summe); } } example { "EXAMPLE:"; echo = 2; ring r1 = 0,(x,y,z),dp; vector pv = [xy,xz,yz,x2,y2,z2]; sum(pv); sum(pv,2..5); matrix M[2][3] = 1,x,2,y,3,z; intvec w=2,4,6; sum(M,w); intvec iv = 1,2,3,4,5,6,7,8,9; sum(iv,2..4); iv = intvec(1..100); sum(iv); ring r2 = 0,(x(1..10)),dp; sum(x(3..7),intvec(1,3,5)); } /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// proc which (command) "USAGE: which(command); command = string expression RETURN: Absolute pathname of command, if found in search path. Empty string, otherwise. NOTE: Based on the Unix command 'which'. EXAMPLE: example which; shows an example " { int rs; int i; string fn = "which_" + string(system("pid")); string pn; string cmd; if( typeof(command) != "string") { return (pn); } if (system("uname") != "ix86-Win") { cmd = "which "; } else { // unfortunately, it does not take -path cmd = "type "; } i = system("sh", cmd + command + " > " + fn); pn = read(fn); if (system("uname") != "ix86-Win") { // TBC: Hmm... should parse output to get rid of 'command is ' pn[size(pn)] = ""; i = 1; while ((pn[i] != " ") and (pn[i] != "")) { i = i+1; } if (pn[i] == " ") {pn[i] = "";} rs = system("sh", "ls " + pn + " > " + fn + " 2>&1 "); } else { rs = 0; } i = system("sh", "rm " + fn); if (rs == 0) {return (pn);} else { print (command + " not found "); return (""); } } example { "EXAMPLE:"; echo = 2; which("sh"); } /////////////////////////////////////////////////////////////////////////////// proc watchdog(int i, string cmd) "USAGE: watchdog(i,cmd); i integer, cmd string RETURN: Result of cmd, if the result can be computed in i seconds. Otherwise the computation is interrupted after i seconds, the string "Killed" is returned and the global variable 'watchdog_interrupt' is defined. NOTE: * the MP package must be enabled * the current basering should not be watchdog_rneu, since watchdog_rneu will be killed * if there are variable names of the structure x(i) all polynomials have to be put into eval(...) in order to be interpreted correctly * a second Singular process is started by this procedure EXAMPLE: example watchdog; shows an example " { string rname=nameof(basering); def rsave=basering; if (defined(watchdog_rneu)) { kill watchdog_rneu; } // If we do not have MP-links, watchdog cannot be used if (system("with","MP")) { if ( i > 0 ) { int j=10; int k=999999; // fork, get the pid of the child and send it the command link l_fork="MPtcp:fork"; open(l_fork); write(l_fork,quote(system("pid"))); int pid=read(l_fork); execute("write(l_fork,quote(" + cmd + "));"); // sleep in small, but growing intervals for appr. 1 second while(j < k) { if (status(l_fork, "read", "ready", j)) {break;} j = j + j; } // sleep in intervals of one second j = 1; if (!status(l_fork,"read","ready")) { while (j < i) { if (status(l_fork, "read", "ready", k)) {break;} j = j + 1; } } // check, whether we have a result, and return it if (status(l_fork, "read", "ready")) { def result = read(l_fork); if (nameof(basering)!=rname) { def watchdog_rneu=basering; setring rsave; if (!defined(result)) { def result=fetch(watchdog_rneu,result); } } if(defined(watchdog_interrupt)) { kill watchdog_interrupt; } close(l_fork); } else { string result="Killed"; if(!defined(watchdog_interrupt)) { int watchdog_interrupt=1; export watchdog_interrupt; } close(l_fork); j = system("sh","kill " + string(pid)); } return(result); } else { ERROR("First argument of watchdog has to be a positive integer."); } } else { ERROR("MP-support is not enabled in this version of Singular."); } } example { "EXAMPLE:"; echo=2; ring r=0,(x,y,z),dp; poly f=x^30+y^30; watchdog(1,"factorize(eval("+string(f)+"))"); watchdog(100,"factorize(eval("+string(f)+"))"); } /////////////////////////////////////////////////////////////////////////////// proc deleteSublist(intvec v,list l) "USAGE: deleteSublist(v,l); intvec v; list l where the entries of the integer vector v correspond to the positions of the elements to be deleted RETURN: list without the deleted elements EXAMPLE: example deleteSublist; shows an example" { list k; int i,j,skip; j=1; skip=0; intvec vs=sort(v)[1]; for ( i=1 ; i <=size(vs) ; i++) { while ((j+skip) < vs[i]) { k[j] = l[j+skip]; j++; } skip++; } if(vs[size(vs)] 2147483647 ) //2147483647 max. integer representation { v = primes(2,p); number m; for( ii=1; ii<=size(v); ii++) { jj=0; while(1) { q = v[ii]; jj = jj+1; m = n/q; //divide n as often as possible if (denominator(m)!=1) { break; } n=m; } if( jj>1 ) { w1 = w1,v[ii]; //primes w2 = w2,jj-1; //powers } if( n <= 2147483647 ) { break; } } } if( n > 2147483647 ) //n is still too big { if( size(w1) >1 ) //at least 1 primefactor was found { w1 = w1[2..size(w1)]; w2 = w2[2..size(w2)]; } else //no primefactor was found { w1 = 1; w2 = 1; } l = w1,w2,n; return(l); } if( n <= 2147483647 ) //n is in inter range { num = int(n); kill n; int n = num; } } // --------------------------- trivial cases -------------------- if( n==0 ) { w1=1; w2=1; w3=0; l=w1,w2,w3; return(l); } if( n==1 ) { w3=1; if( size(w1) >1 ) //at least 1 primefactor was found { w1 = w1[2..size(w1)]; w2 = w2[2..size(w2)]; } else //no primefactor was found { w1 = 1; w2 = 1; } l=w1,w2,w3; return(l); } if ( prime(n)==n ) //note: prime(n) <= 32003 in Singular { //case n is a prime if (p > n) { w1=w1,n; w2=w2,1; w3=1; w1 = w1[2..size(w1)]; w2 = w2[2..size(w2)]; l=w1,w2,w3; return(l); } else { w3=n; if( size(w1) >1 ) //at least 1 primefactor was found { w1 = w1[2..size(w1)]; w2 = w2[2..size(w2)]; } else //no primefactor was found { w1 = 1; w2 = 1; } l=w1,w2,w3; return(l); } } else { if ( p >= n) { v = primes(q,n div 2 + 1); } else { v = primes(q,p); } //------------- search for primfactors <= last entry of v ------------ for(ii=1; ii<=size(v); ii++) { z=0; while( (n mod v[ii]) == 0 ) { z=z+1; n = n div v[ii]; } if (z!=0) { w1 = w1,v[ii]; //primes w2 = w2,z; //multiplicities } } } //--------------- case:at least 1 primefactor was found --------------- if( size(w1) >1 ) //at least 1 primefactor was found { w1 = w1[2..size(w1)]; w2 = w2[2..size(w2)]; } else //no primefactor was found { w1 = 1; w2 = 1; } w3 = n; l = w1,w2,w3; return(l); } example { "EXAMPLE:"; echo = 2; primefactors(7*8*121); ring r = 0,x,dp; primefactors(123456789100); } /////////////////////////////////////////////////////////////////////////////// proc primecoeffs(J, list #) "USAGE: primecoeffs(J[,p]); J any type which can be converted to a matrix e.g. ideal, matrix, vector, module, int, intvec p = integer COMPUTE: primefactors <= min(p,32003) of coeffs of J (default p = 32003) RETURN: a list, say l, of two intvectors:@* l[1] : the different primefactors of all coefficients of J@* l[2] : the different remaining factors NOTE: the procedure works for small integers only, just by testing all primes (not to be considered as serious prime factorization!) EXAMPLE: example primecoeffs; shows an example " { int q,ii,n,mark;; if (size(#) == 0) { q=32003; } else { if( typeof(#[1]) != "int") { ERROR("2nd parameter must be of type int"+newline); } q=#[1]; } if (defined(basering) == 0) { mark=1; ring r = 0,x,dp; } def I = ideal(matrix(J)); poly p = product(maxideal(1)); matrix Coef=coef(I[1],p); ideal id, jd, rest; intvec v,re; list result,l; for(ii=2; ii<=ncols(I); ii++) { Coef=concat(Coef,coef(I[ii],p)); } id = Coef[2,1..ncols(Coef)]; id = simplify(id,6); for (ii=1; ii<=size(id); ii++) { l = primefactors(number(id[ii]),q); jd = jd,l[1]; rest = rest,l[3]; } jd = simplify(jd,6); for (ii=1; ii<=size(jd); ii++) { v[ii]=int(jd[ii]); } v = sort(v)[1]; rest = simplify(rest,6); id = sort(id)[1]; if (mark) { for (ii=1; ii<=size(rest); ii++) { re[ii] = int(rest[ii]); } result = v,re; } else { result = v,rest; } return(result); } example { "EXAMPLE:"; echo = 2; primecoeffs(intvec(7*8*121,7*8));""; ring r = 0,(b,c,t),dp; ideal I = -13b6c3t+4b5c4t,-10b4c2t-5b4ct2; primecoeffs(I); } /////////////////////////////////////////////////////////////////////////////// proc timeFactorize(poly i,list #) "USAGE: timeFactorize(p,d); poly p , integer d RETURN: factorize(p) if the factorization finished after d-1 seconds otherwhise f is considered to be irreducible EXAMPLE: example timeFactorize; shows an example " { def P=basering; if (size(#) > 0) { if (system("with", "MP")) { if ((typeof(#[1]) == "int")&&(#[1])) { int wait = #[1]; int j = 10; string bs = nameof(basering); link l_fork = "MPtcp:fork"; open(l_fork); write(l_fork, quote(system("pid"))); int pid = read(l_fork); write(l_fork, quote(timeFactorize(eval(i)))); // sleep in small intervalls for appr. one second if (wait > 0) { while(j < 1000000) { if (status(l_fork, "read", "ready", j)) {break;} j = j + j; } } // sleep in intervalls of one second from now on j = 1; while (j < wait) { if (status(l_fork, "read", "ready", 1000000)) {break;} j = j + 1; } if (status(l_fork, "read", "ready")) { def result = read(l_fork); if (bs != nameof(basering)) { def PP = basering; setring P; def result = imap(PP, result); kill PP; } kill l_fork; } else { list result; intvec v=1,1; result[1]=list(1,i); result[2]=v; j = system("sh", "kill " + string(pid)); } return (result); } } } return(factorH(i)); } example { "EXAMPLE:"; echo = 2; ring r=0,(x,y),dp; poly p=((x2+y3)^2+xy6)*((x3+y2)^2+x10y); p=p^2; //timeFactorize(p,2); //timeFactorize(p,20); } proc timeStd(ideal i,list #) "USAGE: timeStd(i,d), i ideal, d integer RETURN: std(i) if the standard basis computation finished after d-1 seconds and i otherwhise EXAMPLE: example timeStd; shows an example " { def P=basering; if (size(#) > 0) { if (system("with", "MP")) { if ((typeof(#[1]) == "int")&&(#[1])) { int wait = #[1]; int j = 10; string bs = nameof(basering); link l_fork = "MPtcp:fork"; open(l_fork); write(l_fork, quote(system("pid"))); int pid = read(l_fork); write(l_fork, quote(timeStd(eval(i)))); // sleep in small intervalls for appr. one second if (wait > 0) { while(j < 1000000) { if (status(l_fork, "read", "ready", j)) {break;} j = j + j; } } j = 1; while (j < wait) { if (status(l_fork, "read", "ready", 1000000)) {break;} j = j + 1; } if (status(l_fork, "read", "ready")) { def result = read(l_fork); if (bs != nameof(basering)) { def PP = basering; setring P; def result = imap(PP, result); kill PP; } kill l_fork; } else { ideal result=i; j = system("sh", "kill " + string(pid)); } return (result); } } } return(std(i)); } example { "EXAMPLE:"; echo = 2; ring r=32003,(a,b,c,d,e),dp; int n=6; ideal i= a^n-b^n, b^n-c^n, c^n-d^n, d^n-e^n, a^(n-1)*b+b^(n-1)*c+c^(n-1)*d+d^(n-1)*e+e^(n-1)*a; timeStd(i,2); timeStd(i,20); } proc factorH(poly p) "USAGE: factorH(p) p poly RETURN: factorize(p) NOTE: changes variables to make the last variable the principal one in the multivariate factorization and factorizes then the polynomial EXAMPLE: example factorH; shows an example " { def R=basering; int i,j; int n=1; int d=nrows(coeffs(p,var(1))); for(i=1;i<=nvars(R);i++) { j=nrows(coeffs(p,var(i))); if(d>j) { n=i; d=j; } } ideal ma=maxideal(1); //die letzte Variable ist die Hauptvariable ma[nvars(R)]=var(n); ma[n]=var(nvars(R)); map phi=R,ma; list fac=factorize(phi(p)); list re=phi(fac); return(re); } example { "EXAMPLE:"; echo = 2; system("random",992851144); ring r=32003,(x,y,z,w,t),lp; poly p=y2w9+yz7t-yz5w4-z2w4t4-w8t3; factorize(p); //fast system("random",992851262); //factorize(p); //slow system("random",992851262); factorH(p); }