{-
Summation of Primes
Problem 10

The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17.

Find the sum of all the primes below two million.


https://wiki.haskell.org/Prime_numbers

Compile with `ghc -O3 src/projectEuler/question10.hs`
Run with `time src/projectEuler/question10`
-}
{-# OPTIONS_GHC -O2 #-}

import           Control.Monad            (forM_, when)
import           Data.Array.Base          (IArray (unsafeAt),
                                           MArray (unsafeWrite),
                                           UArray (UArray),
                                           unsafeFreezeSTUArray)
import           Data.Array.ST            (MArray (newArray), readArray,
                                           runSTUArray, writeArray)
import           Data.Array.Unboxed       (UArray, assocs)
import           Data.Bits                (Bits (shiftL, shiftR))
import           Data.Int                 (Int64)
import           Data.List                (genericIndex, genericTake)
import           Data.Word                (Word64)
-- https://hackage.haskell.org/package/arithmoi-0.13.0.0
import           Math.NumberTheory.Primes (Prime (unPrime), primes)

main :: IO ()
main = do
  print "Hello! Welcome to the prime number generator."
  -- print "Please enter which nth prime you'd like: "
  -- n <- getLine
  -- let n' = read n :: Integer
  -- print $ "Finding the " ++ n ++ "th prime."
  -- print $ primes() `genericIndex` (10^ (10 :: Int))
  print $ last $ take (10^ (10 :: Int)) primes

ans :: Integer
ans = sum $ takeWhile (< 2_000_000) $ map unPrime primes


primes1 :: [Integer]
primes1 = sieve [2..]
  where
    sieve (p : xs) = p : sieve [x | x <- xs, x `mod` p > 0]

primes2 :: [Integer]
primes2 = 2:([3..] `minus` composites)
  where
    composites = union [multiples p | p <- primes2]

    multiples n = map (n*) [n..]

    (x:xs) `minus` (y:ys) | x < y = x:(xs `minus` (y:ys))
                          | x == y = xs `minus` ys
                          | x > y = (x:xs) `minus` ys
    union = foldr merge [ ]
      where
        merge (x:xs) ys = x:merge' xs ys
        merge' (x:xs) (y:ys) | x < y = x:merge' xs (y:ys)
                            | x == y = x:merge' xs ys
                            | x > y = y:merge' (x:xs) ys


{-
https://wiki.haskell.org/Prime_numbers
Using Page-Segmented ST-Mutable Unboxed Array
-}
-- type Prime = Word64

-- cSieveBufferLimit :: Int
-- cSieveBufferLimit = 2^17 * 8 - 1 -- CPU L2 cache in bits

-- primes :: () -> [Prime]
-- primes() = 2 : _Y (listPagePrms . pagesFrom 0) where
--   _Y g = g (_Y g) -- non-sharing multi-stage fixpoint combinator
--   listPagePrms pgs@(hdpg@(UArray lwi _ rng _) : tlpgs) =
--     let loop i | i >= fromIntegral rng = listPagePrms tlpgs
--                | unsafeAt hdpg i = loop (i + 1)
--                | otherwise = let ii = lwi + fromIntegral i in
--                              case fromIntegral $ 3 + ii + ii of
--                                p -> p `seq` p : loop (i + 1) in loop 0
--   makePg lwi bps = runSTUArray $ do
--     let limi = lwi + fromIntegral cSieveBufferLimit
--         bplmt = floor $ sqrt $ fromIntegral $ limi + limi + 3
--         strta bp = let si = fromIntegral $ (bp * bp - 3) `shiftR` 1
--                    in if si >= lwi then fromIntegral $ si - lwi else
--                    let r = fromIntegral (lwi - si) `mod` bp
--                    in if r == 0 then 0 else fromIntegral $ bp - r
--     cmpsts <- newArray (lwi, limi) False
--     fcmpsts <- unsafeFreezeSTUArray cmpsts
--     let cbps = if lwi == 0 then listPagePrms [fcmpsts] else bps
--     forM_ (takeWhile (<= bplmt) cbps) $ \ bp ->
--       forM_ (let sp = fromIntegral $ strta bp
--              in [ sp, sp + fromIntegral bp .. cSieveBufferLimit ]) $ \c ->
--         unsafeWrite cmpsts c True
--     return cmpsts
--   pagesFrom lwi bps = map (`makePg` bps)
--                           [ lwi, lwi + fromIntegral cSieveBufferLimit + 1 .. ]




-- sieveUA :: Int -> UArray Int Bool
-- sieveUA top = runSTUArray $ do
--     let m = (top-1) `div` 2
--         r = floor . sqrt $ fromIntegral top + 1
--     sieve <- newArray (1,m) True      -- :: ST s (STUArray s Int Bool)
--     forM_ [1..r `div` 2] $ \i -> do
--       isPrime <- readArray sieve i
--       when isPrime $ do               -- ((2*i+1)^2-1)`div`2 == 2*i*(i+1)
--         forM_ [2*i*(i+1), 2*i*(i+2)+1..m] $ \j -> do
--           writeArray sieve j False
--     return sieve

-- primesToUA :: Int -> [Int]
-- primesToUA top = 2 : [i*2+1 | (i,True) <- assocs $ sieveUA top]