//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.59 2009-01-14 16:07:04 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
"
{
if (system("with","Namespaces"))
{
list @marie=names(Top);
}
else
{
list @marie=names();
}
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 ((system("with","Namespaces")) && 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 if id.
Zero generators of ideal/module are deleted
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)
}
///////////////////////////////////////////////////////////////////////////////
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);
}