function [index1 index2 beta]=precompute(L, alpha)
%precompute indices for M->Mp computation
%precompute sums over alpha and store in beta

%ymax takes values from 0 to x*, where x* is equal to 
%x but with the highest significant bit set to zero
ymax=zeros(2^L-1,1,'uint16'); 
for l=0:L-1  
    for j=0:2^l-1
        ymax(2^l+j)=j;
    end
end

index1=cell(2^L-1,1); index2=cell(2^L-1,1);
for x=uint16(1:2^L-1) %index of M
    y=0:ymax(x); %possible range of indices of first M in product
    index1{x}=y(logical(~bitand(y,x-y))); %indices of first M in product
    index2{x}=x-index1{x}; %indices of second M in product
end

beta=cell(2^L-1,1);
alpha=reshape(alpha,2*ones(1,L));
for x=uint16(1:2^L-1) %index of M
    beta{x}=alpha;
    for bit=1:L
        if ~bitget(x,bit) %if bit index of M' is 0
            beta{x}=sum(beta{x},bit); %sum over index for which bit in M' is 0
        end
    end
    beta{x}=beta{x}(:);
    beta{x}=beta{x}(1:end/2);
end