# Euler #7 : Benchmarked Primes

by Nicolas Wu

Posted on 12 August 2010

This week’s Project Euler question is:

Find the 10001st prime.

We already described a simple algorithm for finding primes in a previous post, so rather than repeat ourselves, in this article we’ll discuss benchmarking using Criterion to find the fastest prime number algorithm that doesn’t require too much magic. I’ll be taking implementations found on the Haskell wiki.

## Imported modules

First we’ll need to import the Criterion modules, which provides us with our benchmarking suite:

import System.Environment (getArgs, withArgs)
import Criterion (bgroup, bench, nf)
import Progression.Main (defaultMain)
import Data.List.Ordered (minus, union)

I’ll actually be using Progression in conjuction with Criterion, which just makes collecting the results of several benchmarks a little easier.

## Prime Algorithms

The prime number generator we used previously was Turner’s sieve, defined as follows:

turner :: [Int]
turner =  sieve [2 .. ]
where
sieve (p:xs) = p : sieve [x | x <- xs, x mod p /= 0]

The Haskell wiki documents a whole range of other algorithms that can be used to generate primes. Here are the definitions that I pulled from the wiki:

postSieve :: [Int]
postSieve = 2 : 3 : sieve (tail postSieve) [5,7..]
where
sieve (p:ps) xs = h ++ sieve ps [x | x <- t, x rem p /= 0]
where (h,~(_:t)) = span (< p*p) xs

trialOdds :: [Int]
trialOdds = 2 : 3 : filter isPrime [5,7..]
where
isPrime n   = all (notDivs n)
$takeWhile (\p-> p*p <= n) (tail trialOdds) notDivs n p = n mod p /= 0 nestedFilters :: [Int] nestedFilters = 2 : 3 : sieve [] (tail nestedFilters) 5 where notDivsBy d n = n mod d /= 0 sieve ds (p:ps) x = foldr (filter . notDivsBy) [x,x+2..p*p-2] ds ++ sieve (p:ds) ps (p*p+2) spansPrimes :: [Int] spansPrimes = 2 : 3 : sieve 0 (tail spansPrimes) 5 where sieve k (p:ps) x = [n | n <- [x,x+2..p*p-2], and [nremp/=0 | p <- fs]] ++ sieve (k+1) ps (p*p+2) where fs = take k (tail spansPrimes) bird :: [Int] bird = 2 : primes' where primes' =  ++ [5,7..] minus foldr union' [] mults mults = map (\p -> let q=p*p in (q,tail [q,q+2*p..]))$ primes'
union' (q,qs) xs = q : union qs xs

wheel :: [Int]
wheel = 2:3:primes'
where
1:p:candidates = [6*k+r | k <- [0..], r <- [1,5]]
primes'        = p : filter isPrime candidates
isPrime n      = all (not . divides n)
$takeWhile (\p -> p*p <= n) primes' divides n p = n mod p == 0 I won’t go into the details of explaining these different algorithms, since I want us to focus on how we might benchmark these implementations. ## Benchmarking In order to compare these different algorithms, we construct a program that takes as its argument the name of the function that should be used to produce prime numbers. Once the user has provided this input, the benchmark is executed using Criterion to produce the first 101, 1001, and 10001 primes. main = do args <- getArgs let !primes = case head args of "turner" -> turner "postSieve" -> postSieve "trialOdds" -> trialOdds "nestedFilters" -> nestedFilters "spansPrimes" -> spansPrimes "bird" -> bird "wheel" -> wheel _ -> error "prime function unkown!" withArgs (("-n" ++ (head args)) : tail args)$ do
defaultMain . bgroup "Primes" $[ bench "101"$ nf (\n -> primes !! n) 101
, bench "1001" $nf (\n -> primes !! n) 1001 , bench "10001"$ nf (\n -> primes !! n) 10001
]

We then run this code with each prime function name as an argument individually, and the Progression library puts the results together. Here’s a bar chart generated from the data:

These results have been normalised against the turner function, and show the results of how long it took for the various algorithms to find the 10001th, 1001th and 100th primes.
euler7 = spansPrimes 10001