我需要计算一个数字的组合。
计算 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 最有效的方法是什么？
问问题
28965 次
11 回答
17
您可以使用动态编程来生成二项式系数
您可以创建一个数组，然后使用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 的最大值
于 2012-06-14T12:16:08.367 回答
16
如果您需要为所有 n 计算它们，Ribtoks 的答案可能是最好的。对于单个 n，您最好这样做：
C[0] = 1
for (int k = 0; k < n; ++ k)
C[k+1] = (C[k] * (n-k)) / (k+1)
如果在乘法之后完成，则除法是精确的。
并注意 C[k] * (nk) 溢出：使用足够大的整数。
于 2012-06-14T12:30:54.450 回答
8
如果您想要对较大的 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
- 从缓冲区中读出系数 - 第一个数字是中心峰，其中头=尾，随后的条目对应于距离中心较远的对称对
准确性
在累积的浮点舍入误差使您无法获得准确的系数整数值之前，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 时间并损害准确性。
于 2013-03-25T21:18:09.833 回答
7
这是我的版本：
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
于 2012-06-14T14:25:36.867 回答
7
我最近写了一段代码，需要调用大约 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))
于 2013-03-20T19:44:19.333 回答
6
如果你真的只需要 n 远大于 p 的情况，一种方法是使用斯特林公式来计算阶乘。（如果 n>>1 和 p 是一阶，斯特林近似 n！和 (np)！，保持 p！原样等）
于 2012-06-14T20:21:12.427 回答
3
我自己的基准测试中最快的合理近似值是 Apache Commons Maths 库使用的近似值：http: //commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/special/Gamma.html #logGamma（双）
我和我的同事试图看看我们是否能打败它，同时使用精确的计算而不是近似值。除了一种方法慢了 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
真正糟糕的方法是使用尾递归实现帕斯卡三角形的各种尝试。
于 2014-08-11T14:48:53.550 回答
2
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;
}
于 2015-01-01T05:11:10.193 回答
1
我一直在寻找相同的东西但找不到它，所以我自己写了一个似乎对于最终结果适合长的任何二项式系数的最佳选择。
// 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;
}
于 2012-10-16T10:40:45.913 回答
1
时间复杂度：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));
}
}
于 2019-03-06T09:58:42.180 回答
1
如果我理解问题中的符号，您不仅想要 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片段：
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
于 2017-08-24T20:05:26.907 回答