生成二项式系数的最快方法答案

【问题标题】：Fastest way to generate binomial coefficients生成二项式系数的最快方法
【发布时间】：2012-06-14 12:10:23
【问题描述】：

我需要计算一个数字的组合。

计算 nCp where n>>p 的最快方法是什么？

我需要一种快速的方法来为多项式方程生成二项式系数，并且我需要获取所有项的系数并将其存储在一个数组中。

(a+b)^n = a^n + nC1 a^(n-1) * b + nC2 a^(n-2) * ...... +nC(n-1) a * b^(n-1) + b^n

计算 nCp 最有效的方法是什么？

【问题讨论】：

我尝试使用传统方法计算组合，但您可能已经猜到它的效率非常低，即使对于小型多项式也需要很多时间
@rao_555：n 和 p 的大小数字是多少？以及什么样的语言？例如，在对小数进行运算从根本上更快的语言中，是否值得尝试计算unsigned long long中的结果（示例取自C）并避免溢出，这对答案产生了影响通过巧妙的技巧，或者结果是否比这大得多，所以尝试没有意义。
@rao_555：好吧，如果 n 是 10^10 并且 p 是 n/2，那么 nCp 大约有 450 亿个十进制数字。就个人而言，我会尝试在那时更改要求。
抱歉，100 亿位数。仍然。
@SteveJessop：你必须接受过物理或工程方面的培训。

标签： algorithm math data-structures

【解决方案1】：

您可以使用动态规划来生成二项式系数

您可以创建一个数组，然后使用 O(N^2) 循环来填充它

C[n, k] = C[n-1, k-1] + C[n-1, k];

在哪里

C[1, 1] = C[n, n] = 1

之后，在您的程序中，您只需在 [n, k] 索引处查看二维数组即可获得 C(n, k) 值

更新之类的

for (int k = 1; k <= K; k++) C[0][k] = 0;
for (int n = 0; n <= N; n++) C[n][0] = 1;

for (int n = 1; n <= N; n++)
   for (int k = 1; k <= K; k++)
      C[n][k] = C[n-1][k-1] + C[n-1][k];

其中 N, K - 你的 n, k 的最大值

【讨论】：

N和K不应该是同一个值吗？我认为记住不成比例的帕斯卡三角形没有意义
@ScottNguyen 代码以最通用的形式显示。您可以根据需要修改 N 和 K
+1。很好的解决方案。只是为了具体而不是使用O（n * 2）内存，它可以使用O（N）内存进行评估，因为根据问题中的问题定义，它只对N的系数感兴趣。
@Dreamer，对于你和 Scott 指出的 OP 提出的问题，是 O(n) 内存，但仍然是 O(n^2) 时间。

【解决方案2】：

如果您需要为所有 n 计算它们，Ribtoks 的答案可能是最好的。对于单个 n，您最好这样做：

C[0] = 1
for (int k = 0; k < n; ++ k)
    C[k+1] = (C[k] * (n-k)) / (k+1)

除法是精确的，如果在乘法之后完成。

并注意 C[k] * (n-k) 溢出：使用足够大的整数。

【讨论】：

【解决方案3】：

如果您想对较大的 n 值进行完全扩展，FFT 卷积可能是最快的方法。在具有相等系数（例如一系列公平的硬币投掷）和偶数顺序（例如投掷次数）的二项式展开的情况下，您可以利用对称性：

理论

用表达式 A + A*cos(Pi*n/N) 表示两次抛硬币的结果（例如正面和反面总数之差的一半）。 N 是缓冲区中的样本数 - 偶数阶 O 的二项式展开将具有 O+1 个系数，并且需要 N >= O/2 + 1 个样本的缓冲区 - n 是正在生成的样本数，A 是比例因子通常为 2（用于生成二项式系数）或 0.5（用于生成二项式概率分布）。

请注意，在频率上，此表达式类似于这两次抛硬币的二项式分布 - 在对应于数字 (heads-tails)/2 的位置存在三个对称尖峰。由于对独立事件的整体概率分布进行建模需要对它们的分布进行卷积，因此我们希望在频域中对我们的表达式进行卷积，这相当于在时域中进行乘法运算。

换句话说，通过将我们的两次抛掷结果的余弦表达式提升为一次幂（例如，模拟 500 次抛掷，将其提升为 250 次幂，因为它已经代表一对），我们可以安排二项式分布大量出现在频域中。既然这都是真实的，那么我们可以用 DCT-I 代替 DFT 来提高效率。

算法

确定缓冲区大小 N，它至少为 O/2 + 1 并且可以方便地进行 DCT 处理
使用表达式 pow(A + A*cos(Pi*n/N),O/2) 对其进行初始化
应用正向 DCT-I
从缓冲区中读出系数 - 第一个数字是中心峰，其中 head=tails，随后的条目对应于远离中心的对称对

准确度

在累积的浮点舍入误差使您无法获得准确的系数整数值之前，O 可以达到多高是有限制的，但我猜这个数字相当高。双精度浮点可以完全准确地表示 53 位整数，我将忽略使用 pow() 所涉及的舍入损失，因为生成表达式将发生在 FP 寄存器中，给我们额外的 11尾数位以吸收英特尔平台上的舍入误差。因此，假设我们使用通过 FFT 实现的 1024 点 DCT-I，这意味着在转换过程中舍入误差会损失 10 位的精度，仅此而已，剩下大约 43 位的清晰表示。我不知道二项式展开的什么顺序会产生这种大小的系数，但我敢说它足以满足您的需求。

不对称扩展

如果您想要 a 和 b 的不等系数的不对称展开，您将需要使用双边（复数）DFT 和复数 pow() 函数。生成表达式 A*A*e^(-Pi*i*n/N) + A*B + B*B*e^(+Pi*i*n/N) [使用复 pow() 函数来提高它的一半扩展命令的力量]和DFT它。同样，缓冲区中的中心点是偏移量为零的中心点（但如果 A 和 B 非常不同，则不是最大值），然后是分布的上半部分。缓冲区的上半部分将包含分布的下半部分，对应于负值的正面-反面值。

请注意，源数据是 Hermitian 对称的（输入缓冲区的后半部分是前半部分的复共轭），因此该算法不是最优的，可以使用所需一半的复数到复数 FFT 来执行尺寸以获得最佳效率。

毋庸置疑，与上述对称分布的纯实算法相比，所有复数幂运算都会占用更多 CPU 时间并损害准确性。

【讨论】：

【解决方案4】：

这是我的版本：

def binomial(n, k):
if k == 0:
    return 1
elif 2*k > n:
    return binomial(n,n-k)
else:
    e = n-k+1
    for i in range(2,k+1):
        e *= (n-k+i)
        e /= i
    return e

【讨论】：

【解决方案5】：

我最近编写了一段代码，需要调用大约 1000 万次二进制系数。所以我做了一个组合查找表/计算方法，仍然不会太浪费内存。您可能会发现它很有用（而且我的代码在公共领域）。代码在

http://www.etceterology.com/fast-binomial-coefficients

有人建议我在这里内联代码。一个大喇叭查找表似乎很浪费，所以这是最后一个函数，以及一个生成表的 Python 脚本：

extern long long bctable[]; /* See below */

long long binomial(int n, int k) {
    int i;
    long long b;
    assert(n >= 0 && k >= 0);

    if (0 == k || n == k) return 1LL;
    if (k > n) return 0LL;

    if (k > (n - k)) k = n - k;
    if (1 == k) return (long long)n;

    if (n <= 54 && k <= 54) {
        return bctable[(((n - 3) * (n - 3)) >> 2) + (k - 2)];
    }
    /* Last resort: actually calculate */
    b = 1LL;
    for (i = 1; i <= k; ++i) {
        b *= (n - (k - i));
        if (b < 0) return -1LL; /* Overflow */
        b /= i;
    }
    return b;
}

#!/usr/bin/env python3

import sys

class App(object):
    def __init__(self, max):
        self.table = [[0 for k in range(max + 1)] for n in range(max + 1)]
        self.max = max

    def build(self):
        for n in range(self.max + 1):
            for k in range(self.max + 1):
                if k == 0:      b = 1
                elif  k > n:    b = 0
                elif k == n:    b = 1
                elif k == 1:    b = n
                elif k > n-k:   b = self.table[n][n-k]
                else:
                    b = self.table[n-1][k] + self.table[n-1][k-1]
                self.table[n][k] = b

    def output(self, val):
        if val > 2**63: val = -1
        text = " {0}LL,".format(val)

        if self.column + len(text) > 76:
            print("\n   ", end = "")
            self.column = 3
        print(text, end = "")
        self.column += len(text)

    def dump(self):
        count = 0
        print("long long bctable[] = {", end="");

        self.column = 999
        for n in range(self.max + 1):
            for k in range(self.max + 1):
                if n < 4 or k < 2 or k > n-k:
                    continue
                self.output(self.table[n][k])
                count += 1
        print("\n}}; /* {0} Entries */".format(count));

    def run(self):
        self.build()
        self.dump()
        return 0

def main(args):
    return App(54).run()

if __name__ == "__main__":
    sys.exit(main(sys.argv))

【讨论】：

【解决方案6】：

如果您真的只需要 n 远大于 p 的情况，一种方法是使用 Stirling's formula 进行阶乘。（如果 n>>1 和 p 是一阶，Stirling 近似 n! 和 (n-p)!，保持 p! 原样等）

【讨论】：

更好的是，使用 lgamma 和添加来避免溢出

【解决方案7】：

我自己的基准测试中最快的合理近似值是 Apache Commons Math 库使用的近似值：http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/special/Gamma.html#logGamma(double)

我和我的同事试图看看我们是否可以击败它，同时使用精确计算而不是近似值。除了一种方法慢了 2-3 倍之外，所有方法都失败了（许多命令慢了）。性能最佳的方法使用https://math.stackexchange.com/a/202559/123948，这是代码（在 Scala 中）：

var i: Int = 0
var binCoeff: Double = 1
while (i < k) {
  binCoeff *= (n - i) / (k - i).toDouble
  i += 1
}
binCoeff

各种尝试使用尾递归实现帕斯卡三角的非常糟糕的方法。

【讨论】：

【解决方案8】：

nCp = n! / ( p! (n-p)! ) =
      ( n * (n-1) * (n-2) * ... * (n - p) * (n - p - 1) * ... * 1 ) /
      ( p * (p-1) * ... * 1     * (n - p) * (n - p - 1) * ... * 1 )

如果我们删减分子和分母的相同项，我们只需要最少的乘法。我们可以用 C 语言编写一个函数来执行 2p 乘法和 1 次除法得到 nCp：

int binom ( int p, int n ) {
    if ( p == 0 ) return 1;
    int num = n;
    int den = p;
    while ( p > 1 ) {
        p--;
        num *= n - p;
        den *= p;
    }
    return num / den;
}

【讨论】：

num 可能超过 32/64 位，即使结果可能不会。只有p <= 12，该功能实际上是安全的。

【解决方案9】：

我一直在寻找相同的东西但找不到它，所以我自己写了一个似乎对于任何最终结果适合 Long 的二项式系数来说都是最佳的。

// Calculate Binomial Coefficient 
// Jeroen B.P. Vuurens
public static long binomialCoefficient(int n, int k) {
    // take the lowest possible k to reduce computing using: n over k = n over (n-k)
    k = java.lang.Math.min( k, n - k );

    // holds the high number: fi. (1000 over 990) holds 991..1000
    long highnumber[] = new long[k];
    for (int i = 0; i < k; i++)
        highnumber[i] = n - i; // the high number first order is important
    // holds the dividers: fi. (1000 over 990) holds 2..10
    int dividers[] = new int[k - 1];
    for (int i = 0; i < k - 1; i++)
        dividers[i] = k - i;

    // for every dividers there is always exists a highnumber that can be divided by 
    // this, the number of highnumbers being a sequence that equals the number of 
    // dividers. Thus, the only trick needed is to divide in reverse order, so 
    // divide the highest divider first trying it on the highest highnumber first. 
    // That way you do not need to do any tricks with primes.
    for (int divider: dividers) {
       boolean eliminated = false;
       for (int i = 0; i < k; i++) {
          if (highnumber[i] % divider == 0) {
             highnumber[i] /= divider;
             eliminated = true;
             break;
          }
       }
       if(!eliminated) throw new Error(n+","+k+" divider="+divider);
    }


    // multiply remainder of highnumbers
    long result = 1;
    for (long high : highnumber)
       result *= high;
    return result;
}

【讨论】：

相当不错，但旨在避免浮点运算的嵌套 for 循环似乎比仅使用浮点运算更昂贵。
好点。但是，您是否尝试过浮点方法来处理非常大的数字？我想知道它是否在某些时候不会失去精度。这种方法仍然可以使用 BigDecimal 而不是 longs 来计算准确的结果。
这个算法不正确。循环for (int divider: dividers) ... 无效。
@Atom 您能否通过说明哪个输入的答案不正确来详细说明为什么该循环无效？我刚刚再次对其进行了测试，它位于 lib.MathTools 的 maven io.github.htools 中，并且似乎对我有用...
@JeroenVuurens 对于n 和k 的某些组合，永远不会到达break 语句，因此不会消除分隔符。

【解决方案10】：

如果我理解问题中的表示法，您不仅需要 nCp，您实际上还需要所有 nC1、nC2、... nC(n-1)。如果这是正确的，我们可以利用以下关系使这变得相当简单：

对于所有 k>0：nCk = prod_{from i=1..k}( (n-i+1)/i )
即对于所有 k>0：nCk = nC(k-1) * (n-k+1) / k

这是一个实现这种方法的python sn-p：

def binomial_coef_seq(n, k):
    """Returns a list of all binomial terms from choose(n,0) up to choose(n,k)"""
    b = [1]
    for i in range(1,k+1):
        b.append(b[-1] * (n-i+1)/i)
    return b

如果您需要直到某个 k > 上限（n/2）的所有系数，您可以使用对称性来减少需要执行的操作数量，方法是在上限（n/2）的系数处停止，然后仅回填只要你需要。

import numpy as np

def binomial_coef_seq2(n, k):
    """Returns a list of all binomial terms from choose(n,0) up to choose(n,k)"""

    k2 = int(np.ceiling(n/2))

    use_symmetry =  k > k2
    if use_symmetry:
        k = k2

    b = [1]
    for i in range(1, k+1):
        b.append(b[-1] * (n-i+1)/i)

    if use_symmetry:
        v = k2 - (n-k)
        b2 = b[-v:]
        b.extend(b2)
    return b

【讨论】：

【解决方案11】：

时间复杂度：O（分母）空间复杂度：O(1)

public class binomialCoeff {
    static double binomialcoeff(int numerator, int denominator) 
    { 
        double res = 1; 
        //invalid numbers
        if (denominator>numerator || denominator<0 || numerator<0) {
            res = -1;
            return res;}
        //default values
        if(denominator==numerator || denominator==0 || numerator==0)
            return res;


        // Since C(n, k) = C(n, n-k) 
        if ( denominator > (numerator - denominator) ) 
            denominator = numerator - denominator;


        // Calculate value of [n * (n-1) *---* (n-k+1)] / [k * (k-1) *----* 1] 
        while (denominator>=1)
        { 

        res *= numerator;
        res = res / denominator; 

        denominator--;
        numerator--;
        } 

        return res; 
    } 

    /* Driver program to test above function*/
    public static void main(String[] args) 
    { 
        int numerator = 120; 
        int denominator = 20; 
        System.out.println("Value of C("+ numerator + ", " + denominator+ ") "
                                + "is" + " "+ binomialcoeff(numerator, denominator)); 
    } 

}

【讨论】：