## Problem Description

Link to Project Euler problem 3

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.