0

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?

Will Ness
  • 70,110
  • 9
  • 98
  • 181

1 Answers1

0

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.

Will Ness
  • 70,110
  • 9
  • 98
  • 181
  • Have you read [The Genuine Sieve of Eratosthenes](https://www.cambridge.org/core/journals/journal-of-functional-programming/article/genuine-sieve-of-eratosthenes/FD3E90871269020CA6C64C25AB8A4FBD)? (In Haskell, but the ideas should be "portable".) – molbdnilo Sep 14 '22 at 06:00
  • Yes. And indeed they are. Except the "sieve of Bird" is a much better basis for development. Btw there's lots of related stuff on Rosettacode's sieve of Eratosthenes page's Scheme (and Haskell) section. As well as haskellwiki's prime numbers page. @molbdnilo I'm focusing on SICP styled approach here. – Will Ness Sep 14 '22 at 07:35
  • I find that article very confusing. It misses this most immediate imorovement of "postponement" of filtering which is applicable there as well, and even claims that "using squares" is no improvement at all (which *is* true, but for random access arrays, not lists). Two valuable parts in the article, for me, were the complexity derivations and Richard Bird's code in the comments. @molbdnilo – Will Ness Sep 14 '22 at 07:56
  • @molbdnilo IOW, the "true reason" that Turner's sieve is so slow (and the original SICP one as well) is not that it uses lists instead of priority queues, but that it does not postpone the filtering properly. – Will Ness Sep 14 '22 at 08:03
  • Excellent points. To be honest, I haven't read that paper in a good while and forgot both the details and that there was something odd about it. Memory like a sieve... – molbdnilo Sep 14 '22 at 08:40