Atkin の篩と，チューンナップした Eratosthenes の篩 速度はほぼ同等！！

タイトルどおり。

エラトステネスの篩も，さらなるチューニングはあるが，まあ程々で。
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