The library realrad.lib contains a static proc subset which
computes a list of all subset of the range 1,...,n.
Each subset is again represented by a list of ints.
Since the proc is static, it cannot be accessed directly from an user
on the Top level but it can be invoked from another proc when
calling it with complete name as Packagename::procname.
Thus you can easily built a proc which returns a list of intvecs,
describing all k-element subsets out of 1,...n.
The code looks like:
Code:
proc choosesubsetsv1(int n,int k)
{
if(k<0 or k>n)
{
return(L);
}
// the list of all non empty subsets out of 1,...,n
list l = Realrad::subset(n);
list L,li;
int i,j;
for(i=1;i<=size(l);i++)
{
li = l[i];
if (size(li)==k) // select those subsets with k elements
{
j++;
L[j]= intvec(li[1..k]); // convert this list to intvec
}
}
return(L);
}
Note the following:
Due to an unnecessary restriction made by the author of the
library, you have to work in a ring of characteristic zero.
The proc subset is be renamed in Singular 3-1-1 to subsets
Code:
LIB "realrad.lib";
ring r0 = 0,x,dp;
> choosesubsetsv1(5,3);
[1]:
1,2,3
[2]:
1,2,4
[3]:
1,3,4
[4]:
2,3,4
[5]:
1,2,5
[6]:
1,3,5
[7]:
2,3,5
[8]:
1,4,5
[9]:
2,4,5
[10]:
3,4,5
Certainly, it is a huge overload to compute first
allm subsets and then to select those you are interested in.
This limits the effectiveness to approximatelty n=13.
There is another approach, closely related to polynomial
computations.
The idea is the following.
Define the linear polynomial f= x_1+... + x_n in n variables
and compute the kth power. Those terms which depends on
exactly k variables represent the subsets of the k-element
subsets -- their exponation vector consists only of 0's and 1's
where the 1's indicate the selected subset.
But this has the same performance as before, in fact this will compute again [b] all [b] subsets.
To improve it, you may do two things:
* if k > n-k, then compute compute only the f^(n-k) and
consider the complementary position of the 0s.
* instead of computing in a single one step f^k, select iteratively
the monomials depending on i variables and multiply again with f.
But it can be done much better. Instead of computing the kth
power of the polynomial x_1+...+x_n, compute the
kth power of the maximal ideal M =(x_1,...,x_n).
This computation uses a builtin command and is much faster
than the methods decribed before.
After computing M^k you have again select the monomials
depending on kvariables. To do this, you just throw away the
monomials which contain some squared or higher variable.
This is, you make Groebner basis computation to
reduce M^k w.r..t (x_1^2,...x_n^2).
The following proc implements these ideas and is much faster,
and not having the above mentioned restrictions.
Code:
proc choosesubsets(int n, int k)
"USAGE: choosesubsets(n,k); n, k int
RETURN: list, all subsets of 1,...,n with k elements
EXAMPLE: example choosesubsets; shows an example
"
{
list L;
int i,j,m;
intvec v,w;
int flipped = (n-k) < k;
if (flipped) { k = n-k;}
ring r= 0,(x(1..n)),dp;
ideal J;
for (i=1;i<=n;i++) {J[i] = x(i)^2;} // J is a std-basis
attrib(J,"isSB",1);
ideal K = reduce(maxideal(k),J)+0; //
for(j=1;j<=size(K);j++) // transform the exponents to the indices
{
v = leadexp(K[j][1]);
m = 0;
for(i=1;i<=n;i++)
{
if(v[i]==!flipped)
{
m++;
w[m]=i;
}
}
L[j] = w;
}
return(L);
}
C. Gorzel
Login
Register
Forum FAQ
Forum Search
