Legendre prime counting function: Difference between revisions

Content added Content deleted
m (use the precalculated values from Nim example)
Line 2,063: Line 2,063:


function countprimes(n)
function countprimes(n)
n < 3 && return n > 1
n < 3 && return typeof(n)(n > 1)
rtlmt = isqrt(n)
rtlmt = isqrt(n)
mxndx = (rtlmt - 1) ÷ 2
mxndx = (rtlmt - 1) ÷ 2
smalls::Array{UInt32} = [i for i in 0:mxndx+1]
rilmt = mxndx + 1
smalls = [i for i in 0:rilmt]
roughs::Array{UInt32} = [i + i + 1 for i in 0:mxndx+1]
roughs = [i + i + 1 for i in 0:rilmt]
larges::Array{Int64} = [(n ÷ (i + i + 1) - 1) ÷ 2 for i in 0:mxndx+1]
larges = [(n ÷ (i + i + 1) - 1) ÷ 2 for i in 0:rilmt]
cmpsts = falses(mxndx + 1)
cullbuf = zeros(Int, (mxndx + 8) ÷ 8)
bp, npc, mxri = 3, 0, mxndx
nbps = 0
while true
for i in 1:typemax(UInt32)
i = bp >> 1
sqri = (i + i) * (i + 1)
sqri = (i + i) * (i + 1)
sqri > mxndx && break
sqri > mxndx && break
cullbuf[i >> 3 + 1] & masks[i & 7 + 1] != 0 && continue
if !cmpsts[i + 1]
cullbuf[i >> 3 + 1] |= masks[i & 7 + 1]
cmpsts[i + 1] = true
bp = i + i + 1
for c in sqri:bp:mxndx
for c in sqri:bp:arrlen-1
cmpsts[c + 1] = true
⚫
cullbuf[c >> 3 + 1] |= masks[c & 7 + 1]
⚫
end
⚫
nri = 0
⚫
for ori in 0:rilmt-1
⚫
r = roughs[ori + 1]
⚫
rci = r >> 1
cullbuf[rci >> 3 + 1] & masks[rci & 7 + 1] != 0 && continue
⚫
d = r * bp
⚫
t = 0
⚫
if d <= rtlmt
⚫
t = larges[smalls[d>>1 + 1] - nbps + 1]
⚫
else
⚫
t = smalls[(n ÷ d - 1) >> 1 + 1]
end
end
larges[nri + 1] = larges[ori + 1] - t + nbps
ri = 0
roughs[nri + 1] = r
for k in 0:mxri
nri += 1
q = roughs[k + 1]
end
qi = q >> 1
si = mxndx
cmpsts[qi + 1] && continue
for pm in (rtlmt ÷ bp - 1) | 1 : -2 : bp
d = UInt64(bp) * UInt64(q)
c = smalls[pm>>1 + 1]
larges[ri + 1] = larges[k + 1] + npc -
e = (pm * bp) >> 1
(d <= rtlmt ? larges[smalls[d >> 1 + 1] - npc + 1]
while si >= e
: smalls[(Int(floor(n / d)) - 1) >> 1 + 1])
smalls[si + 1] -= (c - nbps)
roughs[ri + 1] = q
si -= 1
ri += 1
end
end
⚫
m = mxndx
for k in (rtlmt ÷ bp - 1) | 1 : -2 : bp
⚫
c = smalls[k >> 1 + 1] - npc
⚫
ee = (k * bp) >> 1
⚫
while m >= ee
⚫
smalls[m + 1] -= c
⚫
m -= 1
⚫
end
⚫
end
⚫
mxri = ri - 1
⚫
npc += 1
end
end
rilmt = nri
bp += 2
⚫
nbps += 1
end
end

ans = larges[1] + ((rilmt + 2*(nbps-1)) * (rilmt - 1) ÷ 2)
for ri in 1:rilmt-1
result = larges[1]
ans -= larges[ri + 1]
for i in 2:mxri+1
⚫
result -= larges[i]
end
end
result += (mxri + 1 + 2 * (npc - 1)) * mxri ÷ 2
for ri in 1:typemax(UInt32)

p = roughs[ri + 1]
⚫
for j in 1:mxri
⚫
p = UInt64(roughs[j + 1])
m = n ÷ p
m = n ÷ p
ei = smalls[((m ÷ p) - 1) >> 1 + 1] - nbps
ee = smalls[(Int(floor(m / p)) - 1) >> 1 + 1] - npc
ei <= ri && break
ee <= j && break
ans -= (ei - ri) * (nbps + ri - 1)
for k in j+1:ee
⚫
result += smalls[(Int(floor(m / roughs[k + 1])) - 1) >> 1 + 1]
for sri in ri+1:ei
⚫
ans += smalls[(m ÷ (roughs[sri + 1]) - 1) >> 1 + 1]
end
end
result -= (ee - j) * (npc + j - 1)
end
end
return ans + 1
return result + 1
end
end


Line 2,129: Line 2,129:
println("π(10^$i) = ", countprimes(10^i))
println("π(10^$i) = ", countprimes(10^i))
end
end

@time countprimes(10^14)

</syntaxhighlight>{{out}}
</syntaxhighlight>{{out}}
<pre>
<pre>
Line 2,146: Line 2,149:
π(10^13) = 346065536839
π(10^13) = 346065536839
π(10^14) = 3204941750802
π(10^14) = 3204941750802
4.744442 seconds (13 allocations: 153.185 MiB, 0.16% gc time)
</pre>
</pre>


Retrieved from "https://rosettacode.org/wiki/Legendre_prime_counting_function"