Euler Question 10 - Whole lotta primes

flake.nix

@@ -46,6 +46,10 @@
             random
             neat-interpolation
 
+            # maths
+            primes
+            arithmoi
+
             # graphing libraries!
             Chart
             Chart-cairo

src/projectEuler/question10.hs

@@ -0,0 +1,121 @@
+{-
+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]