overflow - Efficient Matlab implementation of Multinomial Coefficient -
i want calculate multinomial coefficient:
where satisifed n=n0+n1+n2
the matlab implementation of operator can done in function:
function n = nchooseks(k1,k2,k3) n = factorial(k1+k2+k3)/(factorial(k1)*factorial(k2)*factorial(k3)); end
however, when index larger 170, factorial infinite generate nan
in cases, e.g. 180!/(175! 3! 2!) -> inf/inf-> nan
.
in other posts, have solved overflow issue c , python.
- in case of c: "you can make lists out of factorials, find prime factorization of numbers in lists, cancel numbers on top on bottom, until numbers reduced".
- in case of python: "make use of fact factorial(n) = gamma(n+1), use logarithm of gamma function , use additions instead of multiplications, subtractions instead of divisions".
the first solution seems extremely slow, have tried second option:
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
i compare original , log_gamma implementation following code:
% calculate n=100; err = zeros(n,n,n); n1=1:n, n2=1:n, 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_(:)));
however, results different cases, presented in following histogram.
thus, should assume implementation correct or numerical error not justify number divergence?
why not use this? it's fast , doesn't suffer overflow:
n = prod([1:n]./[1:n0 1:n1 1:n2]);
Comments
Post a Comment