非常可靠地计算二项式系数

Question

在 C++ 中计算二项式系数的最佳方法是什么？我见过一些代码片段，但在我看来，它总是只在某些特定区域可行。我需要一个非常、非常、非常可靠的计算。我用伽玛函数试了一下：

unsigned n=N;
unsigned k=2;
number = tgammal(n + 1) / (tgammal(k + 1) * tgammal(n - k + 1));

但它已经在 n=8, k=2 of 1 时有所不同（并且在 n=30, k=2 时它崩溃了）.. 我“只”需要计算至少 n=3000，k=2。

Answer 1

这个

constexpr inline size_t binom(size_t n, size_t k) noexcept
{
    return
      (        k> n  )? 0 :          // out of range
      (k==0 || k==n  )? 1 :          // edge
      (k==1 || k==n-1)? n :          // first
      (     k+k < n  )?              // recursive:
      (binom(n-1,k-1) * n)/k :       //  path to k=1   is faster
      (binom(n-1,k) * n)/(n-k);      //  path to k=n-1 is faster
}

需要 O(min{k,n-k}) 操作，可靠且可以在编译时完成。

解释 二项式系数定义为B(n,k)=k!(n-k)!/n!，从中我们可以看到B(n,k)=B(n,n-k)。我们还可以得到递归关系n*B(n,k)=(n-k)*B(n-1,k)=k*B(n-1,k-1)。此外，k=0,1,n,n-1 的结果是微不足道的。

对于k=2，结果也是平凡的(n*(n-1))/2。

当然，您也可以使用其他整数类型来做到这一点。如果你需要知道一个超过最大可表示整数类型的二项式系数，你必须使用近似方法：使用double。在这种情况下，最好使用 gamma 函数

#include <cmath>
inline double logbinom(double n, double k) noexcept
{
    return std::lgamma(n+1)-std::lgamma(n-k+1)-std::lgamma(k+1);
}
inline double binom(double n, double k) noexcept
{
    return std::exp(logbinom(n,k));
}

Answer 2

您可以使用渐近更有效的循环公式：

constexpr inline size_t binom(size_t n, size_t k) noexcept
{
    return
      (        k> n  )? 0 :          // out of range
      (k==0 || k==n  )? 1 :          // edge
      (k==1 || k==n-1)? n :          // first
      binom(n - 1, k - 1) * n / k;   // recursive
}

这将仅使用 O(k) 运算来计算 C(n, k)。