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!