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