import Control.Monad.ST ( runST, ST )
import Data.Array.Base ( MArray(newArray, unsafeRead, unsafeWrite),
                         IArray(unsafeAt),
                         STUArray, unsafeFreezeSTUArray, assocs )
import Data.Time.Clock.POSIX ( getPOSIXTime ) 

primesTo :: Int -> [Int] 
primesTo limit = runST $ do
  let lmt = limit - 2
  cmpsts <- newArray (2, limit) False 
  cmpstsf <- unsafeFreezeSTUArray cmpsts 
  let getbpndx bp = (bp, bp * bp - 2) 
      cullcmpst i = unsafeWrite cmpsts i True 
      cull4bpndx (bp, si0) = mapM_ cullcmpst [ si0, si0 + bp .. lmt ]
  mapM_ cull4bpndx
        $ takeWhile ((>=) lmt . snd) 
                    [ getbpndx bp | (bp, False) <- assocs cmpstsf ]
  return [ p | (p, False) <- assocs cmpstsf ] 


main :: IO ()
main = do
  putStrLn $ "The primes up to 100 are " ++ show (primesTo 100)
  putStrLn $ "The number of primes up to a million is " ++
               show (length $ primesTo 1000000)
  let top = 1000000000
  start <- getPOSIXTime
  let answr = length $ primesTo top
  stop <- answr `seq` getPOSIXTime 
  let elpsd =  round $ 1e3 * (stop - start) :: Int

  putStrLn $ "Found " ++ show answr ++ " to " ++ show top ++
               " in " ++ show elpsd ++ " milliseconds."