135

math我需要在 Python 中计算组合(nCr) ,但在numpystat 库中找不到执行此操作的函数。类似于以下类型的函数:

comb = calculate_combinations(n, r)

我需要可能组合的数量,而不是实际组合,所以itertools.combinations我不感兴趣。

最后,我想避免使用阶乘,因为我要计算组合的数字可能会变得太大,而阶乘会变得很可怕。

这似乎是一个非常容易回答的问题,但是我被关于生成所有实际组合的问题所淹没,这不是我想要的。

4

19 回答 19

135

请参阅scipy.special.comb(旧版本 scipy 中的 scipy.misc.comb)。当exact为 False 时,它​​使用 gammaln 函数来获得良好的精度,而无需花费太多时间。在确切的情况下,它返回一个任意精度的整数,这可能需要很长时间来计算。

于 2010-06-11T18:29:58.130 回答
128

为什么不自己写呢?这是一个单线或这样的:

from operator import mul    # or mul=lambda x,y:x*y
from fractions import Fraction

def nCk(n,k): 
  return int( reduce(mul, (Fraction(n-i, i+1) for i in range(k)), 1) )

测试 - 打印帕斯卡三角形:

>>> for n in range(17):
...     print ' '.join('%5d'%nCk(n,k) for k in range(n+1)).center(100)
...     
                                                   1                                                
                                                1     1                                             
                                             1     2     1                                          
                                          1     3     3     1                                       
                                       1     4     6     4     1                                    
                                    1     5    10    10     5     1                                 
                                 1     6    15    20    15     6     1                              
                              1     7    21    35    35    21     7     1                           
                           1     8    28    56    70    56    28     8     1                        
                        1     9    36    84   126   126    84    36     9     1                     
                     1    10    45   120   210   252   210   120    45    10     1                  
                  1    11    55   165   330   462   462   330   165    55    11     1               
               1    12    66   220   495   792   924   792   495   220    66    12     1            
            1    13    78   286   715  1287  1716  1716  1287   715   286    78    13     1         
         1    14    91   364  1001  2002  3003  3432  3003  2002  1001   364    91    14     1      
      1    15   105   455  1365  3003  5005  6435  6435  5005  3003  1365   455   105    15     1   
    1    16   120   560  1820  4368  8008 11440 12870 11440  8008  4368  1820   560   120    16     1
>>> 

PS。编辑以替换int(round(reduce(mul, (float(n-i)/(i+1) for i in range(k)), 1)))int(reduce(mul, (Fraction(n-i, i+1) for i in range(k)), 1))大 N/K 不会出错

于 2010-06-12T01:25:40.530 回答
55

对谷歌代码的快速搜索给出(它使用来自@Mark Byers's answer的公式):

def choose(n, k):
    """
    A fast way to calculate binomial coefficients by Andrew Dalke (contrib).
    """
    if 0 <= k <= n:
        ntok = 1
        ktok = 1
        for t in xrange(1, min(k, n - k) + 1):
            ntok *= n
            ktok *= t
            n -= 1
        return ntok // ktok
    else:
        return 0

choose()scipy.misc.comb()比您需要准确答案时快 10 倍(在所有 0 <= (n,k) < 1e3 对上测试) 。

def comb(N,k): # from scipy.comb(), but MODIFIED!
    if (k > N) or (N < 0) or (k < 0):
        return 0L
    N,k = map(long,(N,k))
    top = N
    val = 1L
    while (top > (N-k)):
        val *= top
        top -= 1
    n = 1L
    while (n < k+1L):
        val /= n
        n += 1
    return val
于 2010-06-11T19:12:11.727 回答
45

如果您想要准确的结果速度,请尝试gmpy --gmpy.comb应该完全按照您的要求进行,而且速度非常快(当然,作为gmpy的原作者,我偏见;-)。

于 2010-06-11T21:15:40.350 回答
30

如果您想要准确的结果,请使用sympy.binomial. 这似乎是最快的方法,放下手。

x = 1000000
y = 234050

%timeit scipy.misc.comb(x, y, exact=True)
1 loops, best of 3: 1min 27s per loop

%timeit gmpy.comb(x, y)
1 loops, best of 3: 1.97 s per loop

%timeit int(sympy.binomial(x, y))
100000 loops, best of 3: 5.06 µs per loop
于 2013-11-23T11:00:15.927 回答
28

在很多情况下,数学定义的直译就足够了(记住 Python 会自动使用大数算术):

from math import factorial

def calculate_combinations(n, r):
    return factorial(n) // factorial(r) // factorial(n-r)

reduce对于我测试的一些输入(例如 n=1000 r=500),这比在另一个(当前投票最高的)答案中建议的一个衬里快 10 倍以上。另一方面,它的表现优于@JF Sebastian 提供的代码段。

于 2013-10-01T06:54:11.547 回答
22

从 开始Python 3.8,标准库现在包含math.comb计算二项式系数的函数:

数学梳子(n,k)

这是从 n 个项目中选择 k 个项目而不重复的方法数
n! / (k! (n - k)!)

import math
math.comb(10, 5) # 252
于 2019-06-01T10:27:43.967 回答
10

这是另一种选择。这个最初是用 C++ 编写的,因此可以将它向后移植到 C++ 以获取有限精度整数(例如 __int64)。优点是(1)它只涉及整数运算,以及(2)它通过连续成对的乘法和除法来避免使整数值膨胀。我用 Nas Banov 的帕斯卡三角形测试了结果,它得到了正确的答案:

def choose(n,r):
  """Computes n! / (r! (n-r)!) exactly. Returns a python long int."""
  assert n >= 0
  assert 0 <= r <= n

  c = 1L
  denom = 1
  for (num,denom) in zip(xrange(n,n-r,-1), xrange(1,r+1,1)):
    c = (c * num) // denom
  return c

基本原理:为了最小化乘法和除法的次数,我们将表达式重写为

    n!      n(n-1)...(n-r+1)
--------- = ----------------
 r!(n-r)!          r!

为尽可能避免乘法溢出,我们将按照以下 STRICT 顺序进行计算,从左到右:

n / 1 * (n-1) / 2 * (n-2) / 3 * ... * (n-r+1) / r

我们可以证明以这种顺序进行的整数运算是精确的(即没有舍入误差)。

于 2012-08-24T19:29:45.030 回答
6

您可以编写 2 个简单的函数,这些函数实际上比使用scipy.special.comb快 5-8 倍。实际上,您不需要导入任何额外的包,并且该功能非常易于阅读。诀窍是使用记忆来存储先前计算的值,并使用nCr的定义

# create a memoization dictionary
memo = {}
def factorial(n):
    """
    Calculate the factorial of an input using memoization
    :param n: int
    :rtype value: int
    """
    if n in [1,0]:
        return 1
    if n in memo:
        return memo[n]
    value = n*factorial(n-1)
    memo[n] = value
    return value

def ncr(n, k):
    """
    Choose k elements from a set of n elements - n must be larger than or equal to k
    :param n: int
    :param k: int
    :rtype: int
    """
    return factorial(n)/(factorial(k)*factorial(n-k))

如果我们比较时间

from scipy.special import comb
%timeit comb(100,48)
>>> 100000 loops, best of 3: 6.78 µs per loop

%timeit ncr(100,48)
>>> 1000000 loops, best of 3: 1.39 µs per loop
于 2018-10-25T10:40:27.380 回答
5

使用动态规划,时间复杂度为 Θ(n*m),空间复杂度为 Θ(m):

def binomial(n, k):
""" (int, int) -> int

         | c(n-1, k-1) + c(n-1, k), if 0 < k < n
c(n,k) = | 1                      , if n = k
         | 1                      , if k = 0

Precondition: n > k

>>> binomial(9, 2)
36
"""

c = [0] * (n + 1)
c[0] = 1
for i in range(1, n + 1):
    c[i] = 1
    j = i - 1
    while j > 0:
        c[j] += c[j - 1]
        j -= 1

return c[k]
于 2014-09-12T17:30:10.053 回答
5

如果你的程序有一个上限n(比如n <= N) 并且需要重复计算 nCr (最好是 >>N次),使用lru_cache可以给你带来巨大的性能提升:

from functools import lru_cache

@lru_cache(maxsize=None)
def nCr(n, r):
    return 1 if r == 0 or r == n else nCr(n - 1, r - 1) + nCr(n - 1, r)

构建缓存(隐式完成)需要花费大量O(N^2)时间。任何后续调用都nCr将返回O(1).

于 2017-03-30T03:05:57.160 回答
4

sympy 很容易。

import sympy

comb = sympy.binomial(n, r)
于 2016-07-17T18:21:07.607 回答
3

当 n 大于 20 时,直接公式产生大整数。

所以,又是一个回应:

from math import factorial

reduce(long.__mul__, range(n-r+1, n+1), 1L) // factorial(r)

简短、准确和高效,因为这通过坚持 long 来避免 python 大整数。

与 scipy.special.comb 相比,它更准确、更快:

 >>> from scipy.special import comb
 >>> nCr = lambda n,r: reduce(long.__mul__, range(n-r+1, n+1), 1L) // factorial(r)
 >>> comb(128,20)
 1.1965669823265365e+23
 >>> nCr(128,20)
 119656698232656998274400L  # accurate, no loss
 >>> from timeit import timeit
 >>> timeit(lambda: comb(n,r))
 8.231969118118286
 >>> timeit(lambda: nCr(128, 20))
 3.885951042175293
于 2015-12-04T11:46:38.177 回答
3

仅使用随 Python 分发的标准库

import itertools

def nCk(n, k):
    return len(list(itertools.combinations(range(n), k)))
于 2016-11-10T22:19:51.477 回答
3

这个功能非常优化。

def nCk(n,k):
    m=0
    if k==0:
        m=1
    if k==1:
        m=n
    if k>=2:
        num,dem,op1,op2=1,1,k,n
        while(op1>=1):
            num*=op2
            dem*=op1
            op1-=1
            op2-=1
        m=num//dem
    return m
于 2019-05-13T21:43:15.790 回答
2

对于相当大的输入,这可能与您在纯 python 中执行的速度一样快:

def choose(n, k):
    if k == n: return 1
    if k > n: return 0
    d, q = max(k, n-k), min(k, n-k)
    num =  1
    for n in xrange(d+1, n+1): num *= n
    denom = 1
    for d in xrange(1, q+1): denom *= d
    return num / denom
于 2015-05-24T20:42:25.913 回答
2

这是适合您的有效算法

for i = 1.....r

   p = p * ( n - i ) / i

print(p)

例如 nCr(30,7) = fact(30) / ( fact(7) * fact(23)) = ( 30 * 29 * 28 * 27 * 26 * 25 * 24 ) / (1 * 2 * 3 * 4 * 5 * 6 * 7)

所以只需运行从 1 到 r 的循环就可以得到结果。


在蟒蛇中:

n,r=5,2
p=n
for i in range(1,r):
   p = p*(n - i)/i
else:
   p = p/(i+1)
print(p)
于 2019-10-05T10:49:41.343 回答
1

我从这个线程和这里链接的库中计时了 17 个不同的函数。

因为我觉得这里转储有点多,所以我把函数的代码放在了一个pastebin中。

我做的第一个测试是将帕斯卡三角形建立到第 100 行。我用 timeit 做了 100 次。下面的数字是构建一次三角形所需的平均时间(以秒为单位)。

gmpy2.gmpy2.comb 0.0012259269999998423
math.comb 0.007063110999999935
__main__.stdfactorial2 0.011469491
__main__.scipybinom 0.0120114319999999
__main__.stdfactorial 0.012105122
__main__.scipycombexact 0.012569045999999844
__main__.andrewdalke 0.01825201100000015
__main__.rabih 0.018472497000000202
__main__.kta 0.019374668000000383
__main__.wirawan 0.029312811000000067
scipy.special._basic.comb 0.03221609299999954
__main__.jfsmodifiedscipy 0.04332894699999997
__main__.rojas 0.04395155400000021
sympy.functions.combinatorial.factorials.binomial 0.3233529779999998
__main__.nasbanov 0.593365528
__main__.pantelis300 1.7780402499999999

您可能会注意到这里只有 16 个函数。那是因为该recursive()函数甚至无法在合理的时间内完成一次,所以我不得不将它从 timeit 测试中排除。严重的是,它已经持续了几个小时。

我还计时了并非所有上述功能都支持的各种其他类型的输入。请记住,我只运行了每 10 次测试,因为 nCr 的计算量很大而且我很不耐烦

n 的小数值

__main__.scipybinom 0.011481370000000001
__main__.kta 0.01869513999999999
sympy.functions.combinatorial.factorials.binomial 6.33897291

r 的小数值

__main__.scipybinom 0.010960040000000504
scipy.special._basic.comb 0.03681254999999908
sympy.functions.combinatorial.factorials.binomial 3.2962564499999987

n 和 r 的小数值

__main__.scipybinom 0.008623409999998444
sympy.functions.combinatorial.factorials.binomial 3.690936439999999

n 的负值

gmpy2.gmpy2.comb 0.010770989999997482
__main__.kta 0.02187850000000253
__main__.rojas 0.05104292999999984
__main__.nasbanov 0.6153183200000001
sympy.functions.combinatorial.factorials.binomial 3.0460310799999943

n 的负小数值,r 的小数值

sympy.functions.combinatorial.factorials.binomial 3.7689941699999965

目前实现最大速度和多功能性的最佳解决方案是混合函数,可根据输入在不同算法之间进行选择

def hybrid(n: typing.Union[int, float], k: typing.Union[int, float]) -> typing.Union[int, float]:
    # my own custom hybrid solution
    def is_integer(n):
        return isinstance(n, int) or n.is_integer()
    if k < 0:
        raise ValueError("k cannot be negative.")
    elif n == 0:
        return 0
    elif k == 0 or k == n:
        return 1
    elif is_integer(n) and is_integer(k):
        return int(gmpy2.comb(int(n), int(k)))
    elif n > 0:
        return scipy.special.binom(n, k)
    else:
        return float(sympy.binomial(n, k))

由于sympy.binomial()速度如此之慢,真正理想的解决方案是组合scipy.special.binom()对分数执行良好的代码和gmpy2.comb()对整数执行良好的代码。scipy 的 funcgympy2 的 func都是用我不太熟悉的 C 语言编写的。

于 2021-12-14T04:49:23.447 回答
0

这是使用内置记忆装饰器的@killerT2333 代码。

from functools import lru_cache

@lru_cache()
def factorial(n):
    """
    Calculate the factorial of an input using memoization
    :param n: int
    :rtype value: int
    """
    return 1 if n in (1, 0) else n * factorial(n-1)

@lru_cache()
def ncr(n, k):
    """
    Choose k elements from a set of n elements,
    n must be greater than or equal to k.
    :param n: int
    :param k: int
    :rtype: int
    """
    return factorial(n) / (factorial(k) * factorial(n - k))

print(ncr(6, 3))
于 2018-12-31T08:46:38.343 回答