问二项式系数舍入误差
EN

Stack Overflow用户

提问于 2016-09-30 12:04:11

回答 4查看 440关注 0票数 0

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

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; 
}

与递归公式相比，这段代码在单个核心线程上确实非常快速，尽管后者不太容易出现舍入错误，因为它只涉及和而不涉及除法。因此，我想测试这些标志的伟大价值，并试图评估500选择250 (订单10^160)。我发现“相对误差”小于10^(-19)，所以基本上是相同的数字，尽管它们有10^141的不同之处。

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

我谷歌了一些表的二项式系数，以便我可以复制，但我找到的最好的一个在n=100.

rounding

binomial-coefficients

回答 4

Stack Overflow用户

回答已采纳

发布于 2016-09-30 15:21:43

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

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

用libm的Gnu实现进行了测试，结果一致地在几个ULP内精确值，并且通常比基于乘法和除法的解决方案更好。

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

如果n太大，以至于tgammal(n+1)超过了长双(大于1754)的范围，但不超过分子溢出的范围，那么没有大号库就能得到最好的乘法解。但是，您也可以使用

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时)，但是如果你需要所有的系数，并且它有更大范围的可接受的输入，它会快得多。

票数 1

Stack Overflow用户

发布于 2016-09-30 12:21:17

要获得小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。这对组合学也有帮助。当ln(n!)较大时，Stirling公式是n的一个很好的近似。这完全取决于你想要的准确程度。

票数 1

Stack Overflow用户

发布于 2016-09-30 13:45:10

如果您真的想使用乘法公式，我建议采用一种基于异常的方法。

用大整数实现公式(例如长长)
尽快尝试分工合作(如卓然所建议)
添加代码以检查每个除法和乘法的正确性。
解决不正确的除法或乘法，例如
- 尝试卓然提出的循环除法，但如果失败，则返回初始算法(在den中累加除数乘积)。
- 将未解析的乘子、除数存储在附加的长整数中，并尝试在下一个迭代循环中求解它们。

如果您真的使用大数，那么您的结果可能不适合长整数。然后，在这种情况下，您可以切换到长双或使用您的个人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;
}