c - 二项式系数舍入误差

Question

我必须用 c 计算表达式的二项式系数 (x+y)**n，n 非常大(数量级为 500-1000)。我想到的第一个计算二项式系数的算法是 multiplicative formula .所以我把它编码到我的程序中

long double binomial(int k, int m)
{
    int i,j;
    long double num=1, den=1;
    j=m<(k-m)?m:(k-m);
    for(i=1;i<=j;i++)
    {
        num*=(k+1-i);
        den*=i;
    }
    return num/den; 
}

这段代码在单核线程上非常快，例如与 recursive formula 相比，尽管后者较少受到舍入误差的影响，因为它只涉及总和而不涉及除法。 所以我想测试这些算法的巨大值(value)，并尝试评估 500 选择 250(顺序 10^160)。我发现“相对误差”小于 10^(-19)，所以基本上它们是相同的数字，尽管它们相差 10^141。

所以我想知道:有没有一种方法可以评估计算错误的顺序？有没有比乘法公式更精确的快速计算二项式系数的方法？因为我不知道算法的精度，所以我不知道在哪里截断斯特林级数以获得更好的结果。

我在谷歌上搜索了一些二项式系数表，因此我可以从中复制，但我找到的最好的表在 n=100 时停止...

Answer 1

如果您只是计算单个二项式系数 C(n,k)，其中 n 相当大但不大于 1750，那么使用像样的 C 库的最佳选择是使用 tgammal标准库函数:

tgammal(n+1) / (tgammal(n-k+1) * tgammal(k+1))

使用 libm 的 Gnu 实现进行测试，结果始终在精确值的几个 ULP 范围内，并且通常优于基于乘法和除法的解决方案。

如果k足够小(或大)到二项式系数不会溢出64位精度，那么你可以通过交替乘除得到一个精确的结果。

如果 n 太大以至于 tgammal(n+1) 超出了 long double 的范围(超过 1754)但又没有大到分子溢出，那么在没有 bignum 库的情况下，乘法解决方案是最好的解决方案。但是，您也可以使用

expl(lgammal(n+1) - lgammal(n-k+1) - lgammal(k+1))

不太精确但更容易编码。 (此外，如果系数的对数对您有用，则上述公式适用于相当大的 n 和 k 范围。不必使用 expl 将提高准确性。)

如果您需要具有相同 n 值的一系列二项式系数，那么您最好的选择是迭代加法:

void binoms(unsigned n, long double* res) {
  // res must have (n+3)/2 elements
  res[0] = 1;
  for (unsigned i = 2, half = 0; i <= n; ++i) {
    res[half + 1] = res[half] * 2;
    for (int k = half; k > 0; --k)
      res[k] += res[k-1];
    if (i % 2 == 0)
      ++half;
  }
}

以上仅产生 k 从 0 到 n/2 的系数。它具有比乘法算法稍大的舍入误差(至少当 k 接近 n/2 时)，但如果您需要所有系数并且它具有更大范围的可接受输入，它会更快。