function y = multcoef(n,x)
%#codegen

% MULTCOEF Multinomial coefficient.
%  Y = MULTCOEF(N,X) returns the multinomial coefficient with parameter N at the values in X.
%  Let {X1, X2,  . . .  , Xk}, k > 1, be a set of random variables, each of which can take
%  the values 0, 1,  . . .  , n; such that for every set of k nonnegative integers {n1,
%  . . .  , nk} whose sum is n, the multinomial coefficient is,
%
%                                                        (n1 + n2 + ... nk)!
%        C(n; n1, n2, ..., nk) = (n1, n2, ... nk)! =  -------------------------  .
%                                                     n1! × n2! ×  . . .  × nk!
%
%  It is possible to work with large factorials.
%
%  Syntax: y = multcoef(n,x)
%
%  Inputs:
%       n - number of trials.
%       x - vector of the interested values.
%  Outputs:
%       y - multinomial coefficient.
%
%  Example. Assume that a die is thrown 60 times (n=60) and a record is kept of the number of
%  times a 1, 2, 3, 4, 5, or 6 is observed. Each trial (e.g., throw of a die) has a (1 or 2 or
%  3 or . . . or 6) mutually exclusive outcome. In the die tossing data, k = 6, we are
%  interested to get the multinomial coefficient.
%
%  Calling on Matlab the function:
%             multcoef(n,x)
%
%  where n=60 and x=[13,10,8,10,12,7];
%
%  Answer is:
%
%  ans
%      = 1.0425e+042
%

assert(n>=0)
assert(all(x>=0))
assert(sum(x) == n)
if n==0
    y=1;
else
    
    lfactor1 = sum(log(1:n));
    
    c = length(x);
    
    f = zeros(c,1);
    for i = 1:c,
        f(i) = sum(log(1:x(i)));
    end
    
    lfactor2 = sum(f);
    
    y = exp(lfactor1-lfactor2);
    y = round(y);
end
end
