119 lines
6.2 KiB
Haskell
119 lines
6.2 KiB
Haskell
countPrimes :: Word64 -> Int64
|
|
countPrimes n =
|
|
if n < 3 then (if n < 2 then 0 else 1) else
|
|
let
|
|
{-# INLINE divide #-}
|
|
divide :: Word64 -> Word64 -> Int
|
|
divide nm d = truncate $ (fromIntegral nm :: Double) / fromIntegral d
|
|
{-# INLINE half #-}
|
|
half :: Int -> Int
|
|
half x = (x - 1) `shiftR` 1
|
|
rtlmt = floor $ sqrt (fromIntegral n :: Double)
|
|
mxndx = (rtlmt - 1) `div` 2
|
|
(!nbps, !nrs, !smalls, !roughs, !larges) = runST $ do
|
|
mss <- unsafeNewArray_ (0, mxndx) :: ST s (STUArray s Int Word32)
|
|
let msscst =
|
|
castSTUArray :: STUArray s Int Word32 -> ST s (STUArray s Int Int64)
|
|
mdss <- msscst mss -- for use in adjing counts LUT
|
|
forM_ [ 0 .. mxndx ] $ \ i -> unsafeWrite mss i (fromIntegral i)
|
|
mrs <- unsafeNewArray_ (0, mxndx) :: ST s (STUArray s Int Word32)
|
|
forM_ [ 0 .. mxndx ] $ \ i -> unsafeWrite mrs i (fromIntegral i * 2 + 1)
|
|
mls <- unsafeNewArray_ (0, mxndx) :: ST s (STUArray s Int Int64)
|
|
forM_ [ 0 .. mxndx ] $ \ i ->
|
|
let d = fromIntegral (i + i + 1)
|
|
in unsafeWrite mls i (fromIntegral (divide n d - 1) `div` 2)
|
|
cmpsts <- unsafeNewArray_ (0, mxndx) :: ST s (STUArray s Int Bool)
|
|
|
|
let loop i !cbpi !rlmti =
|
|
let sqri = (i + i) * (i + 1) in
|
|
if sqri > mxndx then do
|
|
fss <- unsafeFreezeSTUArray mss
|
|
frs <- unsafeFreezeSTUArray mrs
|
|
fls <- unsafeFreezeSTUArray mls
|
|
return (cbpi, rlmti + 1, fss, frs, fls)
|
|
else do
|
|
v <- unsafeRead cmpsts i
|
|
if v then loop (i + 1) cbpi rlmti else do
|
|
unsafeWrite cmpsts i True -- cull current bp so not a "k-rough"!
|
|
let bp = i + i + 1
|
|
-- partial cull by current base prime, bp...
|
|
cull c = if c > mxndx then return () else do
|
|
unsafeWrite cmpsts c True; cull (c + bp)
|
|
|
|
-- adjust `larges` according to partial sieve...
|
|
part ri nri = -- old "rough" index to new one...
|
|
if ri > rlmti then return (nri - 1) else do
|
|
r <- unsafeRead mrs ri -- "rough" always odd!
|
|
t <- unsafeRead cmpsts (fromIntegral r `shiftR` 1)
|
|
if t then part (ri + 1) nri else do -- skip newly culled
|
|
olv <- unsafeRead mls ri
|
|
let m = fromIntegral r * fromIntegral bp
|
|
adjv <- if m <= fromIntegral rtlmt then do
|
|
let ndx = fromIntegral m `shiftR` 1
|
|
sv <- unsafeRead mss ndx
|
|
unsafeRead mls (fromIntegral sv - cbpi)
|
|
else do
|
|
sv <- unsafeRead mss (half (divide n m))
|
|
return (fromIntegral sv)
|
|
unsafeWrite mls nri (olv - (adjv - fromIntegral cbpi))
|
|
unsafeWrite mrs nri r; part (ri + 1) (nri + 1)
|
|
!pm0 = ((rtlmt `div` bp) - 1) .|. 1 -- max base prime mult
|
|
|
|
adjc lmti pm = -- adjust smalls according to partial sieve:
|
|
if pm < bp then return () else do
|
|
c <- unsafeRead mss (pm `shiftR` 1)
|
|
let ac = c - fromIntegral cbpi -- correction
|
|
bi = (pm * bp) `shiftR` 1 -- start array index
|
|
adj si = if si > lmti then adjc (bi - 1) (pm - 2)
|
|
else do ov <- unsafeRead mss si
|
|
unsafeWrite mss si (ov - ac)
|
|
adj (si + 1)
|
|
adj bi
|
|
dadjc lmti pm =
|
|
if pm < bp then return () else do
|
|
c <- unsafeRead mss (pm `shiftR` 1)
|
|
let ac = c - fromIntegral cbpi -- correction
|
|
bi = (pm * bp) `shiftR` 1 -- start array index
|
|
ac64 = fromIntegral ac :: Int64
|
|
dac = (ac64 `shiftL` 32) .|. ac64
|
|
dbi = (bi + 1) `shiftR` 1
|
|
dlmti = (lmti - 1) `shiftR` 1
|
|
dadj dsi = if dsi > dlmti then return ()
|
|
else do dov <- unsafeRead mdss dsi
|
|
unsafeWrite mdss dsi (dov - dac)
|
|
dadj (dsi + 1)
|
|
when (bi .&. 1 /= 0) $ do
|
|
ov <- unsafeRead mss bi
|
|
unsafeWrite mss bi (ov - ac)
|
|
dadj dbi
|
|
when (lmti .&. 1 == 0) $ do
|
|
ov <- unsafeRead mss lmti
|
|
unsafeWrite mss lmti (ov - ac)
|
|
adjc (bi - 1) (pm - 2)
|
|
cull sqri; nrlmti <- part 0 0
|
|
dadjc mxndx pm0
|
|
loop (i + 1) (cbpi + 1) nrlmti
|
|
loop 1 0 mxndx
|
|
|
|
!ans0 = unsafeAt larges 0 - -- combine all counts; each includes nbps...
|
|
sum [ unsafeAt larges i | i <- [ 1 .. nrs - 1 ] ]
|
|
-- adjust for all the base prime counts subracted above...
|
|
!adj = (nrs + 2 * (nbps - 1)) * (nrs - 1) `div` 2
|
|
!adjans0 = ans0 + fromIntegral adj
|
|
|
|
loopr ri !acc = -- o final phi calculation for pairs of larger primes...
|
|
let r = fromIntegral (unsafeAt roughs ri)
|
|
q = n `div` r
|
|
lmtsi = half (fromIntegral (q `div` r))
|
|
lmti = fromIntegral (unsafeAt smalls lmtsi) - nbps
|
|
addcnt pi !ac =
|
|
if pi > lmti then ac else
|
|
let p = fromIntegral (unsafeAt roughs pi)
|
|
ci = half (fromIntegral (divide q p))
|
|
in addcnt (pi + 1) (ac + fromIntegral (unsafeAt smalls ci))
|
|
in if lmti <= ri then acc else -- break when up to cube root of range!
|
|
-- adjust for the `nbps`'s over added in the `smalls` counts...
|
|
let !adj = fromIntegral ((lmti - ri) * (nbps + ri - 1))
|
|
in loopr (ri + 1) (addcnt (ri + 1) acc - adj)
|
|
in loopr 1 adjans0 + 1 -- add one for only even prime of two!
|