Project Euler in Haskell #3

2012-04-12 — Article of 800 words - Available in Italiano Programming Haskell Project Euler

Problem Description

The prime factors of 13195 are 5, 7, 13 and 29.

What is the largest prime factor of the number 600851475143?

WARNING Solution ahead. Don’t read more if you want to enjoy the benefits of Project Euler and you haven’t already solved the problem.

Solution

We need a primes generator to find the prime factors of a number. Instead of using a package from Hackage, I prefer to create a generator from scratch.

I will use the sieve of Eratosthenes algorithm, as described by Melissa E. O’Neill (PDF 217KB).

The generator proposed in the article use a priority queue data structure. A simple and effective priority queue, the Pairing Heap, is described by Louis Wasserman (PDF 396KB).

Here is the code, adapted to the need of the sieve algorithm, and without using the Maybe monad for min value extraction (for simplicity’s sake).

data PairingHeap a  =  Empty
                    |  PairNode a [PairingHeap a]
                       deriving (Show)

(+++)              :: Ord a => PairingHeap a -> PairingHeap a -> PairingHeap a
Empty  +++  heap   =  heap
heap   +++  Empty  =  heap
heap1@(PairNode x1 ts1)  +++  heap2@(PairNode x2 ts2)
    | x1 <= x2     =  PairNode x1 (heap2:ts1)
    | otherwise    =  PairNode x2 (heap1:ts2)

-- Construction
empty  :: PairingHeap a
empty  =  Empty

singleton    :: a -> PairingHeap a
singleton x  =  PairNode x []

-- Insertion
insert      :: Ord a => a -> PairingHeap a -> PairingHeap a
insert x q  =  (singleton x) +++ q

-- Priority Queue
minValue                 :: PairingHeap a ->  a
minValue Empty           =  error "Empty queue!"
minValue (PairNode x _)  =  x

deleteMin                  :: Ord a => PairingHeap a -> PairingHeap a
deleteMin Empty            =  error "Empty queue!"
deleteMin (PairNode _ ts)  =  meldChildren ts
    where
      meldChildren (t1:t2:ts)  =  t1 +++ t2 +++ meldChildren ts
      meldChildren [t]         =  t
      meldChildren []          =  Empty

deleteMinAndInsert      :: Ord a => a -> PairingHeap a -> PairingHeap a
deleteMinAndInsert x q  =  insert x $ deleteMin q

This Pairing Heap implementation doesn’t support a key/value node. So I define an Iterator data type to manage the key/value concept necessary to the algorithm.

-- Iterator
data Iterator = Iterator Integer [Integer]
    deriving (Show)

instance Eq Iterator where
    (Iterator x _) == (Iterator y _) = x == y

instance Ord Iterator where
    compare (Iterator x _) (Iterator y _) = compare x y

I even create a Table type alias and some utility functions to hide the priority queue specific implementation from the algorithm.

-- Table
type Table = PairingHeap Iterator

emptyTable  :: Table
emptyTable  =  empty

insertTable         :: Integer -> [Integer] -> Table -> Table
insertTable n is t  =  insert (Iterator n is) t

tableMinValue    :: Table -> Integer
tableMinValue t  =  n
    where
      (Iterator n _) = minValue t

tableMinValueIters    :: Table -> (Integer, [Integer])
tableMinValueIters t  =  (n, is)
    where
      (Iterator n is)  =  minValue t

tableDeleteMinAndInsert         :: Integer -> [Integer] -> Table -> Table
tableDeleteMinAndInsert n is t  =  deleteMinAndInsert (Iterator n is) t

With the priority queue completed, we can define the sieve algorithm:

the sieve function filters the input sequence and find the primes incrementally
the wheel stream and the spin function allow to generate an input sequence without the composites of the primes 2,3,5,7; this optimization makes the sieve seven times faster than with a complete input sequence [2..]
the primes stream is created composing the sieve and the spinned wheel

sieve         :: [Integer] -> [Integer]
sieve []      =  []
sieve (x:xs)  =  x : sieve' xs (insertprime x xs emptyTable)
    where
      insertprime p xs table  =  insertTable (p*p) (map (*p) xs) table
      sieve' [] table  =  []
      sieve' (x:xs) table
          | nextComposite <= x  =  sieve' xs (adjust table)
          | otherwise           =  x : sieve' xs (insertprime x xs table)
          where
            nextComposite  =  tableMinValue table
            adjust table
                | n <= x     =  adjust (tableDeleteMinAndInsert n' ns table)
                | otherwise  =  table
                where
                  (n, n':ns)  =  tableMinValueIters table

wheel2357 :: [Integer]
wheel2357 = 2:4:2:4:6:2:6:4:2:4:6:6:2:6:4:2:6:4:6:8:4:2:4:2:4
            :8:6:4:6:2:4:6:2:6:6:4:2:4:6:2:6:4:2:4:2:10:2:10:wheel2357

spin           :: [Integer] -> Integer -> [Integer]
spin (x:xs) n  =  n : spin xs (n+x)

primes  :: [Integer]
primes  =  2:3:5:7: sieve (spin wheel2357 11)

I use the primes generator to find the prime factors of an integer.

primeFactors    :: Integer -> [Integer]
primeFactors n  =  factors n primes
    where
      factors 1 _                   =  []
      factors m (p:ps) | m < p*p    =  [m]
                       | r == 0     =  p : factors q (p:ps)
                       | otherwise  =  factors m ps
                       where (q,r)  =  quotRem m p

And, finally, the solution:

solution = last $ primeFactors 600851475143

You can find the Literate Haskell code on GitHub and on Bitbucket.