Eratosthenes 筛网速度比较：Python vs Julia答案

【问题标题】：Sieve of Eratosthenes Speed Comparison: Python vs JuliaEratosthenes 筛网速度比较：Python vs Julia
【发布时间】：2015-08-14 22:32:42
【问题描述】：

所以我有一个用 Python 和 Julia 编写的 Eratosthenes 函数的小筛子，我正在比较运行时间。

这是 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 倍。这可能有一个完全合乎逻辑的原因，但我真的不知道我在做什么。

【问题讨论】：

为什么你在 python 中使用 set 而在 Julia 中使用数组？
另外，尝试对第二次运行计时 - 第一次有编译开销
@IainDunning，Python 代码是从其他人那里逐字提取的，我在 Julia 中重写并清理了一下，但没有做任何清理 Python 代码的操作
我认为 Python 代码有问题。 set 没有排序，尽管由于实现细节它可能似乎被排序。筛子非常需要按顺序计算素数。
@user1378571 不符合the documentation：“从集合中移除并返回任意元素。”

标签： python performance julia

【解决方案1】：

主要区别在于，在 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 倍。

【讨论】：

我在 Set 上使用 getindex 时遇到问题。这是最近的变化吗？
@waTeim：我必须将 Set(1:n) 更改为 Set([1:n]...) 并将 p = numbers[1] 更改为 p = minimum(numbers) 才能运行代码，对我而言，这使得代码现在只比 python 慢 10 倍。第一个更改可能是（非常）最近的构建问题，但似乎 get_index 不应该是集合数据结构所具有的方法。
嗯。剪切和粘贴代码时我一定做错了什么，因为现在这对我不起作用。正在进行更新。
我想我之前不小心定时了内置的primes函数，之前版本的代码从来没有工作过。相反，我现在使用的是IntSet 和shift!，它们确实有效并且相当快。惯用的布尔数组版本现在比 Julia 和 Python 集合版本都好很多。
@user1378571 我也是这么想的，但是貌似IntSet是有序的

【解决方案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 重新散列上——我对未更改版本的时间安排与 OP 相似。所以看起来 setdiff! 可能比 setdiff 快 100 倍，但只接受 Sets 而不是数组.

这仍然比 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)

【讨论】：

【解决方案3】：

我测试了很多方法，我发现这个是 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 代码在我的 PC 上以大约 1'000'000us 计算 n=10000，而这个大约需要 95us，这比我的 PC 快 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 n=100e6 [Execution time 3.694s]
@time get_primes(convert(Int64,100e6))

现在差异变小了。 Julia 比 Python 快 12 倍左右。

【讨论】：