66 lines
2.7 KiB
Nim
66 lines
2.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]))
|
|
|
|
const TinyPhiPrimes = [2, 3, 5, 7, 11, 13]
|
|
const TinyPhiCirc = 3 * 5 * 7 * 11 * 13
|
|
const TinyPhiRes = 2 * 4 * 6 * 10 * 12
|
|
const CC = 6
|
|
proc makeTinyPhiLUT(): array[TinyPhiCirc, uint16] =
|
|
for i in 0 .. TinyPhiCirc - 1: result[i] = 1
|
|
for i in 1 .. 6:
|
|
if result[i] == 0: continue
|
|
result[i] = 0; let bp = i + i + 1
|
|
let sqri = (i + i) * (i + 1)
|
|
for c in countup(sqri, TinyPhiCirc - 1, bp): result[c] = 0
|
|
var acc = 0'u16;
|
|
for i in 0 .. TinyPhiCirc - 1: acc += result[i]; result[i] = acc
|
|
const TinyPhiLUT = makeTinyPhiLUT()
|
|
proc tinyPhi(x: int64): int64 {.inline.} =
|
|
let ndx = (x - 1) div 2; let numtot = ndx div TinyPhiCirc.int64
|
|
return numtot * TinyPhiRes.int64 + TinyPhiLUT[(ndx - numtot * TinyPhiCirc.int64).int].int64
|
|
|
|
proc countPrimes(n: int64): int64 =
|
|
if n < 169: # below 169 whose sqrt is 13 is where TinyPhi doesn't work...
|
|
if n < 3: return if n < 2: 0 else: 1
|
|
# adjust for the missing "degree" base primes
|
|
if n <= 13: return (n - 1) div 2 + (if n < 9: 1 else: 0)
|
|
return 5 + TinyPhiLUT[(n - 1).int div 2].int64
|
|
let rtlmt = n.float64.sqrt.int
|
|
let mxndx = (rtlmt - 1) div 2
|
|
var cmpsts = cast[ptr[UncheckedArray[byte]]](alloc0((mxndx + 8) div 8))
|
|
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
|
|
let 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
|
|
var phi = tinyPhi(n)
|
|
proc lvl(m, mbf: int64; mxa: int) = # recurse from bottom left of "tree"...
|
|
for a in CC .. mxa:
|
|
let p = primes[a].int64
|
|
if m < p * p: phi += mbf * (mxa - a + 1).int64; return # rest of level all ones!
|
|
let nm = (m.float64 / p.float64).int64; phi += mbf * tinyPhi(nm)
|
|
if a > CC: lvl(nm, -mbf, a - 1) # split
|
|
# finished level!
|
|
lvl(n, -1, pisqrt.int - 1); result = phi + pisqrt - 1
|
|
cmpsts.dealloc; primes.dealloc
|
|
|
|
let strt = getMonoTime()
|
|
var pow = 1'i64
|
|
for i in 0 .. 9: echo "π(10^", i, ") = ", pow.countPrimes; pow *= 10
|
|
let elpsd = (getMonoTime() - strt).inMilliseconds
|
|
echo "This took ", elpsd, " milliseconds."
|