Prime sieve results in (me going to) stack overflow

93 views Asked by At

For project euler i was looking for a way to implement the sieve of Eratosthenes because i expected it to be faster than my original implementation (and needed it to be so). I came up with this function:

sieve :: Int -> [Int] -> [Int]
sieve p (x:xs) | rem x p == 0 = sieve p xs
               | otherwise = x : (sieve x $ sieve p xs)
sieve _ [] = []

Which works, but hits the fan very fast, resulting in a stack overflow. I headed here for some advice and immediately struck a spoiler which, to me, looks quite the same, but the difference in performance is outlandish. I still want to continue using my own implementation, and would like to know whether or not my function can be changed easily to use less memory.

1

There are 1 answers

0
Will Ness On BEST ANSWER

Your code has an exponential explosion inside it:

sieve p (x:xs) | rem x p == 0 = sieve p xs
               | otherwise = x : (sieve x $ sieve p xs)
--                                          ^^^^^^^       here!
--                                ^^^^^^^                 and here!

You intended for the inner sieve call to just continue filtering by p, but since you use the same sieve function, it will also start up new filters for new primes as it encounters them - but that's entirely superfluous since the "upper" level call will also start new filter for the same prime!

sieve 2 [3..]
 = 3 : sieve 3 (sieve 2 [4..])
 = 3 : sieve 3 (5 : sieve 5 (sieve 2 [6..]))
 = 3 : 5 : sieve 5 (sieve 3 (sieve 5 (sieve 2 [6..])))
....         -- ^^^               ^^^                      -- !!!

After 7 passes to the top through the four sieves here, each will split in two creating four new sieve 7s, so there will be eight sieves in the chain!! Even worse, when 9 - the composite! - goes through the sieves, it will split the 2, the 7s and the 5, and will only be rejected by the 3. So it's actually worse than exponential in number of primes n, and close to being exponential in the value of the last prime produced (which is ~ n log(n)).

Change it to

sieve p (x:xs) | rem x p == 0 = sieve p xs
               | otherwise = x : (sieve x $ filter ((>0).(`rem` p)) xs)

and you get the code equivalent to the one you cite.

If you prefer to hand code everything, you could introduce a new argument, a Boolean flag controlling whether to start up new filters or not:

primes = 2 : sieve True 2 [3..]
sieve b p (x:xs) | rem x p == 0 = sieve b p xs
                 | b = x : (sieve b x $ sieve False p xs)
                 | otherwise = x : sieve b p xs