Racket Sieve of Eratosthenes using Stream

Question

When I use the following code in Racket:

#lang racket
(define (sieve x)
  (if (stream-empty? x) empty-stream
      (stream-cons (stream-first x)
            (sieve (stream-filter (λ(q)(not (zero? (modulo q (stream-first x)))))
                           (stream-rest x))))))
(define (in-primes)
  (sieve (in-naturals 2)))
(define (nth-prime n)
  (for/last ([i (in-range n)]
             [j (in-primes)]) j))

The largest number for which I can effectively compute the nth prime is 1000. Is this a reasonable implementation of the sieve of Eratosthenes, or is there something I can do to significantly speed up the above?

Answer 1

No, it's not. It's a trial division algorithm, and an extremely inefficient and suboptimal one.

Each candidate here is tested by all its preceding primes, whereas just those not greater than its square root are enough. This translates to immense worsening of complexity. I expect it runs at ~ n² at best, in n primes produced, instead of ~ n^1.45 of an optimal trial division sieve, or ~ n^1.1 of a proper sieve of Eratosthenes implementation.

The creation of filters should be postponed until a prime's square is seen among the candidates, to make it an optimal trial division algorithm.

---

You can greatly improve your code's performance with a minimal edit, following the principle of "do less, get done more": instead of calling stream-first at each step, don't. Just produce the intermediate streams in full, as they are:

(define (sieves x)
  (if (stream-empty? x) 
      empty-stream
      (stream-cons x      ; <-- here
        (sieves (stream-filter 
                   (λ (q) (not (zero? (modulo q (stream-first x)))))
                   (stream-rest x))))))

Now sieves produces a stream of streams. In each interim stream, all the numbers in the initial prefix from the first value up to its square are prime by construction. Now we can stop early, and thus drastically reduce the number of the interim streams.

To actually produce the primes, take first element from each interim stream except the last interim stream, and from the last interim stream take all elements, from the first element up to its square _{(or the desired upper limit - below that square)}. This will have roughly the same overall time complexity as the optimal trial division _{(which, at each step, takes away not just the head element from the current interim stream, but the whole prefix up to the head's square, so the next filter starts from there).}

To estimate the magnitude of n-th prime, you can use formula p_n < n * log(n*log(n)), for n > 6 (according to Wiipedia).

You can find stream-based SoE code here, though in Scheme, not Racket.

see also:

• From Turner's sieve to Bird's -- Haskell gist
• How do I define the sieve function for prime computation using higher–order functions?