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

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 时它崩溃了）。我 "only" 需要计算至少 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 代替。在这种情况下，最好使用伽玛函数

#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)。

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

Calculate binomial coffeficient very reliably

c++

math

binomial-coefficients

gamma-function