What's the ideal implementation for the Sieve of Eratosthenes between Lists, Arrays, and Mutable Arrays?

haskell, performance, primes, sieve-of-eratosthenes

Solution

This

primes = sieve [2..] where
         sieve (p:x) = p : sieve [ n | n <- x, n `mod` p > 0 ]

is not a sieve. It's very inefficient trial division. Don't use that!

I'm curious about how you got your times, there is no way that the Turner "sieve" could produce the primes not exceeding 2,000,000 in mere seconds. Letting it find the primes to 200,000 took

MUT     time    6.38s  (  6.39s elapsed)
GC      time    9.19s  (  9.20s elapsed)
EXIT    time    0.00s  (  0.00s elapsed)
Total   time   15.57s  ( 15.59s elapsed)

on my box (64-bit Linux, ghc-7.6.1, compiled with -O2). The complexity of that algorithm is `O(N² / log² N)`, almost quadratic. Letting it proceed to 2,000,000 would take about twenty minutes.

Your times for the array versions are suspicious too, though in the other direction. Did you measure interpreted code?

Sieving to 2,000,000, compiled with optimisations, the mutable array code took 0.35 seconds to run, and the immutable array code 0.12 seconds.

Now, that still has the mutable array about three times slower than the immutable array.

But, it's an unfair comparison. For the immutable array, you used `Int`s, and for the mutable array `Integer`s. Changing the mutable array code to use `Int`s - as it should, since under the hood, arrays are `Int`-indexed, so using `Integer` is an unnecessary performance sacrifice that buys nothing - made the mutable array code run in 0.15 seconds. Close to the mutable array code, but not quite there. However, you let the mutable array do more work, since in the immutable array code you only eliminate odd multiples of the odd primes, but in the mutable array code, you mark all multiples of all primes. Changing the mutable array code to treat 2 specially, and only eliminate odd multiples of odd primes brings that down to 0.12 seconds too.

But, you're using range-checked array indexing, which is slow, and, since the validity of the indices is checked in the code itself, unnecessary. Changing that to using `unsafeRead` and `unsafeWrite` brings down the time for the immutable array to 0.09 seconds.

Then you have the problem that using

forM_ [x, y .. z]

uses boxed `Int`s (fortunately, GHC eliminates the list). Writing a loop yourself, so that only unboxed `Int#`s are used, the time goes down to 0.02 seconds.

{-# LANGUAGE MonoLocalBinds #-}
import Control.Monad (forM_, when)
import Control.Monad.ST
import Data.Array.ST
import Data.Array.Unboxed
import Data.Array.Base

primeSieve :: Int -> UArray Int Bool
primeSieve top = runSTUArray $ do
    a <- newArray (0,top) True
    unsafeWrite a 0 False
    unsafeWrite a 1 False
    let r = ceiling . sqrt $ fromIntegral top
        mark step idx
            | top < idx = return ()
            | otherwise = do
                unsafeWrite a idx False
                mark step (idx+step)
        sift p
            | r < p     = return a
            | otherwise = do
                prim <- unsafeRead a p
                when prim $ mark (2*p) (p*p)
                sift (p+2)
    mark 2 4
    sift 3

-- Return primes from sieve as list:
primesTo :: Int -> [Int]
primesTo top = [p | (p,True) <- assocs $ primeSieve top]

main :: IO ()
main = print .last $ primesTo 2000000

So, wrapping up, for a Sieve of Eratosthenes, you should use an array - not surprising, its efficiency depends on being able to step from one multiple to the next in short constant time.

You get very simple and straightforward code with immutable arrays, and that code performs decently for not too high limits (it gets relatively worse for higher limits, since the arrays are still copied and garbage-collected, but that's not too bad).

When you need better performance, you need mutable arrays. Writing efficient mutable array code is not entirely trivial, one has to know how the compiler translates the various constructs to choose the right one, and some would consider such code unidiomatic. But you can also use a library (disclaimer: I wrote it) that provides a fairly efficient implementation rather than writing it yourself.

Problem

Original source

Related problems

• Parallel Algorithms for Generating Prime Numbers (possibly using Hadoop's map reduce)