Raj*_*ula 25 algorithm math data-structures
我需要计算一个数字的组合.
什么是计算nCp的最快方法,其中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最有效的方法？
Rib*_*oks 15
您可以使用动态编程来生成二项式系数
您可以创建一个数组,然后使用O(N ^ 2)循环来填充它
C[n, k] = C[n-1, k-1] + C[n-1, k];
哪里
C[1, 1] = C[n, n] = 1
在你的程序中,你可以得到C(n,k)值只是在[n,k]索引处查看你的2D数组
像这样更新 smth
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的最大值
小智 13
如果你需要为所有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):使用足够大的整数.
如果想要对n的大值进行完全展开,则FFT卷积可能是最快的方法.在具有相等系数的二项式扩展(例如,一系列公平的硬币投掷)和偶数顺序(例如投掷数量)的情况下,您可以利用对称性:
理论
代表两个硬币投掷的结果(例如的头和尾的总数之间的差的一半)与表达式A + A*COS(PI*N/N).N是缓冲区中的样本数 - 偶数阶O的二项式扩展将具有O + 1个系数,并且需要N> = O/2 + 1个样本的缓冲区 - n是生成的样本编号,A是a比例因子通常为2(用于生成二项式系数)或0.5(用于生成二项式概率分布).
请注意,在频率上,此表达式类似于这两个硬币投掷的二项分布 - 在与数字(头尾)/ 2对应的位置处存在三个对称尖峰.由于造型独立事件的整体概率分布需要卷积它们的分布,我们要卷积在频域,这相当于在时域乘我们的表达.
换句话说,通过将两次投掷的结果的余弦表达式提高到一个幂(例如,为了模拟500次投掷,将其提升到250的幂,因为它已经代表一对),我们可以安排大的二项式分布要出现在频域中的数字.由于这一切都是真实的,我们可以用DCT-I代替DFT以提高效率.
算法
准确性
还有多高Ø可前滚存的浮点舍入误差抢你的系数精确的整数值的限制,但我猜的数量是相当高的.双精度浮点可以完全准确地表示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非常不同则不是最大值),然后是分布的上半部分.缓冲区的上半部分将包含分布的下半部分,对应于负数的head-minus-tails值.
请注意,源数据是Hermitian对称的(输入缓冲区的后半部分是第一个的复共轭),因此该算法不是最优的,可以使用所需大小的一半的复杂到复杂的FFT来执行,以获得最佳效果.效率.
毋庸置疑,与上面对称分布的纯实数算法相比,所有复指数都会咀嚼更多的CPU时间并损害精度.
这是我的版本:
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
我最近编写了一段代码,需要调用大约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))