二项式系数舍入误差

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 相比，此代码在单个核心线程上确实非常快，尽管后者较少受到舍入误差的影响，因为它只涉及求和而不涉及除法。所以我想测试这些算法的巨大价值，并尝试评估 500 选择 250（顺序 10^160）。我发现 "relative error" 小于 10^(-19)，所以基本上它们是相同的数字，尽管它们有所不同，例如 10^141.

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

我在谷歌上搜索了一些二项式系数表，所以我可以从中复制，但我发现最好的一个在 n=100 处停止...

Answer 1

要获得小 k 和 m 的精确整数结果，更好的解决方案可能是（您的代码略有变化）：

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

每次做完除法后得到一个组合数num/=i，所以你不会被截断。要获得更大 k 和 m 的近似结果，您的解决方案可能不错。但请注意，long double 乘法已经比整数的乘法和除法（unsigned long 或 size_t）慢得多。如果你想得到更大的数字，可能必须编码或从库中包含一个大整数 class。如果对极大整数 n 的 n! 有快速阶乘算法，您也可以 google。这也可能有助于组合数学。当 n 很大时，斯特林公式是 ln(n!) 的一个很好的近似值。这完全取决于您想要的准确度。

Answer 2

这可能不是 OP 正在寻找的东西，但可以通过二元熵函数对大 n 分析近似 nCr。

中提到

Answer 3

如果您真的想使用乘法公式，我会推荐一种基于异常的方法。

用大整数（例如long long）实现公式
尽快尝试除法运算（卓然建议）
添加代码以检查每个除法和乘法的正确性
解决不正确的除法或乘法，例如
- 尝试卓然提出的循环除法，如果失败则返回初始算法（累加den中除数的乘积）
- 将未解决的乘数、除数存储在额外的长整数中，并尝试在下一次迭代循环中解决它们
如果您真的使用大数字，那么您的结果可能不适合长整数。那么在这种情况下，您可以切换到 long double 或使用您的个人 LongInteger 存储。

这是一个框架代码，给你一个思路：

long long binomial_l(int k, int m)
{
    int i,j;
    long long num=1, den=1;
    j=m<(k-m)?m:(k-m);
    for(i=1;i<=j;i++)
    {
        int multiplier=(k+1-i);
        int divisor=i;
        long long candidate_num=num*multiplier;
        //check multiplication
        if((candidate_num/multiplier)!=num)
        {
            //resolve exception...
        }
        else
        {
            num=candidate_num;
        }

        candidate_num=num/divisor;
        //check division
        if((candidate_num*divisor)==num)
        {
            num=candidate_num;
        }
        else
        {
            //resolve exception
            den*=divisor;
            //this multiplication should also be checked...
        }
    }
    long long candidate_result= num/den; 
    if((candidate_result*den)==num)
    {
        return candidate_result;
    }
    // you should not get here if all exceptions are resolved
    return 0;
}

Answer 4

如果您只是计算 n 相当大但不超过 1750 的单个二项式系数 C(n,k)，那么使用像样的 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 时），但如果您需要所有系数并且它具有更大范围的可接受输入，它会更快。

二项式系数舍入误差

Binomial Coefficients rounding error

c

rounding

binomial-coefficients