逼近非常大数的二项式系数的对数

Question

我目前正在努力计算非常大的数字的二项式系数，假设 n

到目前为止，我已经尝试了许多方法来处理这些计算产生的大量数字。但是，问题是我不需要计算这些二项式系数一次而是数十万次。这意味着计算阶乘的常规方法迄今为止过于昂贵，而像long long int 这样的标准数据类型太有限，无法容纳这些数字。

我已经尝试过来自Boost 库的多精度数据类型，但正如我所提到的，计算如此多次会导致性能极慢。我也尝试过使用OpenMP 的多线程，但性能优势仍然太低。结果，我改用计算二项式系数的对数来保持数字很小。虽然这解决了大数的问题，但这并没有加快进程。这就是为什么我尝试了对数二项式系数的斯特林近似。我当前的解决方案如下所示：

#include <math.h>

long double calc_hgeom(unsigned int k, unsigned int n, unsigned int K, unsigned int N)
{
    long double hprob = std::exp((log_C(K, k) + log_C(N-K, n-k)) - log_C(N, n));
    return hprob;
}

long double log_C(unsigned int u, unsigned int m)
{
    long double C = u * std::log(u) - m * std::log(m) - (u-m) * std::log(u-m)) + 0.5 * (std::log(u) - std::log(m) - std::log(u-m) - std::log(2*M_PI));
    return C;
}

但是，结果与实际值相差很大，最高可达 7 %。因此我的问题是：是否有一种有效的方法来计算二项式系数的对数，或者我的近似值是否可以改进以提高准确性？

任何帮助将不胜感激，因为这个计算是我整个算法的基础。

Answer 1

考虑 R 的 lchoose 函数 ...

> choose(10000, 5000) 
[1] Inf
> lchoose(10000, 5000)
[1] 6926.641

基础 R 语言源代码库是解决此类问题的重要思路。

见https://github.com/wch/r-source/blob/trunk/src/nmath/choose.c

这里的技巧是使用 ln 转换的输入来避免溢出。

请注意代码是在 GNU 许可下的。

Answer 2

您应该使用Sterling's approximation formula 表示n！ ，它应用于二项式系数，为您提供：

对于二项式系数本身，对于对数，只需取该等式右侧的对数即可；大多数这些东西很快就会变得简单得多。你仍然有k！但是，对于较大的 k，您将再次需要近似公式。最终你会得到更可行的东西（即在数值上更稳定）。

如果这还不够好，即如果您仍然有几乎相互抵消的条款，请考虑在其中一个变量上应用 Taylor expansion。