PrimeSieves.jl

#==============================================================================#
# Optimized Sieve of Atkin
#==============================================================================#

function atkin(n::Int)
    s = zeros(Bool, n < 2 ? 0 : n)
    n < 2 && return s; s[2] = true
    n < 3 && return s; s[3] = true
    for x = 1:floor(Int,sqrt(n/4))
        j = 4 * x * x
        y, Δy = 1, 2
        if x % 3 == 0
            while (k = j + y * y) ≤ n # k = 4x² + y²
                @inbounds s[k] = !s[k]
                Δy $= 6; y += Δy
            end
        else
            while (k = j + y * y) ≤ n
                @inbounds s[k] = !s[k]
                y += 2
            end
        end
    end
    for x=1:2:floor(Int,sqrt(n/3))
        j = 3 * x * x
        y, Δy = 2, 4
        while (k = j + y * y) ≤ n # k = 3x² + y²
            @inbounds s[k] = !s[k]
            Δy $= 6; y += Δy
        end
    end
    for x = 1:floor(Int,sqrt(n/2))
        j = 3 * x * x
        y, Δy = ifelse(x % 2 == 0, (1, 2), (2, 4))
        while x > y
            (k = j - y * y) ≤ n && (@inbounds s[k] = !s[k]) # k = 3x² - y²
            Δy $= 6; y += Δy
        end
    end
    @inbounds for p = 5:isqrt(n)
        if s[p]; for j = p*p:2p*p:n; s[j] = false; end; end
    end
    return s
end

# Smallest multiple of k that is greater than or equal to both k and l
@inline smallestmultiple(k::Int, l::Int) = ifelse(k ≥ l, k, k * ceil(Int, l / k))

function atkin(lo::Int, hi::Int)
    lo ≤ hi || throw(ArgumentError("the condition lo ≤ hi must be met"))
    lo < 2 && return atkin(hi)
    n = hi - lo + 1
    s = zeros(Bool, n)
    smallprimes = primes(isqrt(hi))
    lo ≤ 2 ≤ hi && (s[3-lo] = true)
    lo ≤ 3 ≤ hi && (s[4-lo] = true)
    for x = 1:floor(Int,sqrt(hi/4))
        j = 4 * x * x
        y, Δy = 1, 2
        if j < lo
            y = 2 * ceil(Int, (sqrt(lo - j) - 1) / 2) + 1
        end
        if x % 3 == 0
            if y % 3 != 1; Δy = 4; if y % 3 == 0; y += 2; end; end
            while (k = j + y * y - lo + 1) ≤ n # k + lo - 1 = 4x² + y²
                @inbounds s[k] = !s[k]
                Δy $= 6; y += Δy
            end
        else
            while (k = j + y * y - lo + 1) ≤ n
                @inbounds s[k] = !s[k]
                y += 2
            end
        end
    end
    for x = 1:2:floor(Int,sqrt(hi/3))
        j = 3 * x * x
        y, Δy = 2, 4
        if j < lo
            y = 2 * ceil(Int, sqrt(lo - j) / 2)
            if y % 3 == 0; y += 2; elseif y % 3 == 1; Δy = 2; end
        end
        while (k = j + y * y - lo + 1) ≤ n # k = 3x² + y²
            @inbounds s[k] = !s[k]
            Δy $= 6; y += Δy
        end
    end
    for x = 1:floor(Int,sqrt(hi/2))
        j = 3 * x * x
        y, Δy = ifelse(x % 2 == 0, (1, 2), (2, 4))
        while x > y
            0 < (k = j - y * y - lo + 1) ≤ n && (@inbounds s[k] = !s[k])
            Δy $= 6; y += Δy
        end
    end
    @inbounds for i = 3:length(smallprimes)
        p = smallprimes[i]
        for j = smallestmultiple(p*p,lo)-lo+1:p*p:n; s[j] = false; end
    end
    return s
end