⏱️
N番目の素数を求める (Julia版)
すぎゃーんさんのメモを読んだので，Juliaで書いてみます (Juliaの書き方が最適化されているとはいえないので，注意です)．アトキンの篩は読んでもよく分からなかったので，Python版の実装を真似して書いてみます．
Juliaの環境はこちらです．計測はBenchmarkTools.jlの**@benchmark**でやります．
julia> versioninfo()
Julia Version 1.5.3
Commit 788b2c77c1* (2020-11-09 13:37 UTC)
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: Intel(R) Core(TM) i7-8550U CPU @ 1.80GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-9.0.1 (ORCJIT, skylake)
実行の仕方としては関数をまとめて書いたスクリプトファイルをRevise.jl経由でincludetし，動作を確認したのにち**@benchmark**を呼び出します．
結果まとめ
以下が結果のまとめ (
main(1000) : 7.540 ms (0 allocations: 0 bytes)
main2(1000) : 1.355 ms (10 allocations: 16.39 KiB)
main3(1000) : 56.400 μs (41 allocations: 60.63 KiB)
main3_max(1000) : 32.500 μs (11 allocations: 24.64 KiB)
main(10000) : 1.026 s (0 allocations: 0 bytes)
main2(10000) : 130.002 ms (14 allocations: 256.64 KiB)
main3(10000) : 897.900 μs (98 allocations: 777.80 KiB)
main3_max(10000) : 431.800 μs (16 allocations: 364.66 KiB)
main3(100000) : 15.866 ms (172 allocations: 10.93 MiB)
main3_max(100000) : 6.010 ms (19 allocations: 3.32 MiB)
内容
例1 順番に割っていく
どうやって書くかを少し悩んだ結果この形にしました．
function main(n)
i = 0
for num in 2:2_000_000
flag = true
for q in 2:num-1
if num % q == 0
flag = false
break
end
end
(flag) && (i += 1)
(i == n) && return num
end
return -1
end
ベンチマークの結果(
julia> @benchmark main(1000)
BenchmarkTools.Trial:
memory estimate: 0 bytes
allocs estimate: 0
--------------
minimum time: 7.539 ms (0.00% GC)
median time: 7.975 ms (0.00% GC)
mean time: 8.107 ms (0.00% GC)
maximum time: 11.950 ms (0.00% GC)
--------------
samples: 617
evals/sample: 1
julia> @benchmark main(10000)
BenchmarkTools.Trial:
memory estimate: 0 bytes
allocs estimate: 0
--------------
minimum time: 1.060 s (0.00% GC)
median time: 1.078 s (0.00% GC)
mean time: 1.084 s (0.00% GC)
maximum time: 1.132 s (0.00% GC)
--------------
samples: 5
evals/sample: 1
例2 試し割り
primesのVectorを外に出しつつ，関数を中身に入れ込んで書きました．
function main2(n)
primes = Int[]
function is_prime(num)
for p in primes
(num % p == 0) && return false
end
push!(primes, num)
return true
end
i = 0
for num in 2:2_000_000
is_prime(num) && (i += 1)
(i == n) && return num
end
return -1
end
ベンチマークの結果(
julia> @benchmark main2(1000)
BenchmarkTools.Trial:
memory estimate: 16.39 KiB
allocs estimate: 10
--------------
minimum time: 1.369 ms (0.00% GC)
median time: 1.394 ms (0.00% GC)
mean time: 1.439 ms (0.03% GC)
maximum time: 2.807 ms (46.59% GC)
--------------
samples: 3475
evals/sample: 1
julia> @benchmark main2(10000)
BenchmarkTools.Trial:
memory estimate: 256.64 KiB
allocs estimate: 14
--------------
minimum time: 129.550 ms (0.00% GC)
median time: 132.639 ms (0.00% GC)
mean time: 135.832 ms (0.00% GC)
maximum time: 169.846 ms (0.00% GC)
--------------
samples: 37
evals/sample: 1
例3 篩
Pythonのフラグ管理版を真似して書いたものです．
function main3(n)
function list_primes(limit)
primes = Int[]
is_prime = [true for _ in 1:limit]
is_prime[1] = false
for p in 1:limit
(!is_prime[p]) && continue
push!(primes, p)
for i in p*p:p:limit
is_prime[i] = false
end
end
return primes
end
limit = 1000
while true
primes = list_primes(limit)
if length(primes) > n
return primes[n]
end
limit *= 2
end
end
ベンチマークの結果(
julia> @benchmark main3(1000)
BenchmarkTools.Trial:
memory estimate: 60.63 KiB
allocs estimate: 41
--------------
minimum time: 66.000 μs (0.00% GC)
median time: 73.800 μs (0.00% GC)
mean time: 81.005 μs (2.54% GC)
maximum time: 2.721 ms (82.98% GC)
--------------
samples: 10000
evals/sample: 1
julia> @benchmark main3(10000)
BenchmarkTools.Trial:
memory estimate: 777.80 KiB
allocs estimate: 98
--------------
minimum time: 1.079 ms (0.00% GC)
median time: 1.152 ms (0.00% GC)
mean time: 1.213 ms (1.27% GC)
maximum time: 3.723 ms (0.00% GC)
--------------
samples: 4120
evals/sample: 1
例3+最大推定
こちらもPythonの実装を真似して書いたものです．
function main3_max(n)
limit = max(100, floor(Int, n * log(n) * 1.2))
primes = Int[]
is_prime = [true for _ in 1:limit]
is_prime[1] = false
for p in 1:limit
(!is_prime[p]) && continue
push!(primes, p)
(length(primes) == n) && return primes[end]
for i in p*p:p:limit
is_prime[i] = false
end
end
end
ベンチマークの結果(
julia> @benchmark main3_max(1000)
BenchmarkTools.Trial:
memory estimate: 24.64 KiB
allocs estimate: 11
--------------
minimum time: 28.300 μs (0.00% GC)
median time: 31.000 μs (0.00% GC)
mean time: 32.821 μs (0.78% GC)
maximum time: 1.516 ms (84.89% GC)
--------------
samples: 10000
evals/sample: 1
julia> @benchmark main3_max(10000)
BenchmarkTools.Trial:
memory estimate: 364.66 KiB
allocs estimate: 16
--------------
minimum time: 366.500 μs (0.00% GC)
median time: 385.950 μs (0.00% GC)
mean time: 413.729 μs (0.53% GC)
maximum time: 2.947 ms (0.00% GC)
--------------
samples: 10000
evals/sample: 1
例4. Primes.jlパッケージ
ついでに適当に見つけてきたPrimes.jlパッケージを使った場合を見てみます．使い方は**prime()**を呼び出すだけです．そんなに高速ではないみたいでしょうか．実装はこちらの行周辺ですが，次の素数を求める関数(nextprime)というやつを使っているみたいですね．実装を眺めてみると，ずらしていくような計算をしているっぽいです．
julia> @benchmark prime(Int,1000)
BenchmarkTools.Trial:
memory estimate: 18.05 KiB
allocs estimate: 1155
--------------
minimum time: 403.400 μs (0.00% GC)
median time: 453.000 μs (0.00% GC)
mean time: 503.771 μs (0.10% GC)
maximum time: 8.601 ms (0.00% GC)
--------------
samples: 9910
evals/sample: 1
julia> @benchmark prime(Int,10000)
BenchmarkTools.Trial:
memory estimate: 265.56 KiB
allocs estimate: 16996
--------------
minimum time: 7.189 ms (0.00% GC)
median time: 9.174 ms (0.00% GC)
mean time: 9.718 ms (0.09% GC)
maximum time: 24.044 ms (0.00% GC)
--------------
samples: 515
evals/sample: 1
終わりです．
Discussion
Primes.jlってそんなに遅かったっけ？と思ったらprime()の実装はあれは遅いですね、なんであんな実装になってるんでしょうね。ということで自分なりに
Primes.jlの機能を使ってオレオレ実装してみました。primesmask(lo, hi)は、lo以上hi以下の数が素数かどうかを表す要素数hi-lo+1のBitVector（<:AbstractVector{Bool}） を返す関数です。ベンチマーク：
手元で
main3_max(n)も実装してみましたが、この記事のベンチマーク結果とほぼ同じで、それと比較して↑の実装は所要時間もメモリ使用量も半分くらいになりました。割と満足。ご参考までに。
ありがとうございます！あれは完全に実装が謎だと思いました，こっちの版は良いですね！勉強になります．