5

This is boring, I know, but I need a little help understanding an implementation of the Sieve of Eratosthenes. It's the solution to this Programming Praxis problem.

(define (primes n)
  (let* ((max-index (quotient (- n 3) 2))
         (v (make-vector (+ 1 max-index) #t)))
    (let loop ((i 0) (ps '(2)))
      (let ((p (+ i i 3)) (startj (+ (* 2 i i) (* 6 i) 3)))
        (cond ((>= (* p p) n)
               (let loop ((j i) (ps ps))
                  (cond ((> j max-index) (reverse ps))
                        ((vector-ref v j)
                          (loop (+ j 1) (cons (+ j j 3) ps)))
                        (else (loop (+ j 1) ps)))))
              ((vector-ref v i)
                (let loop ((j startj))
                  (if (<= j max-index)
                      (begin (vector-set! v j #f)
                             (loop (+ j p)))))
                      (loop (+ 1 i) (cons p ps)))
              (else (loop (+ 1 i) ps)))))))

The part I'm having trouble with is startj. Now, I can see that p is going to be cycling through odd numbers starting at 3, defined as (+ i i 3). But I don't understand the relationship between p and startj, which is (+ (* 2 i i) (* 6 i) 3).


Edit: I understand that the idea is to skip previously sifted numbers. The puzzle definition states that when sifting a number x, sifting should start at the square of x. So, when sifting 3, start by eliminating 9, etc.

However, what I don't understand is how the author came up with that expression for startj (algebraically speaking).

From the puzzle comments:

In general, when sifting by n, sifting starts at n-squared because all the previous multiples of n have already been sieved.

The rest of the expression has to do with the cross-reference between numbers and sieve indexes. There’s a 2 in the expression because we eliminated all the even numbers before we ever started. There’s a 3 in the expression because Scheme vectors are zero-based, and the numbers 0, 1 and 2 aren’t part of the sieve. I think the 6 is actually a combination of the 2 and the 3, but it’s been a while since I looked at the code, so I’ll leave it to you to figure out.


If anyone could help me with this, that'd be great. Thanks!

harto
  • 89,823
  • 9
  • 47
  • 61
  • 3
    Well, expressed algebraically we have `p = 2i + 3` and `startj = 2i^2 + 6i + 3`, so obviously...er...obviously...hm. That is not the clearest implementation of the Sieve that I've ever read. – David Seiler Oct 21 '09 at 00:27
  • Indeed :) A few comments would've been nice. – harto Oct 21 '09 at 00:57
  • 2
    `startj = (i+1)*p+i` It's just skipping some numbers that would have been sieved already by earlier steps, I think. – ephemient Oct 21 '09 at 01:59
  • Well, I'll just saw, not all of us think this stuff is boring. :) – BobbyShaftoe Oct 21 '09 at 04:35
  • @BobbyShaftoe yeah - I guess I thought it might be boring to ask Yet Another Question About The Sieve Of Eratosthenes, though :) – harto Oct 21 '09 at 04:46
  • `((2*i+1)^2-1) / 2 == 2*i*(i+1)` ; `((2*i+3)^2-3) / 2 == 2i^2+6i+3` . I personally prefer the first addressing scheme, `i=0 <--> p=1`, to the second, `i=0 <--> p=3`. One cell is wasted, but addressing calculations are simpler. – Will Ness Aug 25 '12 at 18:07

2 Answers2

4

I think I've figured it out, with the assistance of programmingpraxis' comments on their website.

To restate the problem, startj is defined in the listing as (+ (* 2 i i) (* 6 i) 3), i.e. 2i^2 + 6i + 3.

I didn't initially understand how this expression related to p - since 'sifting' for a number p should start at p^2, I figured that startj should be something relating to 4i^2 + 12i + 9.

However, startj is an index into the vector v, which contains odd numbers starting from 3. Therefore, the index for p^2 is actually (p^2 - 3) / 2.

Expanding the equation: (p^2 - 3) / 2 = ([4i^2 + 12i + 9] - 3) / 2 = 2i^2 + 6i + 3 - which is the value of startj.

I feel it might've been clearer to define startj as (quotient (- (* p p) 3) 2), but anyway - I think that solves it :)

harto
  • 89,823
  • 9
  • 47
  • 61
  • Ah, of course; what I was missing was the relationship between the indices into v and the numbers they represent. Nice work. +1s for all! – David Seiler Oct 21 '09 at 17:49
3

David Seiler: perhaps not the clearest, but in addition to implementing the basic Sieve it also has to implement the three optimizations described in the exercise.

Harto: that was my second exercise. I was still experimenting with my writing style.

Ephemient: correct.

See a more complete explanation in my comment at Programming Praxis.

EDIT: I've added an additional comment at Programming Praxis. And when I actually looked at the code, I was wrong about the derivation of the number 6 in the formula; sorry I mislead you.

  • Thanks for your answer. I read your comments over at the site, but I'm still confused. I don't understand why sifting doesn't start at `p^2`, i.e. `4i^2 + 12i + 9`, as opposed to the given definition of `startj`. – harto Oct 21 '09 at 03:20
  • ... so, hang on - `startj = (p^2 - 3) / 2` - this looks similar to the definition of `max-index`. Am I on the right track? – harto Oct 21 '09 at 03:22
  • I've also left an updated comment at your site. It might be helpful to update your answer here with your remarks? – harto Oct 21 '09 at 04:44