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 处停止......

4

4 回答 4

1

要获得 smallk和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!)近似值n。这完全取决于您想要的准确程度。

于 2016-09-30T20:21:17.947 回答
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 时），但如果您需要所有系数并且它具有更大范围的可接受输入，它会快得多。

于 2016-09-30T23:21:43.860 回答
1

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

  1. 用大整数实现公式（例如 long long）
  2. 尽快尝试分裂操作（如卓然建议）
  3. 添加代码以检查每个除法和乘法的正确性
  4. 解决不正确的除法或乘法，例如
    • 尝试卓然提出的循环除法，如果失败则返回初始算法（在den中累积除数的乘积）
    • 将未解析的乘数、除数存储在额外的长整数中，并尝试在下一次迭代循环中解析它们
  5. 如果你真的使用大数，那么你的结果可能不适合长整数。那么在这种情况下，您可以切换到 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;
}
于 2016-09-30T21:45:10.927 回答
0

这可能不是 OP 所要寻找的，但是对于具有二元熵函数的大 n，我们可以分析近似 nCr。它在中提到

于 2016-09-30T20:22:44.887 回答