多项式系数的高效Matlab实现 [英] Efficient Matlab implementation of Multinomial Coefficient

问题描述

我要计算多项式系数:

满足的地方n=n0+n1+n2

此运算符的Matlab实现可通过以下函数轻松完成:

The Matlab implementation of this operator can be easily done in the function:

function N = nchooseks(k1,k2,k3)
    N = factorial(k1+k2+k3)/(factorial(k1)*factorial(k2)*factorial(k3)); 
end

但是，当索引大于170时，阶乘将是无限的，这在某些情况下会生成NaN，例如180!/(175! 3! 2!) -> Inf/Inf-> NaN.

However, when the index is larger than 170, the factorial would be infinite which would generate NaN in some cases, e.g. 180!/(175! 3! 2!) -> Inf/Inf-> NaN.

在其他帖子中，他们已经解决了 C 和 Python .

In other posts, they have solved this overflow issue for C and Python.

• 对于C:"，您可以从所有阶乘中创建列表，然后找到列表中所有数字的素数分解，然后取消顶部的所有数字和底部的所有数字，直到数量完全减少"..
• 在Python中:利用factorial(n)= gamma(n + 1)的事实，使用gamma函数的对数，并使用加法而不是乘法，减法而不是除法" .

• In the case of C: "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".
• In the case of Python: "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".

第一种解决方案似乎非常慢，因此我尝试了第二种选择:

The first solution seems extremely slow, so I have tried the 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

我将原始和log_gamma实现与以下代码进行比较:

I compare the original and log_gamma implementation with the following 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_(:)));

但是，某些情况下的结果略有不同，如以下直方图所示.

However, the results are slightly different for some cases, as presented in the following histogram.

因此，我应该假设我的实现是正确的，还是数字错误不能证明数字的差异?

Thus, should I assume that my implementation is correct or the numerical error does not justify the number divergence?

推荐答案

为什么不使用它?快速，不会溢出:

Why not use this? It's fast and doesn't suffer from overflow:

N = prod([1:n]./[1:n0 1:n1 1:n2]);

这篇关于多项式系数的高效Matlab实现的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！