Performance improvements with math

.vscode/settings.json

@@ -0,0 +1,7 @@
+{
+  "cSpell.words": [
+    "concat",
+    "coprime",
+    "foldl"
+  ]
+}

flake.nix

@@ -45,6 +45,14 @@
             aeson
             random
             neat-interpolation
+
+            # maths
+            primes
+            arithmoi
+
+            # graphing libraries!
+            Chart
+            Chart-cairo
           ]
         );
       in

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]

src/projectEuler/question11.hs

@@ -0,0 +1,110 @@
+{-
+In the 20x20 grid below, four numbers along a diagonal line have been marked in *red*.
+
+08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08
+49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00
+81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65
+52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91
+22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80
+24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50
+32 98 81 28 64 23 67 10*26*38 40 67 59 54 70 66 18 38 64 70
+67 26 20 68 02 62 12 20 95*63*94 39 63 08 40 91 66 49 94 21
+24 55 58 05 66 73 99 26 97 17*78*78 96 83 14 88 34 89 63 72
+21 36 23 09 75 00 76 44 20 45 35*14*00 61 33 97 34 31 33 95
+78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92
+16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57
+86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58
+19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40
+04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66
+88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69
+04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36
+20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16
+20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54
+01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48
+
+The product of these numbers is 26x63x78x14 = 1788696.
+
+What is the greatest product of four adjacent numbers in the same direction (up, down, left, right, or diagonally) in the 20x20 grid?
+-}
+import           Data.List (genericDrop, genericIndex, genericTake)
+
+grid :: [[Integer]]
+grid =
+  [ [08, 02, 22, 97, 38, 15, 00, 40, 00, 75, 04, 05, 07, 78, 52, 12, 50, 77, 91, 08]
+  , [49, 49, 99, 40, 17, 81, 18, 57, 60, 87, 17, 40, 98, 43, 69, 48, 04, 56, 62, 00]
+  , [81, 49, 31, 73, 55, 79, 14, 29, 93, 71, 40, 67, 53, 88, 30, 03, 49, 13, 36, 65]
+  , [52, 70, 95, 23, 04, 60, 11, 42, 69, 24, 68, 56, 01, 32, 56, 71, 37, 02, 36, 91]
+  , [22, 31, 16, 71, 51, 67, 63, 89, 41, 92, 36, 54, 22, 40, 40, 28, 66, 33, 13, 80]
+  , [24, 47, 32, 60, 99, 03, 45, 02, 44, 75, 33, 53, 78, 36, 84, 20, 35, 17, 12, 50]
+  , [32, 98, 81, 28, 64, 23, 67, 10, 26, 38, 40, 67, 59, 54, 70, 66, 18, 38, 64, 70]
+  , [67, 26, 20, 68, 02, 62, 12, 20, 95, 63, 94, 39, 63, 08, 40, 91, 66, 49, 94, 21]
+  , [24, 55, 58, 05, 66, 73, 99, 26, 97, 17, 78, 78, 96, 83, 14, 88, 34, 89, 63, 72]
+  , [21, 36, 23, 09, 75, 00, 76, 44, 20, 45, 35, 14, 00, 61, 33, 97, 34, 31, 33, 95]
+  , [78, 17, 53, 28, 22, 75, 31, 67, 15, 94, 03, 80, 04, 62, 16, 14, 09, 53, 56, 92]
+  , [16, 39, 05, 42, 96, 35, 31, 47, 55, 58, 88, 24, 00, 17, 54, 24, 36, 29, 85, 57]
+  , [86, 56, 00, 48, 35, 71, 89, 07, 05, 44, 44, 37, 44, 60, 21, 58, 51, 54, 17, 58]
+  , [19, 80, 81, 68, 05, 94, 47, 69, 28, 73, 92, 13, 86, 52, 17, 77, 04, 89, 55, 40]
+  , [04, 52, 08, 83, 97, 35, 99, 16, 07, 97, 57, 32, 16, 26, 26, 79, 33, 27, 98, 66]
+  , [88, 36, 68, 87, 57, 62, 20, 72, 03, 46, 33, 67, 46, 55, 12, 32, 63, 93, 53, 69]
+  , [04, 42, 16, 73, 38, 25, 39, 11, 24, 94, 72, 18, 08, 46, 29, 32, 40, 62, 76, 36]
+  , [20, 69, 36, 41, 72, 30, 23, 88, 34, 62, 99, 69, 82, 67, 59, 85, 74, 04, 36, 16]
+  , [20, 73, 35, 29, 78, 31, 90, 01, 74, 31, 49, 71, 48, 86, 81, 16, 23, 57, 05, 54]
+  , [01, 70, 54, 71, 83, 51, 54, 69, 16, 92, 33, 48, 61, 43, 52, 01, 89, 19, 67, 48]
+  ]
+
+points :: Integer -> [(Integer, Integer)]
+points n =  [ (x,y) | x <- [1..n], y <- [1..n] ]
+
+allDiagonals :: Integer -> [[(Integer, Integer)]]
+allDiagonals n =
+  [ [(i, j), b, c, d]
+  |  (i, j) <- points n
+  ,  i+3 <= n
+  ,  j+3 <= n
+  ,  b <- [(i + 1, j + 1)]
+  ,  c <- [(i + 2, j + 2)]
+  ,  d <- [(i + 3, j + 3)]
+  ]
+
+allUpDown :: Integer -> [[(Integer, Integer)]]
+allUpDown n =
+  [ [(i, j), b, c, d]
+  |  (i, j) <- points n
+  ,  i   <= n
+  ,  j+3 <= n
+  ,  b <- [(i, j + 1)]
+  ,  c <- [(i, j + 2)]
+  ,  d <- [(i, j + 3)]
+  ]
+
+allLeftRight :: Integer -> [[(Integer, Integer)]]
+allLeftRight n =
+  [ [(i, j), b, c, d]
+  |  (i, j) <- points n
+  ,  i+3 <= n
+  ,  j   <= n
+  ,  b <- [(i + 1, j)]
+  ,  c <- [(i + 2, j)]
+  ,  d <- [(i + 3, j)]
+  ]
+
+allSequences :: Integer -> [[(Integer, Integer)]]
+allSequences n = allDiagonals n ++ allUpDown n ++ allLeftRight n
+
+s20x20 :: [[(Integer, Integer)]]
+s20x20 = allSequences 20
+
+getValue :: [[Integer]] -> (Integer, Integer) -> Integer
+getValue grid (x,y) = row `genericIndex` y'
+  where row = (head . genericDrop (x-1) . genericTake x) grid
+        y'  = y-1 -- compensate for 1 index numbering above
+
+-- |
+-- λ> ans
+-- 51267216
+-- (0.13 secs, 38,281,312 bytes)
+ans :: Integer
+ans = maximum $ map (product . map (getValue grid)) s20x20
+
+main :: IO ()
+main = print ans

src/projectEuler/question12.hs

@@ -0,0 +1,32 @@
+-- module Main where
+import           Data.IntSet                           (size)
+import           Data.List                             (unfoldr)
+import           Math.NumberTheory.ArithmeticFunctions (divisorsSmall)
+
+main :: IO ()
+main = print ans
+
+-- triangleNumber :: Integer -> [Integer]
+-- triangleNumber 1 = [1]
+-- triangleNumber n = sum [1..n] : triangleNumber (n-1)
+triangleNumber i = (i * (i + 1)) `div` 2
+triangleNumbers = map triangleNumber [1..]
+-- triangleNumbers :: [Int]
+-- triangleNumbers = unfoldr (\b-> Just (sum [1..b],b+1)) 1
+
+-- isDivisible a b = a `rem` b == 0
+-- divisors n = filter (isDivisible n) [1..(floor . sqrt) n ]
+
+-- Triangle Numbers With More than N divisors
+tnwmtNd :: Int -> [Int]
+tnwmtNd n = filter ((>= n) . size . divisorsSmall) triangleNumbers
+
+ans :: Int
+ans = head $ tnwmtNd 500
+{-
+λ> main
+76576500
+(0.00 secs, 114,944 bytes)
+
+First triangle number to have over 500 divisors is 76576500
+-}

src/projectEuler/question7.hs

@@ -0,0 +1,32 @@
+-- https://projecteuler.net/problem=7
+-- Find the 10_001 prime number
+
+primes1 :: [Integer]
+primes1 = 2:3:prs
+  where
+    1:p:candidates = [6*k+r | k <- [0..], r <- [1,5]]
+    prs            = p : filter isPrime candidates
+    isPrime n      = not (any (divides n) (takeWhile (\p -> p*p <= n) prs))
+    divides n p    = n `mod` p == 0
+
+
+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
+
+
+primes3 = sieve [2..]
+  where
+    sieve (p : xs) = p : sieve [x | x <- xs, x `mod` p > 0]

src/projectEuler/question8.hs

@@ -0,0 +1,66 @@
+{-
+The four adjacent digits in the 1000-digit number that have the greatest product are 9 x 9 x 8 x 9 = 5832.
+
+73167176531330624919225119674426574742355349194934
+96983520312774506326239578318016984801869478851843
+85861560789112949495459501737958331952853208805511
+12540698747158523863050715693290963295227443043557
+66896648950445244523161731856403098711121722383113
+62229893423380308135336276614282806444486645238749
+30358907296290491560440772390713810515859307960866
+70172427121883998797908792274921901699720888093776
+65727333001053367881220235421809751254540594752243
+52584907711670556013604839586446706324415722155397
+53697817977846174064955149290862569321978468622482
+83972241375657056057490261407972968652414535100474
+82166370484403199890008895243450658541227588666881
+16427171479924442928230863465674813919123162824586
+17866458359124566529476545682848912883142607690042
+24219022671055626321111109370544217506941658960408
+07198403850962455444362981230987879927244284909188
+84580156166097919133875499200524063689912560717606
+05886116467109405077541002256983155200055935729725
+71636269561882670428252483600823257530420752963450
+
+Find the thirteen adjacent digits in the 1000-digit number that have the greatest product.
+What is the value of this product?
+-}
+module Main where
+
+import           Data.Char       (ord)
+import           Data.List       (sort)
+import           Data.List.Split (splitOn)
+
+thousandDigitNum :: Integer
+thousandDigitNum = 7316717653133062491922511967442657474235534919493496983520312774506326239578318016984801869478851843858615607891129494954595017379583319528532088055111254069874715852386305071569329096329522744304355766896648950445244523161731856403098711121722383113622298934233803081353362766142828064444866452387493035890729629049156044077239071381051585930796086670172427121883998797908792274921901699720888093776657273330010533678812202354218097512545405947522435258490771167055601360483958644670632441572215539753697817977846174064955149290862569321978468622482839722413756570560574902614079729686524145351004748216637048440319989000889524345065854122758866688116427171479924442928230863465674813919123162824586178664583591245665294765456828489128831426076900422421902267105562632111110937054421750694165896040807198403850962455444362981230987879927244284909188845801561660979191338754992005240636899125607176060588611646710940507754100225698315520005593572972571636269561882670428252483600823257530420752963450
+
+digits :: Integer -> [Integer]
+digits = go []
+  where
+    go acc 0 = acc
+    go acc n = go (digit:acc) rest
+      where
+        digit = n `rem` 10
+        rest  = n `div` 10
+
+main :: IO ()
+main = print ans
+
+ans :: [[Integer]]
+ans = (windowsOf 13 . digits) thousandDigitNum
+-- "9878799272442"
+
+-- what dose it mean to have the greatest product in adjacent digits?
+--  why ask it like that???
+--  doesn't have 0
+--  otherwise is the largest number in the larger 1000-digit number
+chunksWithoutZero :: Int -> [String]
+chunksWithoutZero n = (filter (\str -> length str >= n) . filter (/= "") . splitOn "0" . show) thousandDigitNum
+
+windowsOf :: Int -> [a] -> [[a]]
+windowsOf n = drop (n-1) . go []
+  where
+    go acc [] = acc
+    go acc lst = go (lst':acc) (tail lst)
+      where
+        lst' = take n lst

src/projectEuler/question9.hs

@@ -0,0 +1,105 @@
+import           Debug.Trace                            (trace)
+import           Graphics.Rendering.Chart.Backend.Cairo
+import           Graphics.Rendering.Chart.Easy
+{-
+Special Pythagorean Triplet
+Problem 9
+
+A Pythagorean triplet is a set of three natural numbers, `a < b < c`, for which,
+
+`a^2 + b^2 = c^2`.
+
+For example, `3^2 + 4^2 = 9 + 16 = 25 = 5^2`.
+
+There exists exactly one Pythagorean triplet for which `a + b + c = 1000`.
+
+Find the product `abc`.
+-}
+
+main :: IO ()
+main = print answer
+
+answer :: String
+answer = "I dunno"
+
+
+solve :: Integer -> [(Integer, Integer, Integer)]
+solve x = takeWhile (\(a,b,c) -> a + b + c <= 1000) $ primitiveTriplesUnder x
+
+euclid'sFormula :: Num c => (c, c) -> (c, c, c)
+euclid'sFormula (m, n) = (a,b,c)
+  where
+    a = m^2 - n^2
+    b = 2*m*n
+    c = m^2 + n^2
+
+listOfMNs :: Integer -> [(Integer, Integer)]
+listOfMNs x =
+  [ (m,n)
+      | n <- [2,4..x]     -- one of them is even
+      , m <- [n+1,n+3..x]
+      , gcd m n == 1      -- coprime
+  ]
+
+listOfMNs' :: Integer -> [(Integer, Integer)]
+listOfMNs' x =
+  [ (m,n)
+      | n <- [2,4..]     -- one of them is even
+      , m <- [n+1,n+3..]
+      , gcd m n == 1     -- coprime
+      , a m n + b m n + c m n <= x
+  ] where
+      a m n = m^2 - n^2
+      b m n = 2*m*n
+      c m n = m^2 + n^2
+
+
+primitiveTriplesUnder :: Integer -> [(Integer, Integer, Integer)]
+primitiveTriplesUnder = map euclid'sFormula . listOfMNs
+
+
+test :: [(Integer, Integer, Integer)]
+test = [ (a,b,c)
+          | a <- [3..],
+            b <- take 10 [a+1..],
+            c <- takeWhile (\c -> a^2 + b^2 <= c^2) [b+1..]
+       ]
+
+ls :: [(Integer, Integer)]
+ls = filter (\(m,n) -> gcd m n == 1) $ zip [5,9..] [2,4..]
+
+diags :: Integer -> [(Integer, Integer)]
+diags n = [(x,y) | x<-[0..n], y<-[0..n]]
+
+
+graph :: IO ()
+graph = toFile def "test.png" $ do
+    layout_title .= "Points"
+    plot (line "points" [ [ (x,y) | (x,y) <- listOfMNs 30] ])
+
+
+-- main'' = print ans'
+ans' :: Integer -> [(Integer, Integer, Integer)]
+ans' limit = [(a, b, c)
+              | a <- [1 .. limit]
+              , b <- [a + 1 .. limit]
+              , c <- [limit - a - b]
+              , b < c
+             ]
+
+{- 
+ - Solution to Project Euler problem 9
+ - Copyright (c) Project Nayuki. All rights reserved.
+ - 
+ - https://www.nayuki.io/page/project-euler-solutions
+ - https://github.com/nayuki/Project-Euler-solutions
+ -}
+
+
+-- Computers are fast, so we can implement a brute-force search to directly solve the problem.
+perim = 1000
+main = putStrLn (show ans)
+ans = head [a * b * (perim - a - b) | a <- [1..perim], b <- [a+1..perim], isIntegerRightTriangle a b]
+isIntegerRightTriangle a b = a < b && b < c
+	&& a * a + b * b == c * c
+	where c = perim - a - b