diff --git a/flake.nix b/flake.nix index 9e55396..aeb09a1 100644 --- a/flake.nix +++ b/flake.nix @@ -46,6 +46,10 @@ random neat-interpolation + # maths + primes + arithmoi + # graphing libraries! Chart Chart-cairo diff --git a/src/projectEuler/question10.hs b/src/projectEuler/question10.hs new file mode 100644 index 0000000..1c5b8d5 --- /dev/null +++ b/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]