埃拉托斯特尼筛法速度比较：Python与Julia

Question

我写了一个小的埃拉托斯特尼筛法函数，分别用Python和Julia这两种语言实现，现在我在比较它们的运行时间。

这是Python的代码：

import time

def get_primes(n):
    numbers = set(range(n, 1, -1))
    primes = []
    while numbers:
        p = numbers.pop()
        primes.append(p)
        numbers.difference_update(set(range(p*2, n+1, p)))
    return primes

start = time.time()
get_primes(10000)
print time.time() - start

这是Julia的代码：

function get_primes(n)
        numbers = [2:n]
        primes = Int[]
        while numbers != []
                p = numbers[1]
                push!(primes, p)
                numbers = setdiff(numbers, [p*i for i=1:int(n/p)])
        end
        return primes
end

@time get_primes(10000);

Python的代码大约运行了0.005秒，而Julia的代码则花了大约0.5秒，这意味着Python的运行速度大约快了100倍。可能有很合理的原因解释这个现象，但我其实并不太明白自己在做什么。

Answer 1

我测试了很多方法，发现这个方法在8次不同的测试中是最快的。

# Julia 0.4.0 [Execution time = 95us (After warm up!!)]
function get_primes(n::Int64)
  numbers = fill(true,n)
  numbers[1] = false

  for i = 2:isqrt(n)
    if numbers[i]
      for j = i:div(n,i)
        numbers[i*j] = false
      end
    end
  end

  primes = find(numbers) # t=3-5us
  return primes
end

@time primes = get_primes(10000)
println(get_primes(100))

这个页面上的第一段Julia代码在我的电脑上计算n=10000大约用了1'000'000微秒，而这个方法只用了95微秒，也就是说快了10'000倍。

# Python 3.4 [Execution time = 5ms (Every time)]
def get_primes(n):
  m = n+1
  numbers = [True for i in range(m)]
  for i in range(2, int(math.sqrt(n))):
    if numbers[i]:
      for j in range(i*i, m, i):
        numbers[j] = False
  primes = []
  for i in range(2, m):
    if numbers[i]:
      primes.append(i)
  return primes

start = time.time()
primes = get_primes(10000)
print(time.time() - start)
print(get_primes(100))

Python测试

# Python n=100e6 [Execution time 52.906s]
start = time.time()
get_primes(int(100e6))
print(time.time() - start)

Julia测试

# Julia n=100e6 [Execution time 3.694s]
@time get_primes(convert(Int64,100e6))

现在两者之间的差距小了很多。Julia比Python快大约12倍。

Answer 2

尝试使用 setdiff! 而不是 setdiff

用这段代码

function get_primes(n)
   numbers::Set{Int} = Set(2:n)
   primes::Array{Int64,1} = []
   while !isempty(numbers) 
      p = minimum(numbers)
      push!(primes,p);
      setdiff!(numbers,Set(p:p:n))
   end
   return primes
end

我得到了

julia> @time get_primes(10000);
elapsed time: 0.029556727 seconds (1656448 bytes allocated)

另一个（原始的）版本确实不太好，因为它在计算后大部分时间都花在了重新哈希 setdiff 上——我测试的未修改版本的时间和提问者的差不多。所以看起来 setdiff! 可能比 setdiff 快 100 倍，但它只接受集合而不是数组。

这仍然比 Python 慢 6 倍，但比使用 setdiff 快 13 倍。不过，如果有办法保持一个有序的集合，并且总是取第一个元素，那么它可能会更快，因为几乎 90% 的时间（209/235）都花在了寻找集合中的最小值上。

235 task.jl; anonymous; line: 96
 235 REPL.jl; eval_user_input; line: 54
  235 profile.jl; anonymous; line: 14
   209 /Users/jeffw/src/errata/julia/sive.jl; get_primes; line: 5
    2   reduce.jl; mapfoldl; line: 75
     2 dict.jl; skip_deleted; line: 669
    207 reduce.jl; mapfoldl; line: 81
     1   reduce.jl; mapfoldl_impl; line: 54
      1 dict.jl; skip_deleted; line: 670
     199 reduce.jl; mapfoldl_impl; line: 57
      10  dict.jl; skip_deleted; line: 668
      132 dict.jl; skip_deleted; line: 669
      12  dict.jl; skip_deleted; line: 670
      27  dict.jl; skip_deleted; line: 672
     7   reduce.jl; mapfoldl_impl; line: 58
   1   /Users/jeffw/src/errata/julia/sive.jl; get_primes; line: 6
   25  /Users/jeffw/src/errata/julia/sive.jl; get_primes; line: 7
    14 set.jl; setdiff!; line: 24
     1 dict.jl; skip_deleted; line: 669

更新：改用 IntSet 和 shift!

function get_primes(n)
   numbers::IntSet = IntSet(2:n)
   primes::Array{Int64,1} = []
   while !isempty(numbers)
      p = shift!(numbers)
      push!(primes,p);
      setdiff!(numbers,Set(p:p:n))
   end
   return primes 
end

julia> @time get_primes(10000);
elapsed time: 0.003691856 seconds (1463152 bytes allocated)

Answer 3

主要的区别在于，在Python中你只分配了一组数据，叫做number，然后直接在这组数据上进行修改。而在Julia中，每次循环都会新建一个数组。还有一个区别是，Python使用的是集合，而Julia使用的是数组（在Python中称为“列表”）。我们来改一下Julia的代码，消除这两个区别：

function get_primes(n)
    numbers = IntSet(2:n)
    primes = Int[]
    while !isempty(numbers)
        p = shift!(numbers)
        push!(primes, p)
        setdiff!(numbers, p:p:n)
    end
    return primes
end

做了这个改动后，在我的系统上，Julia的代码比Python的快10倍（0.0005秒对0.0048秒）。值得注意的是，我在setdiff!的第二个参数中使用了一个范围，这样写更简洁，也比用列表推导式构建数组快1.5倍。

在Julia中，更符合习惯的写法是使用一个布尔数组，来控制开关：

function eratosthenes(n)
    primes = fill(true,n)
    primes[1] = false
    for p = 2:n
        primes[p] || continue
        for i = 2:div(n,p)
            primes[p*i] = false
        end
    end
    find(primes)
end

最后的find会返回一个向量中非零（也就是为真的）元素的索引。在我的机器上，这种写法比之前的Julia版本快5倍（0.0001秒），也比Python版本快45倍。