`scipy.misc.comb`是否比自行编写的二项式计算更快？

Question

现在确定了 scipy.misc.comb 比自定义实现更快吗？

根据一个旧答案 Statistics: combinations in Python，当计算组合数 nCr 时，这个自制函数比 scipy.misc.comb 更快:

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

但是在我的电脑上运行了一些测试之后，这似乎不是事实，使用以下脚本：

from scipy.misc import comb
import random, time

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

def timing(f):
    def wrap(*args):
        time1 = time.time()
        ret = f(*args)
        time2 = time.time()
        print '%s function took %0.3f ms' % (f.__name__, (time2-time1)*1000.0)
        return ret
    return wrap

@timing
def test_func(combination_func, nk):
    for n,k in nk:
        combination_func(n, k)

nk = []
for _ in range(1000):
    n = int(random.random() * 10000)
    k = random.randint(0,n)
    nk.append((n,k))

test_func(comb, nk)
test_func(choose, nk)

我得到以下输出：

$ python test.py
/usr/lib/python2.7/dist-packages/scipy/misc/common.py:295: RuntimeWarning: overflow encountered in exp
  vals = exp(lgam(N+1) - lgam(N-k+1) - lgam(k+1))
999
test_func function took 32.869 ms
999
test_func function took 1859.125 ms

$ python test.py
/usr/lib/python2.7/dist-packages/scipy/misc/common.py:295: RuntimeWarning: overflow encountered in exp
  vals = exp(lgam(N+1) - lgam(N-k+1) - lgam(k+1))
999
test_func function took 32.265 ms
999
test_func function took 1878.550 ms

“时间分析测试显示，新的scipy.misc.comb比临时的choose()函数更快吗？我的测试脚本有误导致计时不准确吗？

为什么scipy.misc.comb现在更快了？是因为使用了一些cython/c封装技巧吗？”

---

编辑过的

在@WarrenWeckesser评论后：

使用scipy.misc.comb()时，默认浮点近似值会导致计算出现浮点溢出。

（请参阅http://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.misc.comb.html获取文档）

当使用下面的函数并测试exact=True时，它会使用长整数而不是浮点数进行计算，但计算1000个组合时速度要慢得多：

@timing
def test_func(combination_func, nk):
    for i, (n,k) in enumerate(nk):
        combination_func(n, k, exact=True)

[输出]:

$ python test.py
test_func function took 3312.211 ms
test_func function took 1764.523 ms

$ python test.py
test_func function took 3320.198 ms
test_func function took 1782.280 ms

- alvas

---

默认情况下，scipy的comb函数计算的是浮点值，当参数足够大时会得到一个近似值。你应该使用comb函数中的参数exact=True来比较计算时间。 - Warren Weckesser

哇，使用 exact=True 后速度超慢。那么，有没有理由不使用临时函数而使用 scipy.misc.comb 呢？ - alvas

好问题！如果你感到有动力，你可以添加任何似乎与 https://github.com/scipy/scipy/issues/3449 相关的评论。 - Warren Weckesser

Answer 1

参考scipy.misc.comb的源代码，结果的更新例程如下：

    val = 1
    for j in xrange(min(k, N-k)):
        val = (val*(N-j))//(j+1)
    return val

而你建议的更新程序是：

    ntok = 1
    ktok = 1
    for t in xrange(1, min(k, n - k) + 1):
        ntok *= n
        ktok *= t
        n -= 1
    return ntok // ktok

我猜测SciPy实现较慢的原因是子程序在每次迭代中涉及整数除法，而您的实现仅在返回语句中调用一次除法。

- Victor