Prime numbers

In mathematics, amongst the natural numbers greater than 1, a prime number (or a prime) is such that has no divisors other than itself (and 1).

Prime Number Resources

• At Wikipedia:
  • Prime Numbers
  • Sieve of Eratosthenes

• HackageDB packages:
  • arithmoi: Various basic number theoretic functions; efficient array-based sieves, Montgomery curve factorization ...
  • 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.

• Papers:
  • O'Neill, Melissa E., "The Genuine Sieve of Eratosthenes", Journal of Functional Programming, Published online by Cambridge University Press 9 October 2008 doi:10.1017/S0956796808007004.

Definition

In mathematics, amongst the natural numbers greater than 1, a prime number (or a prime) is such that has no divisors other than itself (and 1). The smallest prime is thus 2. Non-prime numbers are known as composite, i.e. those representable as product of two natural numbers greater than 1.

To find out a prime's multiples we can either a. test each new candidate number for divisibility by that prime, giving rise to a kind of trial division algorithm; or b. we can directly generate the multiples of a prime p by counting up from it in increments of p, resulting in a variant of the sieve of Eratosthenes.

The set of prime numbers is thus

   P  = {n ∈ N₂ : (∀ m ∈ N₂) ((m | n) ⇒ m = n)}

= {n ∈ N₂ : (∀ m ∈ N₂) (m×m ≤ n ⇒ ¬(m | n))}

= N₂ \ {n×m : n,m ∈ N₂}

= N₂ \ ⋃ { {n×m : m ∈ N_n} : n ∈ N₂ }

= N₂ \ ⋃ { {n×n, n×n+n, n×n+n+n, ...} : n ∈ N₂ }

= N₂ \ ⋃ { {p×p, p×p+p, p×p+p+p, ...} : p ∈ P }

  where     N_k = {n ∈ N : n ≥ k}

Thus starting with 2, for each newly found prime we can eliminate from the rest of the numbers all the multiples of this prime, giving us the next available number as next prime. This is known as sieving the natural numbers, so that in the end all the composites are eliminated and what we are left with are just primes. (This is what the last formula is describing, though seemingly impredicative, because it is self-referential. But because N₂ is well-ordered (with the order being preserved under addition), the formula is well-defined.)

In pseudocode, this can be written as

        primes = [2..] \ [[p*p, p*p+p..] | p <- primes]

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 incrementally, one by one.

Short exposition is here.

Sieve of Eratosthenes

Simplest, bounded, very inefficient formulation:

import Data.List (\\)

primesTo m = sieve [2..m]       {- (\\) is set-difference for unordered lists -}
             where 
             sieve (x:xs) = x : sieve (xs \\ [x,x+x..m])
             sieve [] = []

The (unbounded) sieve of Eratosthenes calculates primes as integers above 1 that are not multiples of primes, i.e. not composite — whereas composites are found as enumeration of multiples of each prime, generated by counting up from prime's square in constant increments equal to that prime (or twice that much, for odd primes). This is much more efficient and runs at about n^1.2 empirical orders of growth (corresponding to n log n log log n complexity, more or less, in n primes produced):

import Data.List.Ordered (minus, union, unionAll)

primes = 2 : 3 : minus [5,7..] (unionAll [[p*p, p*p+2*p..] | p <- tail primes])

{- Using `under n = takeWhile (<= n)`, with ordered increasing lists,
   `minus`, `union` and `unionAll` satisfy, for any `n` and `m`:

  under n (minus a b)         == nub . sort $ under n a \\ under n b
  under n (union a b)         == nub . sort $ under n a ++ under n b
  under n . unionAll . take m == under n . foldl union [] . take m
  under n . unionAll          == nub . sort . concat 
                                     . takeWhile (not.null) . map (under n) -}

The definition is primed with 2 and 3 as initial primes, to avoid the vicious circle.

The "big union" unionAll function could be defined as the folding of (\(x:xs) -> (x:) . union xs); or it could use a Bool array as a sorting and duplicates-removing device. The processing naturally divides up into the segments between successive squares of primes.

Stepwise development follows (the fully developed version is here).

Initial definition

First of all, working with ordered increasing lists, the sieve of Eratosthenes can be genuinely represented by

 -- genuine yet wasteful sieve of Eratosthenes
 -- primes = eratos [2.. ]               -- unbounded
primesTo m = eratos [2..m]               -- bounded, up to m
    where
    eratos []     = []
    eratos (p:xs) = p : eratos (xs `minus` [p, p+p..])
 -- eratos (p:xs) = p : eratos (xs `minus` map (p*) [1..])
 -- eulers (p:xs) = p : eulers (xs `minus` map (p*) (p:xs))    
 -- turner (p:xs) = p : turner [x | x <- xs, rem x p /= 0]

This should be regarded more like a specification, not a code. It runs at empirical orders of growth worse than quadratic in number of primes produced. But it has the core defining features of the classical formulation of S. of E. as a. being bounded, i.e. having a top limit value, and b. finding out the multiples of a prime directly, by counting up from it in constant increments, equal to that prime.

The canonical list-difference minus and duplicates-removing union functions (cf. Data.List.Ordered package) are:

 -- 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     []    = xs
union  []     ys    = ys

The name merge ought to be reserved for duplicates-preserving merging operation of the merge sort.

Analysis

So for each newly found prime p the sieve discards its multiples, enumerating them by counting up in steps of p. There are thus $O(m/p)$ multiples generated and eliminated for each prime, and $O(m\log \log (m))$ multiples in total, with duplicates, by virtues of prime harmonic series.

If each multiple is dealt with in $O(1)$ time, this will translate into $O(m\log \log (m))$ RAM machine operations (since we consider addition as an atomic operation). Indeed mutable random-access arrays allow for that. But lists in Haskell are sequential-access, and complexity of minus(a,b) for lists is $O(|a\cup b|)$ instead of $O(|b|)$ of the direct access destructive array update. The lower the complexity of each minus step, the better the overall complexity.

So on k-th step the argument list (p:xs) that the eratos function gets, starts with the (k+1)-th prime, and consists of all the numbers ≤ m coprime with all the primes ≤ p. According to the M. O'Neill's article (p.10) there are $\mathrm{\Phi }(m,p)\in \mathrm{\Theta }(m/\log p)$ such numbers.

It looks like $\sum _{i=1}^{n}1/log(p_i)=O(n/\log n)$ for our intents and purposes. Since the number of primes below m is $n=\pi (m)=O(m/\log (m))$ by the prime number theorem (where $\pi (m)$ is a prime counting function), there will be n multiples-removing steps in the algorithm; it means total complexity of at least $O(mn/\log (n))=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))$.

From Squares

But we can start each elimination step at a prime's square, as its smaller multiples will have been already produced and discarded on previous steps, as multiples of smaller primes. This means we can stop early now, when the prime's square reaches the top value m, and thus cut the total number of steps to around $n=\pi (m^{0.5})=\mathrm{\Theta }(2m^{0.5}/\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}/(\log n{)}^{0.5})$ in n primes produced, a dramatic speedup:

primesToQ m = eratos [2..m] 
    where
    eratos []     = []
    eratos (p:xs) = p : eratos (xs `minus` [p*p, p*p+p..m])
 -- eratos (p:xs) = p : eratos (xs `minus` map (p*) [p..div m p])
 -- eulers (p:xs) = p : eulers (xs `minus` map (p*) (under (div m p) (p:xs)))
 -- turner (p:xs) = p : turner [x | x<-xs, x<p*p || rem x p /= 0]

Its empirical complexity is around $O(n^{1.45})$. This simple optimization works here because this formulation is bounded (by an upper limit). To start late on a bounded sequence is to stop early (starting past end makes an empty sequence – see warning below ^₁), thus preventing the creation of all the superfluous multiples streams which start above the upper bound anyway (note that Turner's sieve is unaffected by this). This is acceptably slow now, striking a good balance between clarity, succinctness and efficiency.

^₁Warning: this is predicated on a subtle point of minus xs [] = xs definition being used, as it indeed should be. If the definition minus (x:xs) [] = x:minus xs [] is used, the problem is back and the complexity is bad again.

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 - because it explicitly stops as soon as possible:

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

(here we also flatly ignore all evens above 2 a priori.) It is now clear that it can't be made unbounded just by abolishing the upper bound m, because the guard can not be simply omitted without changing the complexity back for the worst.

Accumulating Array

So while minus(a,b) takes $O(|b|)$ operations for random-access imperative arrays and about $O(|a|)$ operations here for ordered increasing lists of numbers, using Haskell's immutable array for a one could expect the array update time to be nevertheless closer to $O(|b|)$ if destructive update were used implicitly by compiler for an array being passed along as an accumulating parameter:

{-# OPTIONS_GHC -O2 #-}
import Data.Array.Unboxed

primesToA m = sieve 3 (array (3,m) [(i,odd i) | i<-[3..m]]
                        :: UArray Int Bool)
    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

Indeed for unboxed arrays, with the type signature added explicitly (suggested by Daniel Fischer), the above code runs pretty fast, with empirical complexity of about $O(n^{1.15..1.45})$ in n primes produced (for producing from few hundred thousands to few millions primes, memory usage also slowly growing). But it relies on specific compiler optimizations, and indeed if we remove the explicit type signature, the code above turns very slow.

How can we write code that we'd be certain about? 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 still.

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 prematurely. Fixing this with explicit synchronization won't change complexity but will speed it up some more:

primesPE1 = 2 : sieve [3..] primesPE1 
    where 
    sieve xs (p:pt) | q <- p*p , (h,t) <- span (< q) xs =
                   h ++ sieve (t `minus` [q, q+p..]) pt
                -- h ++ turner [x | x<-t, rem x p>0] pt

Inlining and fusing span and (++) we get:

primesPE = 2 : ops
    where 
    ops = sieve [3,5..] 9 ops                             -- odd primes
    sieve (x:xs) q ps@ ~(p:pt)
       | x < q     = x : sieve xs q ps
       | otherwise =     sieve (xs `minus` [q, q+2*p..]) (head pt^2) pt

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 multiples streams, which were the reason for the original formulation's extreme inefficiency.

Segmented

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

primesSE = 2 : ops
    where
    ops = sieve 3 9 ops []                                -- odd primes
    sieve x q ~(p:pt) fs = 
        foldr (flip minus) [x,x+2..q-2]                   -- chain of subtractions
                           [[y+s, y+2*s..q] | (s,y) <- fs]     -- OR,
        -- [x,x+2..q-2] `minus` foldl union []            -- subtraction of merged
        --                    [[y+s, y+2*s..q] | (s,y) <- fs]            --  lists
        ++ sieve (q+2) (head pt^2) pt
                ((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.

The advantage to working in spans explicitly is that this code is easily amendable to using arrays for the composites marking and removal on each finite span; and memory usage can be kept in check by using fixed sized segments.

Segmented Tree-merging

Rearranging the chain of subtractions into a subtraction of merged streams (see below) and using tree-like folding structure, further speeds up the code and significantly improves its asymptotic time behavior (down to about $O(n^{1.28}empirically)$, space is leaking though):

primesSTE = 2 : ops
    where                  
    ops = sieve 3 9 ops []                                -- odd primes
    sieve x q ~(p:pt) fs = 
        ([x,x+2..q-2] `minus` joinST [[y+s, y+2*s..q] | (s,y) <- fs])
        ++ sieve (q+2) (head pt^2) pt
                   ((++ [(2*p,q)]) [(s,q-rem (q-y) s) | (s,y) <- fs])
    
joinST (xs:t) = (union xs . joinST . pairs) t
    where
    pairs (xs:ys:t) = union xs ys : pairs t
    pairs t         = t
joinST []     = []

Segmented merging via an array

The removal of composites is easy with arrays. Starting points can be calculated directly:

import Data.List (inits, tails)
import Data.Array.Unboxed
 
primesSAE = 2 : sieve 2 4 (tail primesSAE) (inits primesSAE)
         -- (2:) . (sieve 2 4 . tail <*> inits) $ primesSAE
  where
  sieve r q ps (fs:ft) = [n | (n,True) <- assocs (
         accumArray (\ _ _ -> False) True (r+1,q-1)
                    [(m,()) | p <- fs, let s = p * div (r+p) p,
                              m <- [s,s+p..q-1]] :: UArray Int Bool )]
      ++ sieve q (head ps^2) (tail ps) ft

The pattern of iterated calls to tail is captured by a higher-order function tails, which explicitly generates the stream of tails of a stream, making for a bit more readable (even if possibly a bit less efficient) code:

psSAGE = 2 : [n | (r:q:_, fs) <- (zip . tails . (2:) . map (^2) <*> inits) psSAGE,
                  (n,True)    <- assocs (
                                   accumArray (\_ _ -> False) True (r+1, q-1) 
                        [(m,()) | p <- fs, let s = (r+p)`div`p*p, 
                                  m <- [s,s+p..q-1]] :: UArray Int Bool )]

Linear merging

But segmentation doesn't add anything substantially, and each multiples stream starts at its prime's square anyway. What does the Postponed code do, operationally? With each prime's square passed by, there emerges 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 more and more 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? This idea is due to Richard Bird, seen in his code presented in M. O'Neill's article, equivalent to:

primesB = 2 : minus [3..] (foldr (\p r-> p*p : union [p*p+p, p*p+2*p..] r) 
                                 [] primesB)

or,

primesLME1 = 2 : prs
    where
    prs = 3 : minus [5,7..] (joinL [[p*p, p*p+2*p..] | p <- prs])
    
joinL ((x:xs):t) = x : union xs (joinL t)

Here, xs stays near the top, and more frequently odds-producing streams of multiples of smaller primes are above those of the bigger primes, that produce less frequently their multiples which have to pass through more union 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 indefinitely-defined structure, defining joinL (or foldr's combining function) to produce part of its result before accessing the rest of its input (thus making it productive).

Melissa O'Neill introduced double primes feed to prevent unneeded memoization (a memory leak). We can even do multistage. Here's the code, faster still and with radically reduced memory consumption, with empirical orders of growth of around ~ $n^{1.40}$ (initially better, yet worsening for bigger ranges):

primesLME = 2 : _Y ((3:) . minus [5,7..] . joinL . map (\p-> [p*p, p*p+2*p..]))

_Y   :: (t -> t) -> t
_Y g = g (_Y g)                -- multistage, non-sharing,  g (g (g (g ...)))
    -- g (let x = g x in x)    -- two g stages, sharing

_Y is a non-sharing fixpoint combinator, here arranging for a recursive "telescoping" multistage primes production (a tower of producers).

This allows the primesLME stream to be discarded immediately as it is being consumed by its consumer. For prs from primesLME1 definition above it is impossible, as each produced element of prs is needed later as input to the same prs corecursive stream definition. So the prs stream feeds in a loop into itself and is thus retained in memory, being consumed by self much slower than it is produced. With multistage production, each stage feeds into its consumer above it at the square of its current element which can be immediately discarded after it's been consumed. (3:) jump-starts the whole thing.

Tree merging

Moreover, it can be changed into a tree structure. This idea is due to Dave Bayer and Heinrich Apfelmus:

primesTME = 2 : _Y ((3:) . gaps 5 . joinT . map (\p-> [p*p, p*p+2*p..]))

-- joinL ((x:xs):t) = x : union xs (joinL t) 
joinT ((x:xs):t) = x : union xs (joinT (pairs t))    -- set union, ~=
  where  pairs (xs:ys:t) = union xs ys : pairs t     --    nub.sort.concat

gaps k s@(x:xs) | k < x = k:gaps (k+2) s       -- ~= [k,k+2..]\\s, 
                | True  =   gaps (k+2) xs      --   when null(s\\[k,k+2..])

This code is pretty 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).

As an aside, joinT is equivalent to infinite tree-like folding foldi (\(x:xs) ys-> x:union xs ys) []:

tree-like folding

Data.List.Ordered.foldt of the data-ordlist package builds the same structure, but in a lazier fashion, consuming its input at the slowest pace possible. Here this sophistication is not needed (evidently).

Tree merging with Wheel

Wheel factorization optimization can be further applied, and another tree structure can be used which is better adjusted for the primes multiples production (effecting about 5%-10% at the top of a total 2.5x speedup w.r.t. the above tree merging on odds only, for first few million primes):

primesTMWE = [2,3,5,7] ++ _Y ((11:) . tail  . gapsW 11 wheel 
                                    . joinT . hitsW 11 wheel)
    
gapsW k (d:w) s@(c:cs) | k < c     = k : gapsW (k+d) w s    -- set difference
                       | otherwise =     gapsW (k+d) w cs   --   k==c
hitsW k (d:w) s@(p:ps) | k < p     =     hitsW (k+d) w s    -- intersection
                       | otherwise = scanl (\c d->c+p*d) (p*p) (d:w) 
                                       : hitsW (k+d) w ps   --   k==p 

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
  -- cycle $ zipWith (-) =<< tail $ [i | i <- [11..221], gcd i 210 == 1]

The hitsW function is there to find the starting point for rolling the wheel for each prime, but this can be found directly:

primesW = [2,3,5,7] ++ _Y ( (11:) . tail . gapsW 11 wheel . joinT .
                            map (\p->  
                              map (p*) . dropWhile (< p) $
                                scanl (+) (p - rem (p-11) 210) wheel) )

Seems to run about 1.4x faster, too.

Above Limit - Offset Sieve

Another task is to produce primes above a given value:

{-# OPTIONS_GHC -O2 -fno-cse #-}
primesFromTMWE primes m = dropWhile (< m) [2,3,5,7,11] 
                          ++ gapsW a wh2 (compositesFrom a)
    where
    (a,wh2) = rollFrom (snapUp (max 3 m) 3 2)
    (h,p2:t) = span (< z) $ drop 4 primes           -- p < z => p*p<=a
    z = ceiling $ sqrt $ fromIntegral a + 1         -- p2>=z => p2*p2>a
    compositesFrom a = joinT (joinST [multsOf p  a    | p <- h ++ [p2]]
                                   : [multsOf p (p*p) | p <- t] )

snapUp v o step = v + (mod (o-v) step)              -- full steps from o
multsOf p from = scanl (\c d->c+p*d) (p*x) wh       -- map (p*) $
    where                                           --   scanl (+) x wh
    (x,wh) = rollFrom (snapUp from p (2*p) `div` p) --   , if p < from 

wheelNums  = scanl (+) 0 wheel
rollFrom n = go wheelNums wheel 
    where 
    m = (n-11) `mod` 210  
    go (x:xs) ws@(w:ws2) | x < m = go xs ws2
                         | True  = (n+x-m, ws)      -- (x >= m)

A certain preprocessing delay makes it worthwhile when producing more than just a few primes, otherwise it degenerates into simple trial division, which is then ought to be used directly:

primesFrom m = filter isPrime [m..]

Map-based

Runs ~1.7x slower than TME version, but with the same empirical time complexity, ~$n^{1.2}$ (in n primes produced) and same very low (near constant) memory consumption:

import Data.List                  -- based on http://stackoverflow.com/a/1140100
import qualified Data.Map as M

primesMPE :: [Integer]
primesMPE = 2:mkPrimes 3 M.empty prs 9   -- postponed addition of primes into map;
    where                                -- decoupled primes loop feed 
    prs = 3:mkPrimes 5 M.empty prs 9
    mkPrimes n m ps@ ~(p:pt) q = case (M.null m, M.findMin m) of
        (False, (n2, skips)) | n == n2 ->
             mkPrimes (n+2) (addSkips n (M.deleteMin m) skips) ps q
        _ -> if n<q
             then    n : mkPrimes (n+2)  m                  ps q
             else        mkPrimes (n+2) (addSkip n m (2*p)) pt (head pt^2)

    addSkip n m s = M.alter (Just . maybe [s] (s:)) (n+s) m
    addSkips = foldl' . addSkip

Turner's sieve - Trial division

David Turner's (SASL Language Manual, 1983) formulation replaces non-standard minus in the sieve of Eratosthenes by stock list comprehension with rem filtering, turning it into a trial division algorithm, for clarity and simplicity:

  -- unbounded sieve, premature filters
primesT = sieve [2..]
          where
          sieve (p:xs) = p : sieve [x | x <- xs, rem x p > 0]

  -- map head
  --    $ iterate (\(p:xs) -> [x | x <- xs, rem x p > 0]) [2..]

This creates many superfluous implicit filters, because they are created prematurely. 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 here is huge.

Guarded Filters

But this really ought to be changed into bounded and guarded variant, again achieving the "miraculous" complexity improvement from above quadratic to about $O(n^{1.45})$ empirically (in n primes produced):

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

  -- (\(a,b:_) -> map head a ++ b) . span ((< m).(^2).head) $
  --     iterate (\(p:xs) -> [x | x <- xs, rem x p > 0]) [2..m]

Postponed Filters

Or it can remain unbounded, just filters creation must be postponed until the right moment:

primesPT1 = 2 : sieve primesPT1 [3..] 
    where 
    sieve (p:pt) xs = let (h,t) = span (< p*p) xs 
                      in h ++ sieve pt [x | x <- t, rem x p > 0]

  -- fix $ concatMap (fst . fst)
  --     . iterate (\((_,xs), p:pt) -> let (h,t) = span (< p*p) xs in
  --                              ((h, [x | x <- t, rem x p > 0]), pt)) 
  --     . (,) ([2],[3..])

It can be re-written with span and (++) inlined and fused into the sieve:

primesPT = 2 : oddprimes
    where 
    oddprimes = sieve [3,5..] 9 oddprimes
    sieve (x:xs) q ps@ ~(p:pt)
        | x < q = x : sieve xs q ps
        | True  =     sieve [x | x <- xs, rem x p /= 0] (head pt^2) pt

creating here as well the linear filtering nested structure at run-time, (...(([3,5..] >>= filterBy [3]) >>= filterBy [5])...), but unlike the non-postponed code each filter being created at its proper moment, not sooner than the prime's square is seen.

filterBy ds n = [n | noDivs n ds]    -- `ds` assumed to be non-decreasing

noDivs n ds = foldr (\d r -> d*d > n || (rem n d > 0 && r)) True ds

Optimal trial division

The above is algorithmically equivalent to the traditional formulation of trial division,

ps = 2 : [i | i <- [3..],  
              and [rem i p > 0 | p <- takeWhile (\p -> p^2 <= i) ps]]

or,

-- primes = filter (`noDivs`[2..]) [2..]
--        = 2 : filter (`noDivs`[3,5..]) [3,5..]
primesTD  = 2 : 3 : filter (`noDivs` tail primesTD) [5,7..]

isPrime n = n > 1 && noDivs n primesTD

except that this code is rechecking for each candidate number which primes to use, whereas for every candidate number in each segment between the successive squares of primes these will just be the same prefix of the primes list being built.

Trial division is used as a simple primality test and prime factorization algorithm.

Segmented Generate and Test

Next we turn the list of filters into one filter by an explicit list, each one in a progression of prefixes of the primes list. This seems to eliminate most recalculations, explicitly filtering composites out from batches of odds between the consecutive squares of primes.

import Data.List

primesST = 2 : ops                 
    where
    ops = sieve 3 9 ops (inits ops)               -- odd primes
        -- (sieve 3 9 <*> inits) ops              -- [],[3],[3,5],...
    sieve x q ~(_:pt) (fs:ft) =
        filter ((`all` fs) . ((> 0).) . rem) [x,x+2..q-2]
        ++ sieve (q+2) (head pt^2) pt ft

This can also be coded as, arguably more readable,

primesSGT = 2 : ops
    where
    ops = 3 : [n | (r:q:_, px) <- 
                       (zip . tails . (3:) . map (^2) <*> inits) ops,
                   n <- [r+2,r+4..q-2], all ((> 0).rem n) px]

                -- n <- foldl (>>=) [r+2,r+4..q-2]           -- chain of filters
                --                   [filterBy [p] | p <- px]]

                -- n <- [r+2,r+4..q-2] >>= filterBy px]      -- filter by a list

Generate and Test Above Limit

The following will start the segmented Turner sieve at the right place, using any primes list it's supplied with (e.g. TMWE etc.) or itself, as shown, demand computing it just up to the square root of any prime it'll produce:

primesFromST m | m <= 2 = 2 : primesFromST 3
primesFromST m | m > 2 =
         sieve (m`div`2*2+1) (head ps^2) (tail ps) (inits ps)
    where 
    (h,ps) = span (<= (floor.sqrt $ fromIntegral m+1)) ops
    sieve x q ps (fs:ft) =
        filter ((`all` (h ++ fs)) . ((> 0).) . rem) [x,x+2..q-2]
        ++ sieve (q+2) (head ps^2) (tail ps) ft
    ops = 3 : primesFromST 5                      -- odd primes

-- ~> take 3 $ primesFromST 100000001234
-- [100000001237,100000001239,100000001249]

This is usually faster than testing candidate numbers for divisibility one by one which has to re-fetch anew the needed prime factors to test by, for each candidate. Faster is the offset sieve of Eratosthenes on odds, and yet faster the one w/ wheel optimization, on this page.

Conclusions

All these variants being variations of trial division, finding out primes by direct divisibility testing of every candidate number by sequential primes below its square root (instead of just by its factors, which is what direct generation of multiples is doing, essentially), are thus principally of worse complexity than that of Sieve of Eratosthenes.

The initial code is just a one-liner that ought to have been regarded as executable specification in the first place. It can easily be improved quite significantly with a simple use of bounded, guarded formulation to limit the number of filters it creates, or by postponement of filter creation.

Euler's Sieve

Unbounded 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 = 2 : eulers [3,5..]
  where
    eulers (p:xs) = p : eulers (xs `minus` map (p*) (p:xs))
 -- eratos (p:xs) = p : eratos (xs `minus` [p*p, p*p+2*p..])

This code is extremely inefficient, running above $O(n^2)$ empirical complexity (and worsening rapidly), and should be regarded a specification only. Its memory usage is very high, with empirical space complexity just below $O(n^2)$, in n primes produced.

In the stream-based sieve of Eratosthenes we are able to skip along the input stream xs directly to the prime's square, consuming the whole prefix at once, thus achieving the results equivalent to the postponement technique, because the generation of the prime's multiples is independent of the rest of the stream.

But here in the Euler's sieve it is dependent on all xs and we're unable in principle to skip along it to the prime's square - because all xs are needed for each prime's multiples generation. Thus efficient unbounded stream-based implementation seems to be impossible in principle, under the simple scheme of producing the multiples by multiplication.

Wheeled list representation

The situation can be somewhat 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 ([n],1) 
     [n,n+2..]         === fromElt  (n,2)    
                       === fromSpan ([n,n+2],4)     -}

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

{-                   === concat $ iterate (map (+ i)) xs
   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) -}

spanSpecs = 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 spanSpecs 
     ([2],1)
     ([3],2)
     ([5,7],6)
     ([7,11,13,17,19,23,29,31],30)  -}

Generating a list from a span specification is like rolling a wheel as its pattern gets repeated over and over again. For each span 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..1.8})$ complexity, for n primes produced, and also suffers from a severe space leak problem (IOW its memory usage is also very high).

Using Immutable Arrays

Generating Segments of Primes

The sieve of Eratosthenes' removal of multiples on each segment of odds can be done by actually marking them in an array, instead of manipulating ordered lists, and can be further sped up more than twice by working with odds only:

import Data.Array.Unboxed
 
primesSA :: [Int]
primesSA = 2 : oddprimes ()
  where
    oddprimes = (3 :) . sieve 3 [] . oddprimes
    sieve x fs (p:ps) = [i*2 + x | (i,True) <- assocs a]
                        ++ sieve (p*p) ((p,0) :
                             [(s, rem (y-q) s) | (s,y) <- fs]) ps
      where
      q = (p*p-x)`div`2
      a :: UArray Int Bool
      a = accumArray (\ b c -> False) True (1,q-1)
                     [(i,()) | (s,y) <- fs, i <- [y+s, y+s+s..q]]

Runs significantly faster than TMWE and with better empirical complexity, of about $O(n^{1.10..1.05})$ in producing first few millions of primes, with constant memory footprint.

Calculating Primes Upto a Given Value

Equivalent to Accumulating Array above, running somewhat faster (compiled by GHC with optimizations turned on):

{-# OPTIONS_GHC -O2 #-}
import Data.Array.Unboxed

primesToNA 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 a2 else f (head x) a2
      where q = p*p
            a2  :: UArray Int Bool
            a2 = a // [(i,False) | i <- [q, q+2*p..n]]
            x  = [i | i <- [p+2,p+4..n], a2 ! i]

Calculating Primes in a Given Range

primesFromToA 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   -- first odd in the segment
    r  = floor . sqrt $ fromIntegral b + 1
    ar = accumArray (\_ _ -> False) True (o,b)  -- initially all True,
          [(i,()) | p <- [3,5..r]
                    , let q  = p*p      -- flip every multiple of an odd 
                          s  = 2*p                         -- to False
                          (n,x) = quotRem (o - q) s 
                          q2 = if  o <= q  then q
                               else  q + (n + signum x)*s
                    , i <- [q2,q2+s..b] ]

Although sieving by odds instead of by 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 packed) sieve array for the whole sequence to process.

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

sieveUA :: Int -> UArray Int Bool
sieveUA top = runSTUArray $ do
    let m = (top-1) `div` 2
        r = floor . sqrt $ fromIntegral top + 1
    sieve <- newArray (1,m) True      -- :: ST s (STUArray s Int Bool)
    forM_ [1..r `div` 2] $ \i -> do
      isPrime <- readArray sieve i
      when isPrime $ do               -- ((2*i+1)^2-1)`div`2 == 2*i*(i+1)
        forM_ [2*i*(i+1), 2*i*(i+2)+1..m] $ \j -> do
          writeArray sieve j False
    return sieve
 
primesToUA :: Int -> [Int]
primesToUA top = 2 : [i*2+1 | (i,True) <- assocs $ sieveUA top]

Its empirical time complexity is improving with n (number of primes produced) from above $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 quasilinearithmic $O(n\log (n)\log (\log n))$ time 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

See Implicit Heap.

Prime Wheels

See Prime Wheels.

Using IntSet for a traditional sieve

See Using IntSet for a traditional sieve.

Testing Primality, and Integer Factorization

See Testing primality:

• Primality Test and Integer Factorization
• Miller-Rabin Primality Test

One-liners

See primes one-liners.

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 older versions of that file for your projects: it's broken and works for primes only by a lucky chance. The revised version of the file fixes the bug, as pointed out by Eugene Kirpichov on February 24, 2009 on the haskell-cafe mailing list, and fixed by Bertram Felgenhauer.

• test entries for (some of) the above code variants.

• Speed/memory comparison table for sieve of Eratosthenes variants.

• Empirical orders of growth on Wikipedia.