74 lines
3.3 KiB
Forth
74 lines
3.3 KiB
Forth
let masks = Array.init 8 ((<<<) 1uy) // quick bit twiddling
|
|
|
|
let countPrimes (lmt: uint64): int64 =
|
|
if lmt < 3UL then (if lmt < 2UL then 0L else 1L) else
|
|
let inline half x = (x - 1) >>> 1
|
|
let inline divide nm d = (float nm) / (float d) |> int
|
|
let sqrtlmt = lmt |> float |> sqrt |> uint64
|
|
let mxndx = (sqrtlmt - 1UL) / 2UL |> int in let cbsz = (mxndx + 8) / 8
|
|
let cullbuf = Array.zeroCreate cbsz
|
|
let smalls = Array.init (mxndx + 1) uint32
|
|
let roughs = Array.init (mxndx + 1) <| fun i -> uint32 (i + i + 1)
|
|
let larges = Array.init (mxndx + 1) <| fun i ->
|
|
int64 (lmt / uint64 (i + i + 1) - 1UL) / 2L
|
|
|
|
let rec loopbp bp nobps rilmt =
|
|
let i = int (bp - 1UL) >>> 1 in let sqri = (i + i) * (i + 1)
|
|
if sqri > mxndx then nobps, rilmt else
|
|
if (cullbuf.[i >>> 3] &&& masks.[i &&& 7]) <> 0uy then
|
|
loopbp (bp + 2UL) nobps rilmt else
|
|
let w = i >>> 3 in cullbuf.[w] <- cullbuf.[w] ||| masks.[i &&& 7] // cull bp
|
|
{ sqri .. int bp .. mxndx } |> Seq.iter (fun c -> // cull multiples of bp...
|
|
let w = c >>> 3 in cullbuf.[w] <- cullbuf.[w] ||| masks.[c &&& 7] )
|
|
|
|
// adjust `larges for last partial loop pass;
|
|
// compress larges/roughs for current partial sieve pass...
|
|
let rec loopri iri ori =
|
|
if iri > rilmt then ori - 1 else
|
|
let r = uint64 roughs.[iri] in let sri = int (r >>> 1)
|
|
if (cullbuf.[sri >>> 3] &&& masks.[sri &&& 7]) <> 0uy then
|
|
loopri (iri + 1) ori else // skip for roughs culled this pass!
|
|
let d = bp * r
|
|
larges.[ori] <- larges.[iri] -
|
|
( if d <= sqrtlmt then
|
|
larges.[int smalls.[int (d >>> 1)] - nobps]
|
|
else let ndx = (half << divide lmt) d
|
|
int64 smalls.[ndx] ) + int64 nobps
|
|
roughs.[ori] <- uint32 r; loopri (iri + 1) (ori + 1)
|
|
|
|
// adjust `smalls` for last partial loop pass...
|
|
let rec loopbpm bpm mxsi =
|
|
if bpm < bp then () else
|
|
let c = smalls.[int (bpm >>> 1)] - uint32 nobps
|
|
let ei = (bpm * bp) >>> 1 |> int
|
|
let rec loopsi si =
|
|
if si < ei then si else smalls.[si] <- smalls.[si] - c; loopsi (si - 1)
|
|
loopbpm (bpm - 2UL) (loopsi mxsi)
|
|
|
|
let nrilmt = loopri 0 0
|
|
loopbpm ((sqrtlmt / bp - 1UL) ||| 1UL) mxndx
|
|
loopbp (bp + 2UL) (nobps + 1) nrilmt
|
|
|
|
// accumulate result so far; compensate for over subtraction...
|
|
let numobps, mxri = loopbp 3UL 0 mxndx
|
|
let rec smr i acc =
|
|
if i > mxri then // adjust accumulated answer!
|
|
acc + (int64 mxri + 1L + 2L * int64 (numobps - 1)) * int64 mxri / 2L
|
|
else smr (i + 1) (acc - larges.[i])
|
|
let ans0 = smr 1 larges.[0]
|
|
|
|
// finally, add result from pairs of rough primes up to cube root of range,
|
|
// where they are two different primes; compensating for over addition...
|
|
let rec loopri ri acc =
|
|
let p = uint64 roughs.[ri] in let q = lmt / p
|
|
let ei = int smalls.[(half << divide q) p] - numobps
|
|
if ei <= ri then acc else
|
|
let rec loopori ori oacc =
|
|
if ori > ei then oacc else
|
|
let ndx = (half << divide q) (uint64 roughs.[ori])
|
|
loopori (ori + 1) (oacc + int64 smalls.[ndx])
|
|
let nacc = loopori (ri + 1) acc // subtract over addition of base primes:
|
|
loopri (ri + 1) (nacc - int64 (ei - ri) * (int64 numobps + int64 ri - 1L))
|
|
|
|
loopri 1 ans0 + 1L // add one for only even prime of two!
|