タイトルどおり。
エラトステネスの篩も，さらなるチューニングはあるが，まあ程々で。
C で書いたプログラムをコールすれば結果は変わる。
Julia の primes は atkin の篩によるものらしいが，速い（コンパイラ言語のプログラムを読んでいるのだろう）。
Julia の Primes.primes による場合
using Primes
@time p = primes(10^9); # 2.205441 seconds
アトキンの篩
#=
The MIT License (MIT)
Copyright (c) 2022/01/09 r-de-r
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
=#
function Atkin(limit=10^6)
""" n 以下の素数をアトキンの篩で生成 """
sqrt_limit = isqrt(limit)
isprime = falses(limit)
for z in [1, 5]
for y in z:6:sqrt_limit
y² = y*y
max_x = floor(Int, sqrt((limit - y²) / 4))
for i in 1:max_x
n = 4i*i + y²
isprime[n] = !isprime[n]
end
x = y + 1
max_x = floor(Int, sqrt((limit + y²) / 3))
if x <= max_x
for i in x:2:max_x
n = 3i*i - y²
isprime[n] = !isprime[n]
end
end
end
end
for z in [2, 4]
for y in z:6:sqrt_limit
y² = y*y
max_x = floor(Int, sqrt((limit - y²) / 3))
for i in 1:2:max_x
n = 3i*i + y²
isprime[n] = !isprime[n]
end
x = y + 1
max_x = floor(Int, sqrt((limit + y²) / 3))
if x <= max_x
for i in x:2:max_x
n = 3i*i - y²
isprime[n] = !isprime[n]
end
end
end
end
for y in 3:6:sqrt_limit
y² = y*y
for x in 1:2
max_x = floor(Int, sqrt((limit - y²) / 4))
for i in x:3:max_x
n = 4i*i + y²
isprime[n] = !isprime[n]
end
end
end
for n in 5:2:sqrt_limit
if isprime[n]
for i in n*n:n*n:limit
isprime[i] = false
end
end
end
isprime[2] = isprime[3] = true
(1:limit)[isprime]
end
@time a2 = Atkin(10^9); # 3.080654 seconds
エラトステネスの篩
#=
The MIT License (MIT)
Copyright (c) 2022/01/09 r-de-r
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
=#
# https://qiita.com/peria/items/54499b9ce9d5c1e93e5a を参考に
""" 偶数は対象としないバージョン """
function Eratosthenes(n)
tbl = trues(ceil(Int, n / 2))
tbl[1] = false
sqrt_x = isqrt(n)+1
sqrt_xi = floor(Int, sqrt_x / 2)
for i in 2:sqrt_xi
if tbl[i]
for mult in 2i * (i - 1) + 1:2i - 1:length(tbl)
tbl[mult] = false
end
end
end
pushfirst!(((1:length(tbl))[tbl]) .* 2 .- 1, 2)
end
@time S = Eratosthenes(10^9); # 3.084727 seconds
