有效地计算组合和排列

Question

我有一些代码来计算排列和组合，我正在努力使它更好地适用于大数。

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

到目前为止，我已经提出了一个特殊的例子来反映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)

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

在scipy中有一个函数还没有提到：scipy.special.comb。基于您的doctest的一些快速计时结果（约0.004秒用于comb(100000, 1000, 1) == comb(100000, 99000, 1)），它似乎是有效的。

[虽然这个特定的问题似乎是关于算法的，但是问题is there a math ncr function in python被标记为这个的副本…]

Answer 3

两个相当简单的建议：

1. 为避免溢出，请在日志空间中执行所有操作。使用以下事实：log（a*b）=log（a）+log（b），log（a/b）=log（a）-log（b）。这使得使用非常大的阶乘很容易：log（n！/m！）=对数（n！）-对数（m！），等等。

2. 使用gamma函数而不是阶乘。你可以在scipy.stats.loggamma中找到一个。这是一种比直接求和更有效的计算对数因子的方法。loggamma(n) == log(factorial(n - 1))，类似地，gamma(n) == factorial(n - 1)。