⏱️

N番目の素数を求める (Julia版)

7 min read 2

すぎゃーんさんのメモを読んだので,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**を呼び出します.

結果まとめ

以下が結果のまとめ (n=1000,10000,100000).

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

ベンチマークの結果(n=1000,10000)です.

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

ベンチマークの結果(n=1000,10000)です.

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

ベンチマークの結果(n=1000,10000)です.

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

ベンチマークの結果(n=1000,10000)です.

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

終わりです.