function [P]=multinomcoeff(m_vec)


% Computes the multinomial coefficient (N over m1 m2 ... mK) where
% N=sum m_k and where the m:s are found as elements in the vector m_vec. 
% The formula: P=N!/(m1!m2!...mK!). 
% Since computing factorials of large numbers deals with extremely large number, a straightforward
% implementation is unwise since it leads to overflow.
% I here use prime factorization of the factors appearing in the factorials
% and remove all factors in the numerator and denomiator that are identical.
% Tomas Olofsson, Jan 2011.

m_vec=m_vec(:);
K=length(m_vec);

N=sum(m_vec);

pr=primes(N); % A list of prime numbers smaller or equal to N. Prime numbers larger that that should not be necessary in these calculations.

pr=[1 pr];
Npr=length(pr);

%Nf=factor(N); % Factorize N

% Accumulate these prime numbers in the primenumber bin:
Nbin=zeros(Npr,1);
for n=1:N
    nf=factor(n);
    for l=1:length(nf)
        I=find(pr==nf(l),1);
        Nbin(I)=Nbin(I)+1;
    end
end

Mbin=zeros(Npr,K);
for k=1:K
    mk=m_vec(k);
    for k2=1:mk
        k2f=factor(k2);
        for l=1:length(k2f)
            I=find(pr==k2f(l),1);
            Mbin(I,k)=Mbin(I,k)+1;
        end
    end
end

% Remove common factors:
Nresult=Nbin;
for k=1:K
    Nresult=Nresult-Mbin(:,k);
end

% Finally, compute the product of all prime factors in Nresult
P=1;
for n=1:Npr
    P=P*pr(n)^Nresult(n);
end