Pritchard’s Wheel Sieve

January 6, 2012

We begin at the end with the main function that drives the others; setup is done in the let*, then the loop processing updates the set S until k is greater than m:

(define (pritchard n)
  (let* ((ps (eratosthenes (isqrt n)))
         (m (length ps))
         (k (compute-k n ps)))
    (let loop ((k (+ k 1)) (s (initial-s n k ps)))
      (if (< m k) (append ps (cdr s))
        (loop (+ k 1) (next n s))))))

The next function derives a new version of S from its predecessor:

(define (next n s)
  (let* ((p (cadr s)) (pg (make-pg p)))
    (let loop ((ss s) (ps (list p)))
      (let ((p (+ (car ps)
                  (vector-ref pg (/ (- (cadr ss) (car ss) 2) 2)))))
        (if (< n p) (list-minus s (reverse ps))
          (loop (cdr ss) (cons p ps)))))))

Next calls make-pg to build a lookup-table of extended gaps, converting a multiplication to a lookup, and list-minus subtracts one list from another:

(define (make-pg p)
  (let ((pg (make-vector (- p 1) 0)))
    (do ((i 0 (+ i 1))) ((= i (- p 1)) pg)
      (vector-set! pg i (* 2 p (+ i 1))))))

(define (list-minus xs ys)
  (let loop ((xs xs) (ys ys) (zs (list)))
    (cond ((null? xs) (reverse zs))
          ((null? ys) (append (reverse zs) xs))
          ((= (car xs) (car ys))
            (loop (cdr xs) (cdr ys) zs))
          (else (loop (cdr xs) ys (cons (car xs) zs))))))

Compute-k and initial-s are used in the setup phase:

(define (compute-k n ps)
  (let loop ((k 0) (m 2) (ps (cdr ps)))
    (if (< (/ n (log n)) m) k
      (loop (+ k 1) (* m (car ps)) (cdr ps)))))

(define (initial-s n k ps)
  (let ((s (make-vector (+ n 1) #t)))
    (vector-set! s 0 #f)
    (let loop ((k k) (ps ps))
      (if (zero? k) (enlist s)
        (do ((i (car ps) (+ i (car ps))))
            ((< n i) (loop (- k 1) (cdr ps)))
          (vector-set! s i #f))))))

The set S is initially represented as a bit vector locally in initial-s then translated to a list by the enlist function:

(define (enlist s)
  (let loop ((i (- (vector-length s) 1)) (ss (list)))
    (cond ((negative? i) ss)
          ((vector-ref s i)
            (loop (- i 1) (cons i ss)))
          (else (loop (- i 1) ss)))))

We also use the sieve of Eratosthenes from a previous exercise to initialize the small list of primes, and isqrt from the Standard Prelude. Here’s the program in action:

> (pritchard 100)
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)
> (length (pritchard 1000))
168

You can run the program at http://programmingpraxis.codepad.org/22VkAlSX. You will find that Pritchard’s sieve is quite slow compared to Eratosthenes’ sieve. Though that is true in general, our implementation is also rather inefficient due to its use of the set data structure. A better implementation uses arrays, but if you decide to do that, be careful; a naive translation from sets to array indexed front-to-back will fail because some of the values of the S set are changed (stricken from the set) before the sweep through the array catches them, causing errors. There are various ways to solve that problem (Pritchard swept the array back-to-front, other programmers came up with other methods), but all add to the problem that Pritchard’s sieve is slower in practice than Eratosthenes’, even with a better asymptotic complexity. Other versions of Pritchard’s sieve improve on this version, but are still slower in practice than Eratosthenes’. In fact, all of the advanced sieves that we have studied are slower than the original algorithm from the ancient Greek mathematician.

About these ads

Pages: 1 2

3 Responses to “Pritchard’s Wheel Sieve”

  1. Giorgio Valoti said

    If I understand correctly, the wheel for 2, 3, 5, 7 is [10 2 4 2 4 6 2 6 4 2 4 6 6 2 6 4 2 6 4 6 8 4 2 4 2 4 8 6 4 6 2 4 6 2 6 6 4 2 4 6 2 6 4 2 4 2 10 2]. However, in this paper http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf, the wheel starts from 2: [2 4 2 4 6 2 6 4 2 4 6 6 2 6 4 2 6 4 6 8 4 2 4 2 4 8 6 4 6 2 4 6 2 6 6 4 2 4 6 2 6 4 2 4 2 10 2 10].

    Why the difference?

  2. Sanjay Giri said

    There is an error here The second time through the loop, p = 7, the gaps are 7−1=6, 11−7=4, 13−11=6, 17−13=4, and 19−17=2, and the numbers that are stricken are 7, 7+6×7=49, 49+4×7=77, and 77+2×7=91, leaving S4 = {1 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97}.

    i.e 13-11 is 2

  3. nondescript said

    Thanks for the algorithm!

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 630 other followers

%d bloggers like this: