Prime numbers

In mathematics, a prime number (or a prime) is a natural number which has exactly two distinct natural number divisors: 1 and itself. The smallest prime is thus 2.

Prime Number Resources

• At Wikipedia:
  • Prime Numbers
  • Sieve of Eratosthenes

• HackageDB packages:
  • Numbers: An assortment of number theoretic functions.
  • NumberSieves: Number Theoretic Sieves: primes, factorization, and Euler's Totient.
  • primes: Efficient, purely functional generation of prime numbers.

Any natural number is representable as a product of powers of its prime factors, and so a prime has no prime divisors other than itself. That means that starting with 2, for each newly found prime we can eliminate from the rest of the numbers all such numbers that have this prime as their divisor, giving us the next available number as next prime. This is known as sieving the natural numbers, so that in the end what we are left with are just primes.

To eliminate a prime's multiples from the result we can either a. plainly test each new candidate number for divisibility by that prime with a direct test, giving rise to a kind of "trial division" algorithm; or b. we can find out multiples of a prime p by counting up from it by p numbers at a time, resulting in a variant of a "genuine sieve" as it was reportedly originally conceived by Eratosthenes in ancient Greece.

Having a direct-access mutable arrays indeed enables easy marking off of these multiples on pre-allocated array as it is usually done in imperative languages; but to get an efficient list-based code we have to be smart about combining those streams of multiples of each prime - which gives us also the memory efficiency in generating the results one by one.

Sieve of Eratosthenes

Sieve of Eratosthenes is genuinely represented by

 -- genuine yet wasteful Eratosthenes
 primesTo m = 2 : sieve [3,5..m]  where
    sieve []     = []
    sieve (p:xs) = p : sieve (xs `minus` [p,p+2*p..m])

This should be regarded more like a specification, not a code. It is extremely slow, running at empirical time complexities worse than quadratic in number of primes produced. But it has the core defining features of S. of E. as a) being bounded, i.e. having a top limit value, and b) finding out the multiples of a prime by counting up from it. Yes, this is exactly how Eratosthenes defined it (Nicomachus, Introduction to Arithmetic, I, pp. 13, 31).

The canonical list-difference minus and duplicates-removing list-union union functions dealing with ordered increasing lists - infinite as well as finite - are simple enough to define (using compare has an effect of comparing the values only once, unlike when using (<) etc):

 -- ordered lists, difference and union
 minus (x:xs) (y:ys) = case (compare x y) of 
           LT -> x : minus  xs  (y:ys)
           EQ ->     minus  xs     ys 
           GT ->     minus (x:xs)  ys
 minus  xs     _     = xs
 union (x:xs) (y:ys) = case (compare x y) of 
           LT -> x : union  xs  (y:ys)
           EQ -> x : union  xs     ys 
           GT -> y : union (x:xs)  ys
 union  xs     ys    = xs ++ ys

(the name merge ought to be reserved for duplicates-preserving merging as used by mergesort - that's why we use union here, following Leon P.Smith's Data.List.Ordered package).

Analysis

So what does it do, this sieve code? For each found prime it removes its odd multiples from further consideration. It finds them by counting up in steps of 2p. There are thus $O(m/p)$ multiples for each prime, and $O(m\log(\log m))$ multiples total, with duplicates, by virtue of prime harmonic series, where $\sum _{p_{i}<m}{1/p_{i}}=O(\log(\log m))$.

If each multiple is dealt with in $O(1)$ time, this will translate into $O(m\log(\log m))$ RAM-mashine operations (since we consider addition as an atomic operation). Indeed, mutable random-access arrays allow for that. But lists in Haskell don't, and complexity of minus(a,b) is about $O(|union(a,b)|)$, i.e. here it is $O(|a|)$ which is $O(m/\log p)$ according to the notice about the Φ-function in Melissa O'Neill's article.

It looks like $\sum _{i=1}^{k}{1/log(p_{i})}=O(k/\log k)$. Since the number of primes below m is $k=\pi (m)=O(m/\log m)$ by the prime number theorem (where $\pi (m)$ is a prime counting function), there will be k multiples-removing steps in the algorithm; it means total complexity of $O(mk/\log k)=O(m^{2}/(\log m)^{2})$, or $O(n^{2})$ in n primes produced - much much worse than the optimal $O(n(\log n)\log(\log n))$.

Squared

But we can start each step at the prime's square, as its smaller multiples will be already processed on previous steps. This means we can stop early, when the prime's square reaches the top value m, and thus cut the total number of steps to around $k=\pi ({\sqrt {m}})=O({\sqrt {m}}/\log m)$. This does not in fact change the complexity of random-access code, but for lists it makes it $O(m^{1.5}/(\log m)^{2})$, or $O(n^{1.5}/{\sqrt {\log n}})$ in n primes produced, showing an enormous speedup in practice:

 primesToQ m = 2 : sieve [3,5..m]  where
    sieve []     = []
    sieve (p:xs) = p : sieve (xs `minus` [p*p,p*p+2*p..m])

Its empirical complexity is about $O(n^{1.45})$. This simple optimization of starting from a prime's square has dramatic effect here because our formulation is bounded, in accordance with the original algorithm. This has the desired effect of stopping early, turning the unacceptably slow initial specification into a code, yet-far-from-optimal and slow, but acceptably so.

Guarded

This ought to be explicated (improving on clarity though not on time complexity) as in the following, for which it is indeed a minor optimization whether to start from p or p*p - but only after we've went to the necessary trouble of explicitly stopping early on:

 primesToG m = 2 : sieve [3,5..m]  where
    sieve (p:xs) 
       | p*p > m = p : xs
       | True    = p : sieve (xs `minus` [p*p,p*p+2*p..m])

Did Eratosthenes himself achieve the optimal complexity? It rather seems doubtful, as he probably counted boxes in the table by 1 to go from one number to the next, as in 3,5,7,9,11,13,15, ... for he had no access even to Hindu numerals, using Greek alphabet for writing down numbers instead. Was he performing a genuine sieve of Eratosthenes then? We'll leave that as an open question.

Accumulating Array

So while minus(a,b) takes $O(|b|)$ operations for random-access imperative arrays and about $O(|a|)$ operations for lists here, using Haskell's immutable array for a one could expect the array update time to be indeed $O(|b|)$ but it looks like it's not so:

 primesToA m = sieve 3 (array (3,m) [(i,True) | i<-[3..m]] //
                                    [(i,False) | i<-[4,6..m]])
  where
   sieve p a 
    | p*p > m   = 2 : [i | (i,True) <- assocs a]
    | a!p       = sieve (p+2) $ a // [(i,False) | i <- [p*p,p*p+2*p..m]]
    | otherwise = sieve (p+2) a

It's much slower than the above, though it should be running at optimal complexity on implementations which are able to properly use the destructive update here, for an array being passed along as an accumulating parameter, it seems.

How this implementation deficiency is to be overcome? One way is to use explicitly mutable monadic arrays (see below), but we can also think about it a little bit more on the functional side of things.

Postponed

Going back to guarded Eratosthenes, first we notice that though it works with minimal number of prime multiples streams, it still starts working with each a bit prematurely. Fixing this with explicit synchronization won't change complexity but will speed it up some more:

 primesPE = 2 : primes'
  where 
    primes' = sieve [3,5..] primes' 9
    sieve (p:xs) ps@ ~(_:t) q
       | p < q  = p : sieve xs ps q
       | True   =     sieve (xs `minus` [q,q+2*p..]) t (head t^2)

Since the removal of a prime's multiples here starts at the right moment, and not just from the right place, the code could now finally be made unbounded. Because no multiples-removal process is started prematurely, there are no extraneous mutiples streams, which were the reason for the extreme wastefulness and thus inefficiency of the original formulation.

Segmented

With work done segment-wise between the successive primes squares it becomes

  primes = 2: 3: sieve (tail primes) 3 []
   where
    sieve (p:ps) x fs = let q=p*p in
        foldr (flip minus) [x+2,x+4..q-2] 
                           [[y+s,y+2*s..q] | (s,y) <- fs]
        ++ sieve ps q 
               ((2*p,q):[(s,q-rem (q-y) s) | (s,y) <- fs])

This "marks" the odd composites in a given range by generating them - just as a person performing the original sieve of Eratosthenes would do, counting one by one the multiples of the relevant primes. These composites are independently generated so some will be generated multiple times. Rearranging the chain of subtractions into a subtraction of merged streams and using tree-like folding further speeds up the code and improves its asympthotic behaviour.

The advantage in working with spans explicitly is that this code is easily amendable to using arrays for the composites marking and removal on each finite span.

Linear merging

But what does the Postponed code do, operationally? For each prime's square passed over, it builds up a nested linear left-deepening structure, (...((xs-a)-b)-...), where xs is the original odds-producing [3,5..] list, so that each odd it produces must go through `minus` nodes on its way up - and those odd numbers that eventually emerge on top are prime. Thinking a bit about it, wouldn't another, right-deepening structure, (xs-(a+(b+...))), be better?

Here, xs would stay near the top, and more frequently odds-producing streams of multiples of smaller primes would be above those of the bigger primes, which produce less frequently their candidates which have to pass through more `minus` nodes on their way up. Plus, no explicit synchronization is necessary anymore because the produced multiples of a prime start at its square anyway - just some care has to be taken to avoid a runaway access to the infinitely-defined structure (specifically, if each (+) operation passes along unconditionally its left child's head value before polling the right child's head value, then we are safe).

Here's the code, faster yet but still with about same time complexity of $O(n^{1.4})$, which is essentially equivalent to the Richard Bird's code from Melissa O'Neill's article:

 {-# OPTIONS_GHC -O2 -fno-cse #-}
 primesLME = 2 : ([3,5..] `minus` fold [[p*p,p*p+2*p..] | p <- primes']) 
  where
   primes' = 3 : ([5,7..] `minus` fold [[p*p,p*p+2*p..] | p <- primes'])   
   fold ((x:xs):t) = x : union xs (fold t)

The double primes feed is introduced here to prevent unneeded memoization and thus save memory, as per Melissa O'Neill's code, and is dependent on no expression sharing being performed by a compiler - space complexity here is probably $O(\pi ({\sqrt {m}}))=O(2{\sqrt {n/\log n}})$ but in practice it's negligible compared to run-time system memory itself and thus the total memory is very nearly $O(1)$.

Tree merging

Moreover, it can be changed into a tree structure:

 {-# OPTIONS_GHC -O2 -fno-cse #-}
 primesTME = 2 : ([3,5..] `minus` foldt [[p*p,p*p+2*p..] | p <- primes']) 
  where
   primes' = 3 : ([5,7..] `minus` foldt [[p*p,p*p+2*p..] | p <- primes'])   
   foldt ((x:xs):t)    = x : union xs (foldt (pairs t))
   pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t

It is very fast, running at speeds and empirical complexities comparable with the code from Melissa O'Neill's article (about $O(n^{1.2})$ in number of primes n produced).

Tree merging with Wheel

Furthermore, wheel factorization optimization can be applied to it, and another tree structure can be used, better adjusted for the primes multiples production (effecting an additional 5%~10% speedup):

 {-# OPTIONS_GHC -O2 -fno-cse #-}
 primesTMWE = 2:3:5:7: gaps 11 wheel (fold3t $ roll 11 wheel primes')
  where
   primes' = 11: gaps 13 (tail wheel) (fold3t $ roll 11 wheel primes')
   fold3t ((x:xs): ~(ys:zs:t)) = x : union xs (union ys zs)    
                                      `union` fold3t (pairs t)  
   pairs ((x:xs):ys:t)         = (x : union xs ys) : pairs t
   wheel = 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:wheel
   gaps k ws@(w:t) cs@(c:u) | k==c  = gaps (k+w) t u 
                            | True  = k : gaps (k+w) t cs  
   roll k ws@(w:t) ps@(p:u) | k==p  = scanl (\c d->c+p*d) (p*p) ws 
                                        : roll (k+w) t u 
                            | True  = roll (k+w) t ps

The original idea of the right-deepening unbounded linear merging of multiples is due to Richard Bird (presented in Melissa O'Neill's article), and the original idea of using a tree-like folding structure is due to Dave Bayer, on haskell-cafe mailing list.

Turner's sieve - Trial division

David Turner's original 1975 formulation (SASL Language Manual, 1975) replaces non-standard `minus` with stock list comprehension with rem filtering in the sieve of Eratosthenes, turning it into a kind of trial division algorithm:

 -- unbounded sieve, premature filters
 primesT = 2 : sieve [3,5..]  where
    sieve (p:xs) = p : sieve [x | x<-xs, rem x p /= 0]
                          -- filter ((/=0).(`rem`p)) xs

This creates an immense number of superfluous implicit filters in extremely premature fashion. To be admitted as prime, each number will be tested for divisibility here by all its preceding primes, while just those not greater than its square root would suffice. To find e.g. the 1001st prime (7927), 1000 filters are used, when in fact just the first 24 are needed (up to 89's filter only). Operational overhead is enormous here.

Guarded

When we force ourselves away from the Quest for a Mythical One-Liner it really ought to be written at least as bounded and guarded variety (if not abandoned right away in favor of algorithmically superior sieve of Eratosthenes), yet again achieving the miraculous complexity improvement from above quadratic to about $O(n^{1.45})$ empirically (in n primes produced):

 primesToGT m = 2 : sieve [3,5..m]  where
    sieve (p:xs) 
       | p*p > m = p : xs
       | True    = p : sieve [x | x<-xs, rem x p /= 0]
                          -- filter ((/=0).(`rem`p)) xs

Postponed

or, better, as a postponed, unbounded variety:

 primesPT = 2 : primes'
  where 
    primes' = sieve [3,5..] primes' 9
    sieve (p:xs) ps@ ~(_:t) q
       | p < q   = p : sieve xs ps q
       | True    =     sieve [x | x<-xs, rem x p /= 0] t (head t^2)
                          -- filter ((/=0).(`rem`p)) xs

Segmented

A certain further speedup (though not in complexity which stays the same) can be achieved by explicating the nested filters structure created here as well at run-time, turning a list of filters into a filter by an explicit list of primes (which is just a prefix of the primes list itself that is being built):

 primesST = 2 : primes' 
  where
    primes' = 3 : sieve primes' 3 0
    sieve (p:ps) x k = 
       let pfx = take k primes'
       in  [n | n <- [x+2,x+4..p*p-2], and [rem n f /= 0 | f <- pfx]]
           -- filter (\n-> all ((/=0).(n`rem`)) pfx) [x+2,x+4..p*p-2]
           ++ sieve ps (p*p) (k+1)

This seems to be closest to optimal trial division, with most recalculations eliminated, explicitly filtering composites out from batches of odds between the consequent squares of primes. All these variants of course being variations of trial division, finding out primes by direct divisibility testing of every odd number by all primes below its square root (instead of just its factors, which is what eliminating the multiples is doing, essentially), are thus principally of worse complexities than the sieve of Eratosthenes; but one far worse than another, yet fixed easily from a wasteful monstrosity to almost acceptible performance (at least for the first few hunderd thousand primes, when compiled) just by following the proper definition of the sieve as being bounded, simply with the guarded formulation.

Another task is to produce primes above a given number (not knowing their ordinal numbers this time), which will demand compute just the primes below its square root:

  primesFrom m   = sieve ps (m`div`2*2+1) (length h)
    where 
      (h,(_:ps)) = span (<= (floor.sqrt$fromIntegral m+1)) primesST

Conclusion

So it really pays off to analyse the code properly instead of just labeling it "naive". BTW were this divisibility testing somehow turned into an O(1) operation, e.g. by some kind of massive parallelization, the overall complexity would drop to O(n). It's the sequantiality of testing that's the culprit. Though of course the proper multiples-removing S. of E. is a better candidate for parallelization.

So the initial Turner code is just a one-liner that ought to have been regarded as specification only, in the first place, not a code to be executed as such. The reason it was taught that way is probably so that it could provoke this discussion among the students. To regard it as code is what's been naive all along.

Euler's Sieve

With each found prime Euler's sieve removes all its multiples in advance so that at each step the list to process is guaranteed to have no multiples of any of the preceding primes in it (consists only of numbers coprime with all the preceding primes) and thus starts with the next prime:

  primesEU = euler [2..] where
    euler (p:xs) = p : euler (xs `minus` map (p*) (p:xs))

This code is extremely inefficient, running at about O(${n^{2.4}}$) complexity, and should be regarded a specification only. It works very hard trying to avoid double hits on multiples, which are relatively few and cheap to deal with in the first place.

Wheeled list representation

The situation can be improved using a different list representation, for generating lists not from a last element and an increment, but rather a last span and an increment, which entails a set of helpful equivalences:

  {- fromElt (x,i) = x : fromElt (x + i,i)
                       === iterate (+ i) x
     [n..]             === fromElt (n,1)     -}

  fromSpan (xs,i)  = xs ++ fromSpan (map (+ i) xs,i)

  {-                   === concat $ iterate (map (+ i)) xs
     [n..]             === fromSpan ([n],1)
     fromSpan (p:xt,i) === p : fromSpan (xt ++ [p + i], i)  
     fromSpan (xs,i) `minus` fromSpan (ys,i) 
                       === fromSpan (xs `minus` ys, i)  
     map (p*) (fromSpan (xs,i)) 
                       === fromSpan (map (p*) xs, p*i)
     fromSpan (xs,i)   === forall (p > 0).
       fromSpan (concat $ take p $ iterate (map (+ i)) xs, p*i) -}

  eulerS () = iterate eulerStep ([2],1)
  eulerStep (xs@(p:_), i) = 
       ( (tail . concat . take p . iterate (map (+ i))) xs
          `minus` map (p*) xs, p*i )

  {- > mapM_ print $ take 4 $ eulerS ()
     ([2],1)
     ([3],2)
     ([5,7],6)
     ([7,11,13,17,19,23,29,31],30)  -}

This is where the notion of wheel comes from naturally, all by itself. For each wheel specification w@((p:_),_) produced by eulerStep, the numbers in (fromSpan w) up to ${p^{2}}$ are all primes too, so that

  eulerPrimesTo m = if m > 1 then go ([2],1) else []
    where
      go w@((p:_), _) 
         | m < p*p = takeWhile (<= m) (fromSpan w)
         | True    = p : go (eulerStep w)

This runs at about O(${n^{1.5}}$) complexity, for n primes produced.

Using Immutable Arrays

Generating Segments of Primes

The removal of multiples on each segment of odds can be done by actually marking them in an array instead of manipulating the ordered lists, and can be further sped up more than twice by working with odds only, represented as their offsets in segment arrays:

  primes = 2: 3: sieve (tail primes) 3 []
   where 
    sieve (p:ps) x fs = [i*2 + x | (i,e) <- assocs a, e] 
                        ++ sieve ps (p*p) fs'
     where
      q     = (p*p-x)`div`2                  
      fs'   = (p,0) : [(s, rem (y-q) s) | (s,y) <- fs]
      a     = accumArray (\ b c -> False) True (1,q-1)
                         [(i,()) | (s,y) <- fs, i <- [y+s,y+s+s..q]]

Apparently, arrays are fast. The above is the fastest code of all presented so far. When run on Ideone.com it is somewhat faster than Prime_numbers#Tree Merging With Wheel in producing first few million primes, but is very much unacceptable in its memory consumption which grows faster than O(${n}$), quickly getting into tens and hundreds of MBs.

Calculating Primes Upto a Given Value

  primesToN n = 2: [i | i <- [3,5..n], ar ! i]
   where
    ar = f 5 $ accumArray (\ a b -> False) True (3,n) 
                        [(i,()) | i <- [9,15..n]]
    f p a | q > n = a
          | True  = if null x then a' else f (head x) a'
      where q = p*p
            a'= a // [(i,False) | i <- [q,q+2*p..n]]
            x = [i | i <- [p+2,p+4..n], a' ! i]

Calculating Primes in a Given Range

  primesFromTo a b = (if a<3 then [2] else []) 
                     ++ [i | i <- [o,o+2..b], ar ! i]
   where 
    o  = max (if even a then a+1 else a) 3
    r  = floor.sqrt.fromInteger $ b+1
    ar = accumArray (\a b-> False) True (o,b) 
          [(i,()) | p <- [3,5..r]
                    , let q  = p*p 
                          s  = 2*p 
                          (n,x) = quotRem (o - q) s 
                          q' = if  o <= q  then q
                               else  q + (n + signum x)*s
                    , i <- [q',q'+s..b] ]

Although using odds instead of primes, the array generation is so fast that it is very much feasible and even preferable for quick generation of some short spans of relatively big primes.

Using Mutable Arrays

Using mutable arrays is the fastest but not the most memory efficient way to calculate prime numbers in Haskell.

Using ST Array

This method implements the Sieve of Eratosthenes, similar to how you might do it in C, modified to work on odds only. It is fast, but about linear in memory consumption, allocating one (though apparently bit-packed) array for the whole sequence produced.

import Control.Monad
import Control.Monad.ST
import Data.Array.IArray
import Data.Array.MArray
import Data.Array.ST
import Data.Array.Unboxed

primesToNA :: Int -> UArray Int Bool
primesToNA n = runSTUArray sieve  where
  sieve = do
    let m = (n-1)`div`2
    a <- newArray (1,m) True :: ST s (STUArray s Int Bool)
    let sr = floor . (sqrt::Double->Double) . fromIntegral $ n+1
    forM_ [1..sr`div`2] $ \i -> do
      let ii = 2*i*(i+1)     -- == ((2*i+1)^2-1)`div`2
      si <- readArray a i
      when si $
        forM_ [ii,ii+i+i+1..m] $ \j -> writeArray a j False
    return a
 
primesToN :: Int -> [Int]
primesToN n = 2:[i*2+1 | (i,e) <- assocs . primesToNA $n, e]

Its empirical time complexity is improving with n (number of primes produced) from O(n^1.25) through O(n^1.20) towards O(n^1.16). The reference C++ vector-based implementation exhibits this improvement in empirical time complexity too, from O(n^1.5) gradually towards O(n^1.12), where tested (which might be interpreted as evidence towards the expected $O(n*\log n*\log(\log n))$ complexity).

Bitwise prime sieve with Template Haskell

Count the number of prime below a given 'n'. Shows fast bitwise arrays, and an example of Template Haskell to defeat your enemies.

{-# OPTIONS -O2 -optc-O -XBangPatterns #-}
module Primes (nthPrime) where

import Control.Monad.ST
import Data.Array.ST
import Data.Array.Base
import System
import Control.Monad
import Data.Bits

nthPrime :: Int -> Int
nthPrime n = runST (sieve n)

sieve n = do
    a <- newArray (3,n) True :: ST s (STUArray s Int Bool)
    let cutoff = truncate (sqrt $ fromIntegral n) + 1
    go a n cutoff 3 1

go !a !m cutoff !n !c
    | n >= m    = return c
    | otherwise = do
        e <- unsafeRead a n
        if e then
          if n < cutoff then
            let loop !j
                 | j < m     = do
                        x <- unsafeRead a j
                        when x $ unsafeWrite a j False
                        loop (j+n)
                 | otherwise = go a m cutoff (n+2) (c+1)
            in loop ( if n < 46340 then n * n else n `shiftL` 1)
           else go a m cutoff (n+2) (c+1)
         else go a m cutoff (n+2) c

And place in a module:

{-# OPTIONS -fth #-}
import Primes

main = print $( let x = nthPrime 10000000 in [| x |] )

Run as:

$ ghc --make -o primes Main.hs
$ time ./primes
664579
./primes  0.00s user 0.01s system 228% cpu 0.003 total

Implicit Heap

The following is an original implicit heap implementation for the sieve of Eratosthenes, kept here for historical record. The Prime_numbers#Tree Merging with Wheel section above simplifies it, removing the People a structure altogether, and improves upon it by using a folding tree structure better adjusted for primes processing, and a wheel optimization.

See also the message threads Re: "no-coding" functional data structures via lazyness for more about how merging ordered lists amounts to creating an implicit heap and Re: Code and Perf. Data for Prime Finders for an explanation of the People a structure that makes it work.

data People a = VIP a (People a) | Crowd [a]

mergeP :: Ord a => People a -> People a -> People a
mergeP (VIP x xt) ys                    = VIP x $ mergeP xt ys
mergeP (Crowd xs) (Crowd ys)            = Crowd $ merge  xs ys
mergeP xs@(Crowd (x:xt)) ys@(VIP y yt)  = case compare x y of
    LT -> VIP x $ mergeP (Crowd xt) ys
    EQ -> VIP x $ mergeP (Crowd xt) yt
    GT -> VIP y $ mergeP xs yt

merge :: Ord a => [a] -> [a] -> [a]
merge xs@(x:xt) ys@(y:yt) = case compare x y of
    LT -> x : merge xt ys
    EQ -> x : merge xt yt
    GT -> y : merge xs yt

diff xs@(x:xt) ys@(y:yt) = case compare x y of
    LT -> x : diff xt ys
    EQ ->     diff xt yt
    GT ->     diff xs yt

foldTree :: (a -> a -> a) -> [a] -> a
foldTree f ~(x:xs) = x `f` foldTree f (pairs xs)
    where pairs ~(x: ~(y:ys)) = f x y : pairs ys

primes, nonprimes :: [Integer]
primes    = 2:3:diff [5,7..] nonprimes
nonprimes = serve . foldTree mergeP . map multiples $ tail primes
    where
    multiples p = vip [p*p,p*p+2*p..]

    vip (x:xs)       = VIP x $ Crowd xs
    serve (VIP x xs) = x:serve xs
    serve (Crowd xs) = xs

nonprimes effectively implements a heap, exploiting lazy evaluation.

Prime Wheels

The idea of only testing odd numbers can be extended further. For instance, it is a useful fact that every prime number other than 2 and 3 must be of the form $6k+1$ or $6k+5$. Thus, we only need to test these numbers:

primes :: [Integer]
primes = 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

Here, primes' is the list of primes greater than 3 and isPrime does not test for divisibility by 2 or 3 because the candidates by construction don't have these numbers as factors. We also need to exclude 1 from the candidates and mark the next one as prime to start the recursion.

Such a scheme to generate candidate numbers first that avoid a given set of primes as divisors is called a prime wheel. Imagine that you had a wheel of circumference 6 to be rolled along the number line. With spikes positioned 1 and 5 units around the circumference, rolling the wheel will prick holes exactly in those positions on the line whose numbers are not divisible by 2 and 3.

A wheel can be represented by its circumference and the spiked positions.

data Wheel = Wheel Integer [Integer]

We prick out numbers by rolling the wheel.

roll (Wheel n rs) = [n*k+r | k <- [0..], r <- rs]

The smallest wheel is the unit wheel with one spike, it will prick out every number.

w0 = Wheel 1 [1]

We can create a larger wheel by rolling a smaller wheel of circumference n along a rim of circumference p*n while excluding spike positions at multiples of p.

nextSize (Wheel n rs) p =
  Wheel (p*n) [r' | k <- [0..(p-1)], r <- rs,
                      let r' = n*k+r, r' `mod` p /= 0]

Combining both, we can make wheels that prick out numbers that avoid a given list ds of divisors.

mkWheel ds = foldl nextSize w0 ds

Now, we can generate prime numbers with a wheel that for instance avoids all multiples of 2, 3, 5 and 7.

primes :: [Integer]
primes = small ++ large
    where
    1:p:candidates = roll $ mkWheel small
    small          = [2,3,5,7]
    large          = p : filter isPrime candidates
    isPrime n      = all (not . divides n) 
                       $ takeWhile (\p -> p*p <= n) large
    divides n p    = n `mod` p == 0

It's a pretty big wheel with a circumference of 210 and allows us to calculate the first 10000 primes in convenient time.

A fixed size wheel is fine, but how about adapting the wheel size while generating prime numbers? See Euler's Sieve above, or the functional pearl titled Lazy wheel sieves and spirals of primes for more.

Using IntSet for a traditional sieve

module Sieve where
import qualified Data.IntSet as I

-- findNext - finds the next member of an IntSet.
findNext c is | I.member c is = c
              | c > I.findMax is = error "Ooops. No next number in set."
              | otherwise = findNext (c+1) is

-- mark - delete all multiples of n from n*n to the end of the set
mark n is = is I.\\ (I.fromAscList (takeWhile (<=end) (map (n*) [n..])))
                where
                    end = I.findMax is

-- primes - gives all primes up to n 
primes n = worker 2 (I.fromAscList [2..n])
                where
                    worker x is 
                     | (x*x) > n = is
                     | otherwise = worker (findNext (x+1) is) (mark x is)

(doesn't look like it runs very efficiently).

Testing Primality

See Testing primality.

External links

• http://www.cs.hmc.edu/~oneill/code/haskell-primes.zip

A collection of prime generators; the file "ONeillPrimes.hs" contains one of the fastest pure-Haskell prime generators; code by Melissa O'Neill.

WARNING: Don't use the priority queue from that file for your projects: it's broken and works for primes only by a lucky chance.

• http://ideone.com/willness/primes

test entries for (some of) the above code variants