function partitionskofset(n,k)
%PARTITIONSKOFSET 产生[n]的所有的k分拆,个数即为第二类Stirling数S(n,k)
% partitionskofset( n, k ) 输入正整数n和k,程序输出所有[n]的k分拆。输出
%格式:一个长度为n的行向量p,p(i)属于{1,2,...,k},i=1,2,...,n. 如果p(i)=j, 表
%示数i属于第j部分。
%$Author: WBC$ $Date: 2005/10/19$
if k==1 % 一种特殊情况
disp(ones(1,n));
else
ss=ksubsets(n-1,k-1); %先获得[n-1]的所有k-1子集
c=zeros(1,n); %初始化c
c(1)=1; %总保证c的首项为1
for l=1:size(ss,1) %针对[n-1]的每一个k-1子集可得到相应的[n]的k分拆
c(2:end)=0; %把除首项外的其余项全部重新赋值为0
c(ss(l,:)+1)=1;%把[n-1]的k-1子集的每个元素加1,然后在c对应位置赋1,实际
%上这时c表示[n]的一个总包含元素1的k子集
partitionsofcom(c);%罗列[n]的该k子集相应的所有的[n]的k分拆
end
end
%--------------------------------------------------------------------------
%一个子函数
%--------------------------------------------------------------------------
function partitionsofcom(c)
%罗列[n]的k子集相应的所有的[n]的k分拆,这里c是首项为1的0-1向量。
n=length(c); %求出n,c的元素个数即为n
p=zeros(1,n); %初始化p
p(1)=1;%第一个元素总是让它属于第1部分
i=2; %用于标记p的当前位置的循环变量
while c(i)==0 %一直循环到c的第二个1出现时停止,因为k>=2,所以该循环一定会停止
p(i)=1; %这些0对应位置所代表的元素都属于第1部分
i=i+1; %循环变量加1
end %在循环停止时有c(i)==1
pzero=zeros(1,n); %记录第二个1右边的0在c中的位置
len=0; %表示到当前位置止0的个数(第二个1右边的0)
numone=zeros(1,n);%记录相应于pzero位置处左边的1的个数,与pzero同长度,皆为len
s1=2; %标记自左至右到当前位置止1的个数,到当前位置处有两个1
p(i)=s1; %该位置所代表的元素属于第2部分
i=i+1;
while i<=n
if c(i)>0
s1=s1+1;
p(i)=s1;
else
len=len+1;
pzero(len)=i;
numone(len)=s1;
end
i=i+1;
end
p(pzero(1:len))=1;
disp(p);
while 1
i=1;
while i<=len&&p(pzero(i))>=numone(i)
p(pzero(i))=1;
i=i+1;
end
if i>len
break;
end
p(pzero(i))=p(pzero(i))+1;
disp(p);
end