Sieve of Eratosthenes in Scheme (R7RS)

I've seen many implementations of Sieve of Eratosthenes in Scheme, but I thought I'd try to write one that is both space- and time-efficient:

• Space-efficient: I use R7RS bytevectors as a bitset, which encodes odd numbers starting with 3 (so bit 0 = 3, bit 1 = 5, bit 2 = 7, etc.).
• Time-efficient: Many implementations scan the table looking for the next viable prime. I instead use first-set-bit (analogous to Integer.numberOfTrailingZeros in Java) to avoid scanning each bit individually.

This code generates prime numbers up to at least n, and usually a handful more. It depends on SRFI 60 for the bitwise operations.

(define (primes-up-to-about n)
(define table
(make-bytevector (+ (arithmetic-shift n -4) 1) 255))
(define len (bytevector-length table))
(define bitlen (arithmetic-shift len 3))

(define (clear! i)
(define q (arithmetic-shift i -3))
(bytevector-u8-set! table q
(copy-bit (bitwise-and i 7) (bytevector-u8-ref table q) #f)))
(define (next-marked-from i)
(let loop ((q (arithmetic-shift i -3))
(m (arithmetic-shift -1 (bitwise-and i 7))))
(and (< q len)
(let ((v (bitwise-and (bytevector-u8-ref table q) m)))
(cond ((zero? v) (loop (+ q 1) -1))
(else (+ (arithmetic-shift q 3) (first-set-bit v))))))))
(define (index->value i)
(+ i i 3))

(let loop ((rv '(2)) (i 0))
(define next (next-marked-from i))
(if (not next)
(reverse rv)
(let ((nextval (index->value next)))
(do ((j next (+ j nextval)))
((>= j bitlen))
(clear! j))
(loop (cons nextval rv) next)))))


I'm looking for ways to make the code more compact and/or elegant without sacrificing the space or time efficiency.

In particular, if you can make a SRFI 42 :bitset generator that actually does the equivalent of my next-marked-from, I can just turn the whole function into a list-ec comprehension. What's not to like? :-)

For people who want to test this program in Racket (since there aren't very many R7RS implementations currently around), add the following lines before the function definition:

(require srfi/60)
(define bytevector-length bytes-length)
(define bytevector-u8-ref bytes-ref)
(define bytevector-u8-set! bytes-set!)
(define make-bytevector make-bytes)

• Hard coded numbers do not make the code easier to understand. Comments paying homage [or at least lip service] to HtDP might also make the code more accessible. Even a traditional Lisp style documentation string would help. There's a case to be made that both bullet points in the question should be an introduction in the source. Dec 29, 2014 at 19:02
• Most of the "hard coded numbers" (as used with arithmetic-shift and bitwise-and) are basically quicker ways to divide by 8 or 16 (and get the remainder). I could perhaps replace those with truncate/ (or quotient/remainder in Racket). I didn't actually write the code for "ease of understanding" per se, though I suppose I should make a greater effort to do so. e.g., m (in next-marked-from) could have been named mask, since it masks out the bits we don't care about when calling first-set-bit. Dec 29, 2014 at 19:42
• In any case, I wanted to ensure that I have the best algorithm in place before I write docstrings or doc comments, since those would have to be rewritten if the algorithm substantially changes. Dec 29, 2014 at 19:43
• Head down in the maths and the language documentation, things seem obvious. When the code gets browsed by an idiot like me on Github, it's full of magic numbers - divide-by-16 is a useful abstraction, particularly if the reason I'm looking at the code is because I am not up to speed on a bit-shifting sieve or the various related functions. For what it's worth, looking at the Racket code base, there appears to be a lot of code where comments were an unrealized afterthought. Dec 29, 2014 at 20:21
• Indeed, and I'm actually planning to split out the bitset into its own record type for better encapsulation and abstraction. Man, if this were Racket-specific, I'd just make a class and be done with it, but trying to be portable here. :-) Dec 29, 2014 at 20:26

To be honest, if it weren't for the hint in the first line (define (primes-up-to-about n)), I'd have no idea that this was an implementation of the Sieve of Eratosthenes.

The main problem, I think, is that names like table, clear!, and next-marked-from are all too low-level, not much better than the bytevector primitives, and devoid of contextual meaning.

To start, I think you should be able to rearrange the definitions.

  (define bytelen (+ (arithmetic-shift n -4) 1))
(define primes-table (make-bytevector bytelen 255))
(define bitlen (arithmetic-shift len 3))


Instead of clear!, try a more meaningful name like mark-as-composite!.

To complete the abstraction of the bytevector as a sieve, you would also need a helper function (define (next-prime-above n) …) that combines your next-marked-from and index->value.

(reverse rv) seems suspiciously expensive. I think you should be able to build the list of primes in ascending order using smarter recursion.

• Thanks for your review! 1. The original intent was to implement a bitset library, so the functions were named as bitset operations. I suppose it makes sense to use more domain-specific names. 2. I actually use the index from next-marked-from as the starting place for clearing the bitset flags. So with a next-prime-above refactor, I'd probably have to return two values, the index and the value. 3. reverse is a left-folding operation and is just as efficient as other left-folding operations like for-each and sum, and more efficient than right-folding operations like map and filter. Dec 28, 2014 at 23:05
• 3b. Accumulate-and-reverse is a pretty common idiom for simulating right-folding with a left-fold. The advantage of left-folding is that it's tail-recursive and thus doesn't eat up your call stack. I "could" make the function right-folding, which would eliminate the reverse, but it'd eat up a call frame for each prime number found, which could be problematic for huge numbers of primes. Dec 28, 2014 at 23:07