有效地计算组合和排列
Posted
技术标签:
【中文标题】有效地计算组合和排列【英文标题】：counting combinations and permutations efficiently 【发布时间】：2011-01-06 23:35:30 【问题描述】：我有一些代码来计算排列和组合，我正在努力让它更好地处理大数。
我找到了一种更好的置换算法，可以避免大量中间结果，但我仍然认为我可以在组合方面做得更好。
到目前为止，我已经放入了一个特殊情况来反映 nCr 的对称性，但我仍然想找到一个更好的算法来避免调用 factorial(r)，这是一个不必要的大中间结果。如果没有这种优化，最后一个 doctest 会花费很长时间来计算阶乘 (99000)。
谁能建议一种更有效的组合计数方法？
from math import factorial
def product(iterable):
prod = 1
for n in iterable:
prod *= n
return prod
def npr(n, r):
"""
Calculate the number of ordered permutations of r items taken from a
population of size n.
>>> npr(3, 2)
6
>>> npr(100, 20)
1303995018204712451095685346159820800000
"""
assert 0 <= r <= n
return product(range(n - r + 1, n + 1))
def ncr(n, r):
"""
Calculate the number of unordered combinations of r items taken from a
population of size n.
>>> ncr(3, 2)
3
>>> ncr(100, 20)
535983370403809682970
>>> ncr(100000, 1000) == ncr(100000, 99000)
True
"""
assert 0 <= r <= n
if r > n // 2:
r = n - r
return npr(n, r) // factorial(r)
【问题讨论】：
很久以前就有人问过这个问题，但无论如何...我设计了一个计算 C(n, m) = n 的算法！ / (m! (n-m)!) 只要结果适合整数（很容易成为长整数）。我是用 Java 编写的，但应该很容易将其翻译成 Python 或任何其他程序语言：***.com/questions/50292530/…（查找combinations(int n, int m)）
【参考方案1】：
如果 n 离 r 不远，那么使用组合的递归定义可能会更好，因为 xC0 == 1 你只会有几次迭代：
这里相关的递归定义是：
nCr = (n-1)C(r-1) * n/r
这可以很好地使用尾递归和以下列表进行计算：
[(n - r, 0), (n - r + 1, 1), (n - r + 2, 2), ..., (n - 1, r - 1), (n, r )]
izip(xrange(n - r + 1, n+1), xrange(1, r+1))izip(xrange(n - r + 1, n+1), xrange(1, r+1))当然很容易在 Python 中生成（我们省略了第一个条目，因为 nC0 = 1）请注意，这假设 r 也可以优化使用
现在我们只需要使用带有 reduce 的尾递归来应用递归步骤。我们从 1 开始，因为 nC0 为 1，然后将当前值与列表中的下一个条目相乘，如下所示。
from itertools import izip
reduce(lambda x, y: x * y[0] / y[1], izip(xrange(n - r + 1, n+1), xrange(1, r+1)), 1)
【讨论】：
对于单个 nCr 这更好，但是当您有多个 nCr（大约 N）时，动态编程方法更好，即使它有很长的设置时间，因为它不会t 除非必要，否则溢出到“bignum”中。【参考方案2】：两个相当简单的建议：
为避免溢出，请在日志空间中执行所有操作。使用 log(a * b) = log(a) + log(b) 和 log(a / b) = log(a) - log(b) 的事实。这使得处理非常大的阶乘变得容易：log(n! / m!) = log(n!) - log(m!) 等。
使用 gamma 函数代替阶乘。您可以在scipy.stats.loggamma 中找到一个。这是一种比直接求和更有效的计算对数阶乘的方法。 loggamma(n) == log(factorial(n - 1))，类似地，gamma(n) == factorial(n - 1)。
【讨论】：
在日志空间做事的好建议。不过，不确定您所说的“精确”是什么意思。使用 log-floats 不会导致大数舍入错误吗？ @Gorgapor：我想更清楚的表述方式是：“避免溢出”。已编辑。 请注意，由于浮点数的精度有限，这不会给出准确的结果。 @starblue：但是你知道真正的答案必须是一个整数，所以如果你做类似 round(exp(logFactorial(n))) 的事情，它对于小 n 来说是精确的。对于大 n 它可能是不精确的，但除了（慢）任意精度之外的任何东西都将是完全错误的。 对于小 n 计算这个并不麻烦。关键是要对大 n 进行准确计算，而且我已经在使用任意精度，因为我使用的是 python longs。【参考方案3】：scipy 中有一个尚未提及的函数：scipy.special.comb。根据您的 doctest 的一些快速计时结果，这似乎很有效（comb(100000, 1000, 1) == comb(100000, 99000, 1) 约为 0.004 秒）。
[虽然这个特定问题似乎与算法有关，但问题 is there a math ncr function in python 被标记为与此重复...]
【讨论】：
【参考方案4】：如果您不需要纯 python 解决方案，gmpy2 可能会有所帮助（gmpy2.comb 非常快）。
【讨论】：
感谢您的参考，这是一个非常实用的解决方案。不过，这对我来说更像是一个学习项目，所以我对算法比实际结果更感兴趣。 对于那些在写完这个答案几年后得出这个答案的人，gmpy 现在被称为 gmpy2。【参考方案5】：from scipy import misc
misc.comb(n, k)
应该允许你计算组合
【讨论】：
【参考方案6】：如果您正在计算 N 选择 K（我认为您正在使用 ncr 进行操作），那么有一个动态编程解决方案可能会快很多。这将避免阶乘，而且您可以保留该表以供以后使用。
这是它的教学链接：
http://www.csc.liv.ac.uk/~ped/teachadmin/algor/dyprog.html
不过，抱歉，我不确定如何更好地解决您的第一个问题。
编辑：这是模型。有一些非常搞笑的错误，所以它当然可以经受更多的清理。
import sys
n = int(sys.argv[1])+2#100
k = int(sys.argv[2])+1#20
table = [[0]*(n+2)]*(n+2)
for i in range(1,n):
table[i][i] = 1
for i in range(1,n):
for j in range(1,n-i):
x = i+j
if j == 1: table[x][j] = 1
else: table[x][j] = table[x-1][j-1] + table[x-1][j]
print table[n][k]
【讨论】：
据我所知，这个实现似乎是 O(n^2) 而我布置的尾递归是 O(n)。 似乎使用了不同的递归定义。这里n选k=n-1选k-1+n-1选k，而我用n选k=n-1选k-1 * n/k 确实如此。我将很快编辑这篇文章，以包含该算法的快速 Python 模型。你的速度明显更快。我会把我的帖子留在这里，以防 Gorgapor 有一些奇特的机器，其中乘法需要几个小时。 >.> 这可能是 O(N^2)，但它会预先计算 nCr 的所有组合对，所以如果你要大量使用 nCr 和很多不同的值，这会更快，因为查找是O(1) 并且不太容易溢出。不过，对于一个值，O(N) 算法更好。【参考方案7】：更有效的 nCr 解决方案 - 空间方面和精度方面。
保证中间值 (res) 始终为 int 且永远不会大于结果。空间复杂度为 O(1)（无列表、无拉链、无堆栈），时间复杂度为 O(r) - 正好是 r 次乘法和 r 次除法。
def ncr(n, r):
r = min(r, n-r)
if r == 0: return 1
res = 1
for k in range(1,r+1):
res = res*(n-k+1)/k
return res
【讨论】：
【参考方案8】：如果您的问题不需要知道排列或组合的确切数量，那么您可以使用Stirling's approximation 作为阶乘。
这将导致这样的代码：
import math
def stirling(n):
# http://en.wikipedia.org/wiki/Stirling%27s_approximation
return math.sqrt(2*math.pi*n)*(n/math.e)**n
def npr(n,r):
return (stirling(n)/stirling(n-r) if n>20 else
math.factorial(n)/math.factorial(n-r))
def ncr(n,r):
return (stirling(n)/stirling(r)/stirling(n-r) if n>20 else
math.factorial(n)/math.factorial(r)/math.factorial(n-r))
print(npr(3,2))
# 6
print(npr(100,20))
# 1.30426670868e+39
print(ncr(3,2))
# 3
print(ncr(100,20))
# 5.38333246453e+20
【讨论】：
阶乘的主要问题是结果的大小，而不是计算它的时间。此外，这里的结果值比浮点值准确表示的值要大得多。【参考方案9】：3.7 之前的 Python：
def prod(items, start=1):
for item in items:
start *= item
return start
def perm(n, k):
if not 0 <= k <= n:
raise ValueError(
'Values must be non-negative and n >= k in perm(n, k)')
else:
return prod(range(n - k + 1, n + 1))
def comb(n, k):
if not 0 <= k <= n:
raise ValueError(
'Values must be non-negative and n >= k in comb(n, k)')
else:
k = k if k < n - k else n - k
return prod(range(n - k + 1, n + 1)) // math.factorial(k)
对于 Python 3.8+：
math.perm()
math.comb()
有趣的是，一些手动实现的组合功能可能比math.comb()快：
def math_comb(n, k):
return math.comb(n, k)
def comb_perm(n, k):
k = k if k < n - k else n - k
return math.perm(n, k) // math.factorial(k)
def comb(n, k):
k = k if k < n - k else n - k
return prod(range(n - k + 1, n + 1)) // math.factorial(k)
def comb_other(n, k):
k = k if k > n - k else n - k
return prod(range(n - k + 1, n + 1)) // math.factorial(k)
def comb_reduce(n, k):
k = k if k < n - k else n - k
return functools.reduce(
lambda x, y: x * y[0] // y[1],
zip(range(n - k + 1, n + 1), range(1, k + 1)),
1)
def comb_iter(n, k):
k = k if k < n - k else n - k
result = 1
for i in range(1, k + 1):
result = result * (n - i + 1) // i
return result
def comb_iterdiv(n, k):
k = k if k < n - k else n - k
result = divider = 1
for i in range(1, k + 1):
result *= (n - i + 1)
divider *= i
return result // divider
def comb_fact(n, k):
k = k if k < n - k else n - k
return math.factorial(n) // math.factorial(n - k) // math.factorial(k)
所以实际上comb_perm()（用math.perm() 和math.factorial() 实现）实际上在大多数情况下都比math.comb() 快。
请注意，comb_reduce() 很慢，与@wich's answer 的方法基本相同，而同样相对较慢的comb_iter() 与@ZXX's answer 基本相同。
【讨论】：
【参考方案10】：from numpy import prod
def nCr(n,r):
numerator = range(n, max(n-r,r),-1)
denominator = range(1, min(n-r,r) +1,1)
return int(prod(numerator)/prod(denominator))
【讨论】：
【参考方案11】：使用xrange() 而不是range() 会稍微加快速度，因为没有中间列表被创建、填充、迭代和销毁。另外，reduce() 和 operator.mul。
【讨论】：
对不起，我不清楚，我的代码是 python 3，而不是 python 2。python 3 中的范围与 python 2 中的 xrange 相同。【参考方案12】：对于 N 选择 K，您可以使用帕斯卡三角形。基本上，您需要保留大小为 N 的数组来计算所有 N 选择 K 值。只需添加即可。
【讨论】：
这基本上是 Agor 建议的，但它会是 O(n^2)。由于现在使用乘法和除法实际上不再是问题，因此使用不同的递归关系可以使算法像我描述的那样 O(n)。【参考方案13】：您可以输入两个整数并导入数学库以找到阶乘，然后应用 nCr 公式
import math
n,r=[int(_)for _ in raw_input().split()]
f=math.factorial
print f(n)/f(r)/f(n-r)
【讨论】：
以上是关于有效地计算组合和排列的主要内容，如果未能解决你的问题，请参考以下文章