Saturday, December 21, 2013

Solving Some Project Euler Problems in Haskell

Here are a few solutions to Project Euler problems done in Haskell.  All of these were entered in gchi interactively.  As much as possible I try to keep them as one-liners.

The ">" indicates the beginning of a line at the ghci command line.

Starting with problem 7, this ideal could not be kept since it was becoming unmanageable.  With that problem I started to use modules.  The maintenance of that module is documented below.

One final note; this post is an experiment at producing fewer posts but keeping them alive with data and updates rather than flood with posts.  Quality over quantity seems to be my new mantra.  A good long post completed over the course of a few months.

Problem 1:

This problem was relatively straight-forward.  It's actually a nice one-liner!
> sum $ filter (\x -> x `mod` 3 == 0 || x `mod` 5 == 0) [1..999]

Problem 2:

First, create a fibonacci sequence list.  Make sure it is tail recursive for performance - hence why the first two elements are not included in the recursive list.  The filter and takeWhile should be interchangeable.
> let fib x y = (x+y):(fib y (x+y))> sum $ filter (\x -> x `mod` 2 == 0) $ takeWhile (<= 4000000) ([1,2] ++ fib 1 2)

Problem 3:

A much more interesting problem.  Note, the maximum value of your prime number will be the square root of the target number.

You do not need to enumerate all primes.  Actually, this would be slower as the sieve would take more time filtering numbers.  To crack the number, we first find all integers that divide it.  Then, of those integers, to determine which are prime, we see which divide each other.

To see why this works, consider an integer and all the integers that divide it.  If an integer were not prime, it could be decomposed into prime numbers.  Those prime numbers also divide the original number.

First identifying the dividing the integers being faster is only tractable for relatively small numbers.

"divs" finds all the integers that divide the given value.
"noDivs" is true if there is a value larger than "x" in "xs" which is divisible by "x".
"divs'" combines the previous two functions to find all the primes for the given value.
> let noDivs x xs = (length $ filter (\y -> (x `mod` y == 0) && (x > y)) xs) == 0> let divs x = filter (\y -> x `mod` y == 0) [2..(ceiling $ sqrt $ fromIntegral x)]> let divs' x = filter (\y -> noDivs y z) z where z = divs xmaximum $ divs' 600851475143


Problem 4:

Rather simple problem which can be solved in a single line.  Simply import Data.List:
> filter (\x -> reverse x == x) $ map show $ reverse $ sort $ concat $ map (\x -> map (*x) [1..x] ) [1..999]

Problem 5:

So, let's start with the wrong solution.  We could brute-force the problem as so:
> head $ filter snd $ map (\x -> (x,foldl (&&) True $ map (==0) $ map (\y -> x `mod` y) [2..20])) [100..]

We need a better solution.  Let's suppose we knew the prime factors of the numbers that we're dealing with.  Like so:
> let primesUntil p = filter (\x -> length (filter (\y -> x `mod` y == 0) $ [2..x-1]) == 0 ) [2..p]

Now, we can think of the solution as a permutation of the primes that arise from numbers 0...20.  Then, let us suppose we also had a function to enumerate all the possible permutations:
> let combos v n c = zipWith (\x y -> x + ((v `div` y) `mod` c)) (replicate n 1) $ map (\x -> c^x) [0..n-1]> let primeCombos p = take (3^(length z)) $ map (\x -> foldl (*) 1 $ zipWith (*) (combos x (length z) 3) z)[0..] where z = (primesUntil p)

Finally, we can combine everything together to get:
> minimum $ map fst $ filter snd $ map (\x -> (x,foldl (&&) True $ map (==0) $ map (\y -> x `mod` y) [2..20])) $ primeCombos 20

Not ideal as a solution though.  It takes quite a bit of time to run.  What would be ideal?  Proper analysis of the problems in terms of primes.  The smallest product is a product of the primes of 0...20.  Isolating the primes is a problem that would have added excess code though.

Problem 6:

Something that Haskell excels at, thankfully due to built-in big number support:
> let diffSq x = (foldl (+) 0 x)^2 - (foldl (\z y -> z + y*y) 0 x)> diffSq [1..100]

Problem 7:

For this problem, we will need a very fast function to calculate primes.  Very fast.  The naive solution though would be (assuming we have primesUntil from above):
> primesUntil 1000000000 !! 10000

Notice how the function doesn't return.  primesUntil breaks down when we wish to compute large primes since we test against all previous numbers and not just known primes.  A good prime function will use previous primes to find new primes.  It will also discard even numbers from the start.

Let me cook up a few little functions to help...

isDividedBy is the first.  It checks to see if any value divides any other value in a list.  It is intended to be written as: "10 isDividedBy [3, 4, 7]".  Note that foldr is more apt to lazy evaluation, which is exactly what we seek.
isDividedBy             :: Integer -> [Integer] -> Bool
isDividedBy x xs        =  foldr (\a b -> (x `mod` a == 0) || b) False xs


We desire a slightly more efficient way to enumerate primes.  Note that no prime will be even since they are all divisible by 2.  Also, primes become very sparse as we progress through the number line so discarding them quickly becomes kind of important.  The first primes will be multiples of a number more often than the latter primes.  Let us keep that in note and define a function to enumerate future primes:
nextPrime               :: Integer -> [Integer] -> Integer
nextPrime n xs          =  if (n `isDividedBy` xs)
                           then (nextPrime (n+2) xs)
                           else n
Good!  Let us suppose we had a way to perpetually enumerate future primes:
nextPrimes              :: Integer -> [Integer] -> [Integer]
nextPrimes x xs         =  np:(nextPrimes (x+2) pl)
                           where                              np = nextPrime x xs
                              pl = xs ++ [np]
Problem?  Still doesn't perform as well as needed.  How can we speed this code up?  There are primality tests which we can apply to filter numbers rather rapidly.  However, even with Fermat's Primality Test, we do not get much speed.
quickFermat             :: Integer -> Bool
quickFermat p           =  a^(p-1) `mod` p == 1                           where a = 2 -- p `div` 2
nextPrime               :: Integer -> [Integer] -> Integer
nextPrime n xs          =  if ((quickFermat n) || (n `isDividedBy` xs))
                           then (nextPrime (n+2) xs)
                           else n



Sounds like it's time to rethink our algorithm.  We know that division is expensive, so we should build a sieve.  Sieve of Atkin (according to Wikipedia) is the modern-day best-case.  In such a case, we wish to mark off numbers that we know to not be prime.

We can use Bertrand's Postulate and enumerate a set of values that we believe are prime between the intervals.  I would rather point out the idea of diminishing returns and mention that '2' cuts out half the numbers from primality testing.  '3' cuts out at most a quarter of the numbers.  '5' less, '7' even less, etc.  The problem is that we aren't testing just highly potential primes but the entire number line.

Let us try to enumerate all numbers that are multiples of the first few primes...  And Eureka!  A much better solution has emerged!  Unfortunately, it does not have the runtime that I'd desire.

The idea is simple, build a list of known composite numbers and flip it to get anything that is prime.  Yes - it is another implementation of the sieve of Eratosthenes.  Where it excels is at reducing the number of divisions.  There is a list being built that does not contain any prime.  By comparing it with a list of all numbers, we get a list of potential primes.  From that list, we can do simple primality tests.

Here is the code, in completion:
import Data.List
knownComposites         :: Integer -> Integer -> Integer -> [Integer]
knownComposites p s m   =  if ( m > 0 ) then p:(knownComposites (p+s) s (m-s))
                           else []
mergeLists              :: [Integer] -> [Integer] -> [Integer]
mergeLists [] []        = []
mergeLists [] b         =  b
mergeLists a []         =  a
mergeLists (a:aa) (b:bb)=  if (a == b) then a:(mergeLists aa bb)
                              else if (a<b) then a:(mergeLists aa (b:bb))
                                       else b:(mergeLists (a:aa) bb)
mergeComposites         :: [Integer] -> Integer -> [Integer]
mergeComposites [] _    = []
mergeComposites (p:pp)m = mergeLists (knownComposites p p m) (mergeComposites pp m)
compositesForPrimes     :: [Integer] -> [Integer]
compositesForPrimes p   =  mergeComposites p (foldr (*) 1 p)
missingFrom                :: [Integer] -> [Integer] -> [Integer]
missingFrom [] _           = []
missingFrom (a:aa) (b:bb)  = if (a > b) then b:(missingFrom (a:aa) bb)
                              else missingFrom aa bb
potentialPrimes         :: [Integer] -> [Integer]
potentialPrimes p       =  missingFrom (compositesForPrimes p) [2..]
quickFermat             :: Integer -> Bool
quickFermat p           =  a^(p-1) `mod` p == 1                           where a = 2 -- p `div` 2
isDividedBy             :: Integer -> [Integer] -> Bool
isDividedBy _ []        = False
isDividedBy x xs        =  foldr (\a b -> (x `mod` a == 0) || b) False xs
nextPrime               :: Integer -> [Integer] -> Integer
nextPrime n xs          =  if ((quickFermat n) || (n `isDividedBy` xs))
                           then (nextPrime (n+2) xs)
                           else n
nextPrimes              :: Integer -> [Integer] -> [Integer]
nextPrimes x xs         =  np:(nextPrimes (x+2) pl)
                           where                              np = nextPrime x xs
                              pl = xs ++ [np]
allPrimes               :: [Integer]
allPrimes               = 1:2:3:5:(nextPrimes 5 [2,3,5])
primeFilter'            :: [Integer] -> [Integer] -> [Integer]
primeFilter' p (m:mm)   =  if (m `isDividedBy` p)
                           then primeFilter' p mm
                           else m:(primeFilter' (p ++ [m]) mm)
allPrimes'              = (1:ip) ++ (primeFilter' [] $ potentialPrimes ip)
                        where ip = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233]

Just run the command "allPrimes' !! 10002" to get the answer.  I counted 1 as a prime.

Problem 8:

This one is a cute problem.  We'll solve it by first defining a variable called number by manually concatenating the given number together onto one line:
> let number = "7316717653133062491922511967442657474235534919493496983520312774506326239578318016984801869478851843858615607891129494954595017379583319528532088055111254069874715852386305071569329096329522744304355766896648950445244523161731856403098711121722383113622298934233803081353362766142828064444866452387493035890729629049156044077239071381051585930796086670172427121883998797908792274921901699720888093776657273330010533678812202354218097512545405947522435258490771167055601360483958644670632441572215539753697817977846174064955149290862569321978468622482839722413756570560574902614079729686524145351004748216637048440319989000889524345065854122758866688116427171479924442928230863465674813919123162824586178664583591245665294765456828489128831426076900422421902267105562632111110937054421750694165896040807198403850962455444362981230987879927244284909188845801561660979191338754992005240636899125607176060588611646710940507754100225698315520005593572972571636269561882670428252483600823257530420752963450"
Good!  Good!  Now, let us split up this string into groups of 5 digits!
fives                   :: [x] -> [[x]]
fives v                 = if ((length v) < 5) then ([])
                        else ((take 5 v):(fives $ drop 1 v))

Wonderful!  Just one last part...
> import Data.Char
> maximum $ map (\x -> foldr (*) 1 $ map digitToInt x) $ fives number
What am I doing?  Taking the big string called "number" and splitting it into substrings of "5" in length (look at the definition of fives).  Then I'm mapping this lambda to each of the strings that multiplies the digits of the number together.

The "map digitToInt" will convert a string of character digits to a list of numbers.  "234" becomes [2,3,4].  The foldr multiplies it all together.  This is encapsulated in a map which applies the operation to the list from the fives function. 

Finally, maximum gives us our solution.

Problem 9:

Finally!  Another one-liner!
> map (\(x,y) -> x*y*(1000-x-y)) $ filter (\(x,y) -> x*x + y*y == (1000-x-y)*(1000-x-y) ) [(i,j) | i <- [1..1000], j <- [1..1000-i]]

There!  Now read it backwards!  We generate a tuple of numbers.  The last one must be chosen so that the sum is 1000.  Filter out any value that does not match the pythagorean inequality, and then map the remaining values to the product.  And voila!  Done!

Problem 10:

Ah!  Another problem involving prime numbers!  And we must enumerate an even greater set of primes than before!  As a hint, we wish to be able to write something as such into ghci and get a result:
> foldr (+) (-1) $ takeWhile (<2000000) allPrimes'

Until next time...

Problem 11:

(To do)

Problem 12:

First, let us define a function that computes the triangle 'n':
Prelude> let triangle x = if x > 0 then x + triangle (x - 1) else 0
Second, we wish to get a list of all triangles:
Prelude> let triangles x = (triangle x):(triangles (x + 1))
Good, with this list of triangles, we need to get all the factors.  The following returns 1 if divisible:
Prelude> let factors x y = if (y `mod` x == 0) then 1 else 0
Now

Problem n:

To be written.

No comments:

Post a Comment