计算有和没有 SciPy 的 k 组合的数量

Question

comb of SciPy 函数似乎比简单的 Python 实现要慢，这让我感到困惑。这是解决Problem 53 of Project Euler 的两个等效程序的测量时间。

---

使用 SciPy：

%%timeit
from scipy.misc import comb

result = 0
for n in range(1, 101):
    for k in range(1, n + 1):
        if comb(n, k) > 1000000:
            result += 1
result

输出：

1 loops, best of 3: 483 ms per loop

---

没有 SciPy：

%%timeit
from math import factorial

def comb(n, k):
    return factorial(n) / factorial(k) / factorial(n - k)

result = 0
for n in range(1, 101):
    for k in range(1, n + 1):
        if comb(n, k) > 1000000:
            result += 1
result

输出：

10 loops, best of 3: 86.9 ms per loop

---

第二个版本大约快 5 倍（在两台 Mac 上测试，python-2.7.9-1，IPython 2.3.1-py27_0）。这是 SciPy (source code) 的 comb 函数的预期行为吗？我做错了什么？

---

编辑（来自 Anaconda 3.7.3-py27_0 发行版的 SciPy）：

import scipy; print scipy.version.version
0.12.0

---

编辑（IPython 之外的相同差异）：

$ time python with_scipy.py
real    0m0.700s
user    0m0.610s
sys     0m0.069s

$ time python without_scipy.py
real    0m0.134s
user    0m0.099s
sys     0m0.015s

Answer 1

看起来您可能错误地运行了计时并测量了将scipy 加载到内存中所需的时间。当我跑步时：

import timeit
from scipy.misc import comb
from math import factorial

def comb2(n, k):
    return factorial(n) / factorial(k) / factorial(n - k)

def test():
    result = 0
    for n in range(1, 101):
        for k in range(1, n + 1):
            if comb(n, k) > 1000000:
                result += 1

def test2():
    result = 0
    for n in range(1, 101):
        for k in range(1, n + 1):
            if comb2(n, k) > 1000000:
                result += 1

T = timeit.Timer(test)
print T.repeat(3,50)

T2 = timeit.Timer(test2)
print T2.repeat(3,50)

我明白了：

[2.2370951175689697, 2.2209839820861816, 2.2142510414123535]
[2.136591911315918, 2.138144016265869, 2.1437559127807617]

这表明 scipy 比 python 版本稍微快。

Answer 2

回答我自己的问题。 It seems 在 SciPy 中有两个不同的函数用于同一事物。我不太清楚为什么。但是另一个binom比comb快8.5倍，比我的快1.5倍，这让人放心：

%%timeit
from scipy.special import binom
result = 0
for n in range(1, 101):
    for k in range(1, n + 1):
        if binom(n, k) > 1000000:
            result += 1
result

输出：

10 loops, best of 3: 56.3 ms per loop

SciPy 0.14.0 伙计们，这对你也有用吗？

Answer 3

我相信这比提到的其他纯 python 方法更快：

from math import factorial
from operator import mul
def n_choose_k(n, k):
    if k < n-k:
        return reduce(mul, xrange(n-k+1, n+1)) // factorial(k)
    else:
        return reduce(mul, xrange(k+1, n+1)) // factorial(n-k)

与 numpy 解决方案相比，它的速度很慢，但请注意，NumPy 不能像 Python 那样在“无限整数”中工作。这意味着，虽然速度很慢，但 Python 将设法返回正确的结果。 NumPy 的默认行为对于组合数学并不总是理想的（至少在涉及非常大的整数时不是这样）。

例如

In [1]: from scipy.misc import comb

In [2]: from scipy.special import binom

In [3]: %paste
    from math import factorial
    from operator import mul
    def n_choose_k(n, k):
        if k < n-k:
            return reduce(mul, xrange(n-k+1, n+1)) // factorial(k)
        else:
            return reduce(mul, xrange(k+1, n+1)) // factorial(n-k)

## -- End pasted text --

In [4]: n_choose_k(10, 3), binom(10, 3), comb(10, 3)
Out[4]: (120, 120.0, 120.0)

In [5]: n_choose_k(1000, 250) == factorial(1000)//factorial(250)//factorial(750)
Out[5]: True

In [6]: abs(comb(1000, 250) - n_choose_k(1000, 250))
Out[6]: 3.885085558125553e+230

In [7]: abs(binom(1000, 250) - n_choose_k(1000, 250))
Out[7]: 3.885085558125553e+230

错误如此之大并不奇怪，这只是截断的影响。为什么要截断？ NumPy 将自身限制为使用 64 位来表示整数。然而，对这么大的整数的要求要多得多：

In [8]: from math import log

In [9]: log((n_choose_k(1000, 250)), 2)  
Out[9]: 806.1764820287578

In [10]: type(binom(1000, 250))
Out[10]: numpy.float64