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.
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
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
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:
sievefunction filters the input sequence and find the primes incrementally
wheelstream and the
spinfunction 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
primesstream 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