Sieve of Eratosthenes

The problem we are trying to solve here is to find out all primes up to a given natural number. Eratosthenes came up with an algorithm, famously known as the Sieve of Eratosthenes. The gist of the algorithm is to pick up a prime and then cross off all its multiples.

We start with writing the numbers from 2 to \(n\). Let's call the first uncrossed number \(p\). Declare \(p\) as a prime and cross all of its multiples starting at \(p^2\). Starting at \(p^2\) is a slight optimization because the smaller numbers must have been crossed by the primes we discovered earlier. Repeat this procedure until you reach beyond \(\sqrt{n}\). Declare the numbers that remain uncrossed as primes.

The home page of haskell.org has the following definition to illustrate the ease with which Haskell lets us define infinite structures like the sequence of prime numbers.

primes = filterPrime [2..]
  where filterPrime (p:xs) =
          p : filterPrime [x | x <- xs, x `mod` p /= 0]

But the above algorithm isn't The Sieve of Eratosthenes even though it is given that name very often. The paper talks about the real sieve pretty well.

The above algorithm and its real counterpart can be improved by providing a better list of numbers to the function \(filterPrime\). For example, it would be better if we avoid checking even numbers and all multiples of 3. This can be achieved by generating the list to be passed using what is known as a wheel.

Imagine yourself standing on the number line. You are initially on the number 1. Now, I am supposed to give you the number of steps you must take so that you end up on the next number that isn't a multiple of 2 or 3; or in general not a multiple of any of the numbers in a set. What would be the steps that you would take if that set of numbers happens to be empty? It's trivially a never-ending sequence of 1's. You just walk the whole number line stepping upon each number.

Now, let's see what the sequence is when we have one element in the set, say \(p\). Since you are standing on 1, the first multiple is yet to be encountered and is at \(p\). All numbers between 1 and \(p\) are numbers you must step upon. So the sequence must contain \((p - 2)\) values equal to 1 followed by a single 2. The 2 in the sequence is for the point in time when I ask you to jump over \(p\). After this, the same numbers just repeat themselves indefinitely.

genWheel :: Int -> Wheel
genWheel x = let ys = replicate (x - 2) 1 ++ [2] ++ ys in Wheel ys

Eventually, we would want our wheels to be generated for larger sets of numbers because we would like to avoid as many multiples of already known primes as possible. So, we want to combine two wheels, say genWheel 2 and genWheel 3. We definitely expect the combine operation to be associative and the wheel generated for the empty set to be its identity element. Yes, it is a Monoid.

data Wheel = Wheel { getWheel :: [Int] }
           deriving Show

instance Monoid Wheel where
  mempty = Wheel (repeat 1)
  mappend (Wheel xs) (Wheel ys) = Wheel (go xs ys)
    where go as@(x:xs) bs@(y:ys) | x == y = x : go xs ys
                                 | x > y = go as (y + head ys : tail ys)
                                 | otherwise = go (x + head xs : tail xs) bs

Our combine operation is defined above as mappend. The idea is to merge consecutive elements of the wheel when we encounter a mismatch between the two wheels to be combined. It makes sense because the sum of consecutive elements on a wheel constitutes a valid step size. So, in essence we are trying to morph the smaller wheel to look more like the larger wheel whenever we see a local difference in the two. So, the identity wheel, i.e. the one generated from the empty set, morphs trivially into any other wheel by always growing to make its elements equal to the elements in the other wheel.

Now, all we need is a function that would help us generate the sequence of numbers that we must check for primality given a wheel. So, let's spin the wheel.

spin :: Wheel -> [Int]
spin (Wheel (x:xs)) = spin' xs (x + 1)
  where spin' (y:ys) n = n : spin' ys (n + y)

We can either start with 1 or \(x + 1\) where \(x\) is where the wheel starts. Starting with 1 makes sure that we look at numbers smaller than \(x\). In the code snippet above, we have generated numbers starting with \(x + 1\) and hence skipped numbers [2..x].

Here's an example showing it in action.

λ> take 10 $ spin (genWheel 2)
[3,5,7,9,11,13,15,17,19,21]
λ> take 10 $ spin (genWheel 2 `mappend` genWheel 3)
[5,7,11,13,17,19,23,25,29,31]
λ> take 10 $ spin (genWheel 2 `mappend` genWheel 3)
[5,7,11,13,17,19,23,25,29,31]