Can we improve upon this primes sieve code from SICP

Question

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?

Answer 1

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 these primes, which we need for producing our primes:

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

(define (psieve pr-sqrs numbers)       ;; prime+square pairs
  (if (< (stream-car numbers)          
         (cdr (stream-car pr-sqrs)))   ;; prime's square
    (cons-stream 
       (stream-car numbers)
       (psieve pr-sqrs                 ;; same prime+square pair
               (stream-cdr numbers)))  ;;   for the next number
    (psieve 
       (stream-cdr pr-sqrs)            ;; advance prime+square's stream
       (stream-filter                  ;; and start filtering 
          (let ((p  (car (stream-car pr-sqrs)))) ;; by this prime now
               (lambda (x)
                   (not (divisible? x p))))
          (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 original quadratic algorithm, too, and in practice will run much faster than it in absolute terms as well.

By the way, changing (not (divisible? x n)) to ((not-divisible-by n) x) opens our code up to the whole new line of algorithmic improvements, using additions only (and no attempted divisions) just like in the original, from ~2000 years ago.