多项式系数的高效Matlab实现

问题描述 投票:2回答:3

我想计算多项系数:

它满足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-> NaNC解决了这个溢出问题。

  • 在C的情况下:“你可以从所有因子中做出列表,然后找到列表中所有数字的素数因子分解,然后用底部的数字取消顶部的所有数字,直到数字完全减少”。
  • 在Python的情况下:“利用factorial(n)= gamma(n + 1)这一事实,使用gamma函数的对数并使用加法而不是乘法,减法而不是除法”。

第一个解决方案似乎非常慢,所以我尝试了第二个选项:

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

但是,对于某些情况,结果略有不同,如下面的直方图所示。

因此,我应该假设我的实现是正确的还是数值误差不能证明数字偏差是正确的?

matlab overflow factorial multinomial
3个回答
3
投票

为什么不用这个?它速度快,不会溢出:

% 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_(:)));

2
投票

很抱歉复活旧帖子,但对于未来的搜索者,你几乎可以肯定只是将你的多项式系数写成二项式系数的乘积,并使用内置方法计算二项式系数(或使用Pascal三角形或其他方式编写自己的系数)方法)。相关公式出现在N = prod([1:n]./[1:n0 1:n1 1:n2]); 的第一段。 (我在这里写,但似乎没有办法渲染LaTeX。)

这种方法的另一个好处是它可以获得溢出,因为因子都是整数。计算多项式系数时,没有内在的需要划分。


0
投票

使用@jemidiah提供的提示,

Wikipedia section on multinomial coefficients

这是代码

enter image description here

和一些用法示例:

function c = multicoeff (k), 
    c = 1; 
    for i=1:length(k), 
      c = c* bincoeff(sum(k(1:i)),k(i)); 
    end; 
end
© www.soinside.com 2019 - 2024. All rights reserved.