高效的Matlab实现多项式系数

Question

我想要计算多项式系数：

在满足条件 n=n0+n1+n2 的情况下。

该运算符的Matlab实现可以轻松地在以下函数中完成：

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。

在其他帖子中，他们解决了C和Python的溢出问题。

对于C语言情况： "您可以将所有阶乘制成列表，然后找到列表中所有数字的质因数分解，然后取消顶部与底部的所有数字，直到数字完全减少"。

对于Python情况："利用阶乘(n) = gamma(n+1)这个事实，使用伽玛函数的对数，使用加法而不是乘法，使用减法而不是除法"。

第一个解决方案似乎非常慢，因此我尝试了第二个选项：

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实现：

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

然而，对于某些情况，结果略有不同，如下直方图所示。

因此，我应该假设我的实现是正确的，还是数值误差不足以证明数字发散？

- tashuhka

---

计算 gamma 的对数，您可以直接使用 gammaln。 - Luis Mendo

相关链接：https://dev59.com/UWUp5IYBdhLWcg3wY29V#15302394 - tashuhka

Answer 1

为什么不使用这个呢？它速度快且不会出现溢出问题：

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

- Luis Mendo

---

我尝试了你的巧妙解决方案，得到了类似的错误分布。因此，我认为这三种实现是相似的。然而，你的解决方案是最快的（阶乘：57秒，log_gamma：15秒，乘积：6秒）。因此，我选择了你的解决方案。谢谢 :) - tashuhka

Answer 2

抱歉，重新激活旧帖子，但对于未来的搜索者，您几乎肯定只需将多项式系数写成二项式系数的乘积，并使用内置方法计算二项式系数（或编写自己的代码，可以使用Pascal三角形或其他方法）。有关公式出现在维基百科关于多项式系数的部分的第一段中。这种方法的另一个好处是，它是关于溢出最好的方法，因为所有因子都是整数。在计算多项式系数时没有固有的需要进行除法运算。

- Joshua P. Swanson

Answer 3

使用@jemidiah提供的提示，

这里是代码

function c = multicoeff (k), 
    c = 1; 
    for i=1:length(k), 
      c = c* bincoeff(sum(k(1:i)),k(i)); 
    end; 
end

以及一些使用示例：

octave:88> multicoeff([2 2 2])
ans =  90
octave:89> factorial(6)/(factorial(2)*factorial(2)*factorial(2))
ans =  90
octave:90> multicoeff([5 4 3])
ans =  27720
octave:91> factorial(12)/(factorial(5)*factorial(4)*factorial(3))
ans =  27720

- LEo

Answer 4

另一种方法是使用Yannis Manolopoulos的迭代方法。 假设我们有一个带有多项式条目的向量k。

function N = multicoeff (k),
  n=sum(k); 
  [_,imax]=max(k); 
  num=[n:-1:n-k(imax)-1]; 
  den=[]; k(imax)=[]; 
  for i=1:length(k), den=[den 1:k(i)]; endfor; 
  N=prod(num./den);
endfunction

例子

octave:2> k = [5 4 3];
octave:3> multicoeff (k)
ans =  27720

参考文献： Yannis Manolopoulos. 二项式系数计算。ACM SIGCSE通报，34（4）：65，2002年12月。doi：10.1145/820127.820168。URL https: //doi.org/10.1145/820127.820168。

- LEo