Segmentation fault & heap errors when using Remove to generate primes

Question

I am trying to write a sieve to generate prime numbers for a Project Euler problem.

The code looks like this:

(defn sieve [n]
  (reduce (fn [memo f]
              (remove #(and (= 0 (rem % f))
                            (not= % f))
                      memo))
          (range 2 (+ 1 n))
          (range 2 (+ 1 n))))

Until 500000 it runs very fast, under 1 sec, from 600000 and up it starts seg faulting and crashing with memory errors.

I imagine it has something to do with remove and lazyness, i searched a bit, tried using (doall (remove ...)) instead of (remove) but it becomes incredibly slow...

I am a bit at a loss with this, thanks in advance for any help...

It should be stack over flow because of the lots of lazy calls of `remove` being executed when the resulting seq is asked to materialize — Ankur, Mar 20 '12 at 16:55
Which JVM implementation are you using? (use "java -version" to find out) — raek, Mar 20 '12 at 17:42

score 3 · Answer 1 · answered Mar 20 '12 at 17:39

Segfault? That sounds scary! I assume you mean a stack overflow error.

When I run the function I start to get stack overflow errors at about 1000. So why the stack overflows? It has to do with laziness. The fix is to wrap the call to remove in a doall.

Reduce will iterate through each element in the sequence given as the third argument and keep a state along the way. This state is initialized as the rage of integers from 2 to n+1. The state is updated at each iteration using remove. However, since remove is lazy it won't actually do anything. Remove will return an object that can generate a sequence on demand, based on the sequence it was given. I will try to explain this thorugh an example:

(reduce (fn [xs x] (remove #{x} xs)) coll (range 4))

The above expression will return a sequence of the elements of coll, but with numbers 0 to 3 filtered out. To explain what happens at run time I will invent a new notation: I will write the runtime value of (remove #{x} xs) as «I.O.U. a seq like xs but with x removed». The value of the reduce expression is then:

«I.O.U. a seq like
  «I.O.U. a seq like
    «I.O.U. a seq like
      «I.O.U. a seq like 'coll' but with 0 removed»
    but with 1 removed»
  but with 2 removed»
but with 3 removed»

Each call to remove adds a wrapping "I.O.U.". This causes problems when you finally try to access the first value of the resulting seq (the outermost "I.O.U." value). When a lazy-seq object is traversed, it is first checked whether its value has been calculated. If the value is already done, the value is simply returned. If it's not, then it's calculated, stored, and returned.

The problem arises when one lazy-seq ("I.O.U." value) needs to force another lazy-seq to be able to perform its work, because one stack frame is needed per lazy-seq realization. (A stack frame is needed in order to remember where to return to when the second lazy-seq is done.) If you have 1000 nested "I.O.U." values, you need 1000 stack frames to realize them (assuming all of them were unrealized initiallly).

The solution is to use the doall function from the Clojure standard library. It takes a seq and causes it to be fully realized. If you wrap the call to remove in a doall, the reduce state will always contain a fully realized seq between each iteration, and no cascade of "I.O.U." values will build up. Instead of saving all computations for later, you do them incrementally.

Adrian Mouat · Answer 2 · 2012-03-20T16:55:05.927

I'm not surprised - you're working with a huge list there.

I think you might need a fundamentally different algorithm. For example; to test if a number n is prime, you need to check if it can be divided by any prime <= the square root of n. Using this knowledge we can start building a list of primes by testing numbers in sequential order, adding each new prime to the list.

This is a slow algorithm but can be sped up by using a "wheel" that skips obvious non-primes (e.g. numbers that are divisible 2 or 3).

This is all of the top of my head, so apologies for any inaccuracies.

score 2 · Answer 3 · edited May 23 '17 at 12:04

In principal, prime calculation is better suited to associateive structures like sets than iterative structures like lists. A few other StackOverflowers have contributed useful clojure answers:

The answer to this quiestion,
this question: why is this prime sieve implementation slower,
this question on fast prime generation in lisp

in general I like to break this sort of problem into individual steps so here is my answer using sets:

(defn next-not-in-sieve [sieve start]
  (first (drop-while #(sieve %) (range 1 (inc start)))))

(defn add-multiples [sieve max n] (into sieve (range n (inc max) n)))

(defn primes [max sieve answer]
  (let [next (next-not-in-sieve sieve max)]
    (if (nil? next) answer
        (recur max (add-multiples sieve max next) (conj answer next)))))

it runs sufficiently quickly for basic use, the point here is to learn Clojure, not find primes quickly of course :)

user> (time (def foo (sieve 60000)))
"Elapsed time: 63167.82765 msecs"

user> (time (def foo (primes 60000 #{1} [])))
"Elapsed time: 33272.157235 msecs"

And what would be the point of learning Clojure if we didn't make it into a lazy seq of primes:

(defn primes [max sieve]
  (if-let [next (next-not-in-sieve sieve max)]
    (lazy-seq (cons next (primes max (add-multiples sieve max next))))))

and check the time:

(time (def foo (doall (primes 60000 #{1}))))
"Elapsed time: 33714.880871 msecs"

and of course for thous who want some back ground check out the wikipedia page on prime sieves

Segmentation fault & heap errors when using Remove to generate primes

3 Answers3