haskelloptimizationtime-complexityprimessieve-algorithm

Efficient (O(n^2)) Sieve of Sundaram in Haskell


There are many answers on SO explaining how to implement the Sieve of Sundaram in Haskell but they're all... really inefficient?

All solutions that I have seen work like this:

  1. Figure out the numbers <= n to exclude
  2. Filter these from [1..n]
  3. Modify the remaining numbers * 2 + 1

Here for example is my implementation that finds all primes between 1 and 2n+2:

sieveSundaram :: Integer -> [Integer]
sieveSundaram n = map (\x -> 2 * x + 1) $ filter (flip notElem toRemove) [1..n]
  where toRemove = [i + j + 2*i*j | i <- [1..n], j <- [i..n], i + j + 2*i*j <= n]

The problem I have with this, is that filter has to traverse the entire toRemove list for every element of [1..n] and thus this has complexity O(n^3) whereas a straightforward iterative implementation has complexity O(n^2). How can I achieve that in Haskell?


Solution

  • As per the comments, base should not be considered a complete standard library for Haskell. There are several key packages that every Haskell developer knows and uses and would consider part of Haskell's de facto standard library.

    By "straightforward iterative implementation", I assume you mean marking and sweeping an array of flags? It would be usual to use a Vector or Array for this. (Both would be considered "standard".) An O(n^2) Vector solution looks like the following. Though it internally uses a mutable vector, the bulk update operator (//) hides this fact, so you can write it in a typical Haskell immutable and stateless style:

    import qualified Data.Vector as V
    
    primesV :: Int -> [Int]
    primesV n = V.toList                           -- the primes!
      . V.map (\x -> (x+1)*2+1)                    -- apply transformation
      . V.findIndices id                           -- get remaining indices
      . (V.// [(k - 1, False) | k <- removals n])  -- scratch removals
      $ V.replicate n True                         -- everyone's allowed
    
    removals n = [i + j + 2*i*j | i <- [1..n], j <- [i..n], i + j + 2*i*j <= n]
    

    Another possibility that's a little more straightforward is IntSet which is basically a set of integers with O(1) insertion/deletion and O(n) ordered traversal. (This is like the HashSet suggested in the comments, but specialized to integers.) This is in the containers packages, another "standard" package that's actually bundled with the GHC source, even though it's distinct from base. It gives an O(n^2) solution that looks like:

    import qualified Data.IntSet as I
    
    primesI :: Int -> [Int]
    primesI n = I.toAscList               -- the primes!
      . I.map (\x -> x*2+1)               -- apply transformation
      $ I.fromList [1..n]                 -- integers 1..n ...
        I.\\ I.fromList (removals n)      -- ... except removals
    

    Note that another important performance improvement is to use a better removals definition that avoids filtering all n^2 combinations. I believe the following definition produces the same list of removals:

    removals :: Int -> [Int]
    removals n = [i + j + 2*i*j | j <- [1..(n-1) `div` 3], i <- [1..(n-j) `div` (1+2*j)]]
    

    and does so in what I believe is O(n log(n)). If you use it with either primesV or primesI above, it's the bottleneck, so the resulting overall algorithm should be O(n log(n)), I think.