我想计算多项系数:
它满足qazxsw poi
该函数的Matlab实现可以在函数中轻松完成:
n=n0+n1+n2
然而,当指数大于170时,因子将是无穷大的,这将在某些情况下产生function N = nchooseks(k1,k2,k3)
N = factorial(k1+k2+k3)/(factorial(k1)*factorial(k2)*factorial(k3));
end
,例如NaN
。
在其他帖子中,他们已经为180!/(175! 3! 2!) -> Inf/Inf-> NaN
和C解决了这个溢出问题。
第一个解决方案似乎非常慢,所以我尝试了第二个选项:
Python
我将原始和log_gamma实现与以下代码进行比较:
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
但是,对于某些情况,结果略有不同,如下面的直方图所示。
因此,我应该假设我的实现是正确的还是数值误差不能证明数字偏差是正确的?
为什么不用这个?它速度快,不会溢出:
% 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_(:)));
很抱歉复活旧帖子,但对于未来的搜索者,你几乎可以肯定只是将你的多项式系数写成二项式系数的乘积,并使用内置方法计算二项式系数(或使用Pascal三角形或其他方式编写自己的系数)方法)。相关公式出现在N = prod([1:n]./[1:n0 1:n1 1:n2]);
的第一段。 (我在这里写,但似乎没有办法渲染LaTeX。)
这种方法的另一个好处是它可以获得溢出,因为因子都是整数。计算多项式系数时,没有内在的需要划分。
使用@jemidiah提供的提示,
Wikipedia section on multinomial coefficients
这是代码
和一些用法示例:
function c = multicoeff (k),
c = 1;
for i=1:length(k),
c = c* bincoeff(sum(k(1:i)),k(i));
end;
end