Search code examples
arraysmatlabloopsnested-loops

Compute all d-dimensional monomials of degree less than k


I want to compute all d-dimensional monomials of degree less than k and put them into a cell array - Pbase - in an ordered way (I use Matlab but this ploblem applies for other languages too). k and d are provided by the user but arbitrary.

What I have done so far,(the array called degreeindex show the location of monomials of a certain degree in Pbase):

n=nchoosek(d+k,k);%i know there are n over k possibilities for such monomials
Pbase=cell(1,n);
degreeindex=zeros(k+1,3);%showing [degree, start, end]
x=sym('x',[d,1]);

%the polynomials of degree 0
Pbase{1}=1;
degreeindex(1,1:3)=[0,1,1];%initialized for degree 0

%degree 1 monomials
degree=2;
degreeindex(degree,:)=[degree-1,degreeindex(degree-1,3)+1,degreeindex(degree-1,3)+1];
for i=1:d
    Pbase{degreeindex(degree,2)+i-1}=x(i);
    degreeindex(degree,3)=degreeindex(degree,3)+1;
end

%degree2 monomials
degree=3;
degreeindex(degree,:)=[degree-1,degreeindex(degree-1,3),degreeindex(degree-1,3)+1];
for i=1:d
    for j=i:d
        Pbase{degreeindex(degree,3)-1}=x(i).*x(j);
        degreeindex(degree,3)=degreeindex(degree,3)+1;
    end
end

%degree3 monomials
degree=4;
degreeindex(degree,:)=[degree-1,degreeindex(degree-1,3),degreeindex(degree-1,3)+1];
for i1=1:d
    for i2=i1:d
        for i3=i2:d
        Pbase{degreeindex(degree,3)-1}=x(i1).*x(i2)*x(i3);
        degreeindex(degree,3)=degreeindex(degree,3)+1;
        end
    end
end
...

The problem is that I don't find a way to do this for an arbitrary degree k. In the solution above I would have to include a new (and deeper) nested loop for every degree.

I know this seems like a faily trivial problem but I can't get my head around it. I appreciate all advice.


Solution

  • You are trying to do a type of partitioning, specifically you want to find all integer partitions of k into m or fewer parts (where order matters).

    Using nchoosek you can easily achive this (adapted from the Matlab File Exchange,here):

    d = 3; % num dims
    k = 4; % degree
    x=sym('x',[d,1]);
    m = nchoosek(k+d-1,d-1); 
    dividers = [zeros(m,1),nchoosek((1:(k+d-1))',d-1),ones(m,1)*(k+d)]; 
    a = diff(dividers,1,2)-1;
    PBase = cell(1, size(a,1));
    for i = 1:size(a,1)
        PBase{i} = prod(x.' .^ a(i,:));
    end
    

    Loop over whatever values of k you need.