title: Efficient Matlab implementation of Multinomial Coefficient
body: <p>I want to calculate the multinomial coefficient:</p>
  
  <p><img src="https://i.stack.imgur.com/tbVq5.gif" alt="enter image description here"></p>
  
  <p>where it is satisifed <code>n=n0+n1+n2</code></p>
  
  <p>The Matlab implementation of this operator can be easily done in the function:</p>
  
  <pre><code>function N = nchooseks(k1,k2,k3)
      N = factorial(k1+k2+k3)/(factorial(k1)*factorial(k2)*factorial(k3)); 
  end
  </code></pre>
  
  <p>However, when the index is larger than 170, the factorial would be infinite which would generate <code>NaN</code> in some cases, e.g. <code>180!/(175! 3! 2!) -&gt; Inf/Inf-&gt; NaN</code>.</p>
  
  <p>In other posts, they have solved this overflow issue for <a href="https://stackoverflow.com/questions/8591633/permutation-with-repetition-avoiding-overflow/8596462#8596462">C</a> and <a href="https://stackoverflow.com/questions/3037113/python-calculate-multinomial-probability-density-functions-on-large-dataset">Python</a>. </p>
  
  <ul>
  <li>In the case of C: <em>"you can make lists out of all the factorials,  then find the prime factorization of all the numbers in the lists, then cancel all the numbers on the top with those on the bottom, until the numbers are completely reduced"</em>.</li>
  <li>In the case of Python: <em>"make use of the fact that factorial(n) = gamma(n+1), use the logarithm of the gamma function and use additions instead of multiplications, subtractions instead of divisions"</em>.</li>
  </ul>
  
  <p>The first solution seems extremely slow, so I have tried the second option:</p>
  
  <pre><code>function N = nchooseks(k1,k2,k3)
      N = 10^(log_gamma(k1+k2+k3)-(log_gamma(k1)+log_gamma(k2)+log_gamma(k3))); 
  end
  function y = log_gamma(x),  y = log10(gamma(x+1));  end
  </code></pre>
  
  <p>I compare the original and log_gamma implementation with the following code:</p>
  
  <pre><code>% Calculate
  N=100; err = zeros(N,N,N);
  for n1=1:N,
      for n2=1:N,
          for n3=1:N,
              N1 = factorial(n1+n2+n3)/(factorial(n1)*factorial(n2)*factorial(n3)); 
              N2 = 10^(log10(gamma(n1+n2+n3+1))-(log10(gamma(n1+1))+log10(gamma(n2+1))+log10(gamma(n3+1)))); 
              err(n1,n2,n3) = abs(N1-N2); 
          end
      end
  end
  % Plot histogram of errors
  err_ = err(~isnan(err));
  [nelements,centers] = hist(err_(:),1e2);
  figure; bar(centers,nelements./numel(err_(:)));
  </code></pre>
  
  <p>However, the results are slightly different for some cases, as presented in the following histogram.  </p>
  
  <p><img src="https://i.stack.imgur.com/3r73a.png" alt="enter image description here"></p>
  
  <p>Thus, should I assume that my implementation is correct or the numerical error does not justify the number divergence?</p>