2

我想计算多项式系数：

在此处输入图像描述

满意的地方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）这一事实，使用 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

我将原始和 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_(:)));

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

在此处输入图像描述

因此，我应该假设我的实现是正确的还是数字错误不能证明数字分歧是合理的？

4

4 回答 4

3

为什么不使用这个？它速度快，不会溢出：

N = prod([1:n]./[1:n0 1:n1 1:n2]);
于 2014-03-10T13:18:57.510 回答
2

抱歉恢复旧帖子，但对于未来的搜索者，您几乎可以肯定只将多项式系数写为二项式系数的乘积，并使用内置方法来计算二项式系数（或使用帕斯卡三角形或其他方法编写您自己的方法）。相关公式出现在关于多项式系数的维基百科部分的第一段中。（我会在这里写，但似乎没有办法渲染 LaTeX。）

这种方法的另一个好处是它尽可能好地处理溢出，因为因子都是整数。在计算多项式系数时，没有内在的需要除法。

于 2015-09-28T05:10:08.877 回答
0

使用@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
于 2019-04-08T14:28:01.687 回答
0

另一种方法是使用 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。网址https://doi.org/10.1145/820127.820168。

于 2020-02-28T13:53:59.223 回答