function Mp=M_to_Mp(M,index1,index2,beta)
%computes recombination step on marginal frequencies M(1:2^L)
%indices are shifted up by one so that P_x=P(1+x), M_x=M(1+x)
M=[1;M]; %add zero index
N=length(M); %N=2^L
Mp=zeros(N,1); Mp(1)=1;
%bug for L=16 if don't double, 1+(2^L - 1) = 2^L-1 in unint16
for x=1:N-1 %index of M
     Mp(1+x) = (beta{x}.*M(1+double(index1{x})))' * M(1+double(index2{x})); %sum of product
end
Mp=Mp(2:end); %remove zero index