有效地计算组合和排列

Question

我有一些代码可以计算排列和组合,我正在努力让它更适合大数字.

我已经找到了一个更好的排列算法,避免了大的中间结果,但我仍然认为我可以做更好的组合.

到目前为止,我已经提出了一个特殊情况来反映nCr的对称性,但我仍然希望找到一种更好的算法来避免调用阶乘(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)

Answer 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)]

这当然很容易在Python中生成(我们省略了自nC0 = 1以来的第一个条目)izip(xrange(n - r + 1, n+1), xrange(1, r+1))注意这假定r <= n你需要检查它并交换它们,如果它们不是.如果r <n/2则优化使用,则r = n - 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)

Answer 2

两个相当简单的建议:

1. 为避免溢出,请在日志空间中执行所有操作.使用该记录的事实(A*B)=日志的(a)+日志(b)中,和log(A/B)=日志(一) - 日志(b)中.这使得使用非常大的因子很容易:log(n!/ m!)= log(n!) - log(m!)等.

2. 使用伽玛函数代替阶乘.你可以找到一个scipy.stats.loggamma.这是一种比直接求和更有效的计算对数因子的方法. loggamma(n) == log(factorial(n - 1)),同样地,gamma(n) == factorial(n - 1).

Answer 3

对于 Python 3.7 之前的版本：

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()n=256kk = n // 2

请注意，它相当慢，本质上与@wich 的答案comb_reduce()相同，而也相对较慢，本质上与@ZXX 的答案相同。comb_iter()

此处进行部分分析（没有comb_math()，comb_perm()因为截至上次编辑，Colab 的 Python 版本 - 3.7 - 不支持它们）。

Answer 4

scipy中有一个功能尚未提及:scipy.special.comb.基于你的doctest的一些快速计时结果(~0.004秒comb(100000, 1000, 1) == comb(100000, 99000, 1)),它看起来很有效.

[虽然这个特定问题似乎与算法有关,但问题是python中的数学ncr函数被标记为此复制...]

Answer 5

如果你不需要纯python解决方案,gmpy2可能会有所帮助(gmpy2.comb非常快).

Answer 6

from scipy import misc
misc.comb(n, k)

应该允许您计算组合

Answer 7

如果您正在计算 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]

Answer 8

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