45 lines
1.7 KiB
Nim
45 lines
1.7 KiB
Nim
# compile with: nim c -d:danger -t:-march=native --gc:arc
|
|
|
|
from std/monotimes import getMonoTime, `-`
|
|
from std/times import inMilliseconds
|
|
from std/math import sqrt
|
|
|
|
let masks = [ 1'u8, 2, 4, 8, 16, 32, 64, 128 ] # faster than bit twiddling
|
|
let masksp = cast[ptr[UncheckedArray[byte]]](unsafeAddr(masks[0]))
|
|
|
|
proc countPrimes(n: int64): int64 =
|
|
if n < 3:
|
|
return if n < 2: 0 else: 1
|
|
else:
|
|
let rtlmt = n.float64.sqrt.int
|
|
let mxndx = (rtlmt - 1) div 2
|
|
let sz = (mxndx + 8) div 8
|
|
var cmpsts = cast[ptr[UncheckedArray[byte]]](alloc0(sz))
|
|
for i in 1 .. mxndx:
|
|
if (cmpsts[i shr 3] and masksp[i and 7]) != 0: continue
|
|
let sqri = (i + i) * (i + 1)
|
|
if sqri > mxndx: break
|
|
let bp = i + i + 1
|
|
for c in countup(sqri, mxndx, bp):
|
|
let w = c shr 3; cmpsts[w] = cmpsts[w] or masksp[c and 7]
|
|
var pisqrt = 0'i64
|
|
for i in 0 .. mxndx:
|
|
if (cmpsts[i shr 3] and masksp[i and 7]) == 0: pisqrt += 1
|
|
var primes = cast[ptr[UncheckedArray[uint32]]](alloc(sizeof(uint32) * pisqrt.int))
|
|
var j = 0
|
|
for i in 0 .. mxndx:
|
|
if (cmpsts[i shr 3] and masksp[i and 7]) == 0: primes[j] = (i + i + 1).uint32; j += 1
|
|
proc phi(x: int64; a: int): int64 =
|
|
if a <= 1:
|
|
return if a < 1: x else: x - (x shr 1)
|
|
let p = primes[a - 1].int64
|
|
if x <= p: return 1 # very simple one-line optimization that limits exponential growth!
|
|
return phi(x, a - 1) - phi((x.float64 / p.float64).int64, a - 1)
|
|
result = phi(n, pisqrt.int) + pisqrt - 1
|
|
cmpsts.dealloc; primes.dealloc
|
|
|
|
let nstrt = getMonoTime()
|
|
var pow = 1'i64
|
|
for i in 0 .. 9: echo "π(10^", i, ") = ", pow.countPrimes; pow *= 10
|
|
let nelpsd = (getMonoTime() - nstrt).inMilliseconds
|
|
echo "This took ", nelpsd, " milliseconds."
|