Can we improve upon this primes sieve code from SICP

Will Ness 0
Viewed 35

A recent Q&A entry showcased the following primes generating code from SICP, using lazy streams:

(define (sieve stream)
  (cons-stream
   (stream-car stream)
   (sieve (stream-filter
            (lambda (x)
              (not (divisible? x (stream-car stream))))
            (stream-cdr stream)))))

(define primes (sieve (integers-starting-from 2)))

An answer there showed primes to be equivalent, among other possibilities, to the following:

  (cons-stream 2
   (cons-stream 3
    (cons-stream 5
     (cons-stream 7
       (sieve 
         (stream-filter (lambda (x) (not (divisible? x 7)))
          (stream-filter (lambda (x) (not (divisible? x 5)))
           (stream-filter (lambda (x) (not (divisible? x 3)))
            (stream-filter (lambda (x) (not (divisible? x 2)))
             (integers-starting-from 9))))))))))

It seems there are too many filter streams here -- for instance 7 was produced by filtering the input numbers by 2, 3 and 5, whereas it only really had to be tested by 2 alone -- only the numbers above 9 need really be test divided by 3, let alone by 5 etc.

This problem becomes more and more pronounced as we go along producing this stream of primes. Overall, producing first n primes takes O(n^2) with this code.

Can we do better?

sicp lazy-sequences sieve primes scheme
1 Answers

Indeed, we need to only start filtering out multiples of a prime after its square is encountered in the input. For that, we shall use the primes and their squares. And we'll use the same code to produce those primes we need to produce our primes:

(define (pprimes) 
  (cons-stream 2
    (psieve (pprimes)      ;; here           ;; v---- and here too
            (stream-map (lambda (x) (* x x)) (pprimes))
            (integers-starting-from 3))))

(define (psieve primes psquares numbers)
  (if (< (stream-car numbers) (stream-car psquares))
    (cons-stream 
       (stream-car numbers)
       (psieve primes            ;; use the same prime and its 
               psquares          ;;  square for the next number
               (stream-cdr numbers)))
    (psieve 
       (stream-cdr primes)       ;; advance the primes stream
       (stream-cdr psquares)     ;;    and their squares, in sync
       (stream-filter            ;; and start filtering by this prime now
               (lambda (x)
                   (not (divisible? x (stream-car primes))))
               (stream-cdr numbers)))))

Now this leads to

  (pprimes)
=
  ....
=
  (cons-stream 2
   (cons-stream 3
    (cons-stream 5
     (cons-stream 7
      (cons-stream 11
       (cons-stream 13
        (cons-stream 17
         (cons-stream 19
           (psieve (cons-stream 5 ... )
                   (cons-stream 25 ... )
              (stream-filter (lambda (x) (not (divisible? x 3)))
               (stream-filter (lambda (x) (not (divisible? x 2)))
                  (integers-starting-from 20))))))))))))
=
  ....

which, undoubtedly, is much better. No number below 25 will be tested by 5, etc.

This is still trial division, and runs in about n^1.5. True sieve of Eratosthenes should run at n log n log log n which is empirically usually close to n^1.1..1.2 or thereabouts. But this n^1.5 is a great improvement over the quadratic algorithm too, and in practice will run much faster than it in absolute terms as well.

Will Ness0
Related