Purely functional Sieve of Eratosthenes

Question

Many implementations of Sieve of Eratosthenes (used to find prime numbers up to a given n) use a temporary mutable array to keep track of which numbers are composites. I'm looking to write a purely functional version that has the same runtime complexity:

(define (primes-up-to n)
  (cons 2 (unfold null?
                  car
                  (compose remove-odd-multiples car+cdr)
                  (iota (quotient (- n 1) 2) 3 2))))

(define (remove-odd-multiples n xs)
  (let recur ((xs xs)
              (i (* n n)))
    (cond ((null? xs) '())
          ((< (car xs) i) (cons (car xs) (recur (cdr xs) i)))
          ((< i (car xs)) (recur xs (+ i n n)))
          (else (recur (cdr xs) (+ i n n))))))

(Should work on any Scheme implementation that supports SRFI 1. If your implementation doesn't provide compose, use (lambda (x) (remove-multiples (car x) (cdr x))) instead.)

What I'm hoping reviewers can help with: Is there a way to make remove-odd-multiples¹ less verbose (perhaps using higher order functions) or more efficient², without increasing the runtime complexity (so any kind of filter or remove involving non-constant-time predicates is no go)?

---

¹ Since xs does not contain even numbers, there's no point dropping even multiples, hence remove-odd-multiples. The only difference from a more general remove-multiples is that it advances i by 2n instead of n each time.

² I've done some casual timing tests on (primes-up-to 100000) on DrRacket 6.5 and the right-folding version you see above is faster than a left-folding version. (I mention this preemptively because many Scheme answerers on Stack Overflow lean towards left-folding solutions.) Here's the left-folding version I time-tested with:

(define (primes-up-to n)
  (reverse (unfold-right null?
                         car
                         (compose remove-odd-multiples car+cdr)
                         (iota (quotient (- n 1) 2) 3 2)
                         '(2))))

(define (remove-odd-multiples n xs)
  (let loop ((result '())
             (xs xs)
             (i (* n n)))
    (cond ((null? xs) (reverse result))
          ((< (car xs) i) (loop (cons (car xs) result) (cdr xs) i))
          ((< i (car xs)) (loop result xs (+ i n n)))
          (else (loop result (cdr xs) (+ i n n))))))

Answer 1

Here's a small refinement on the technique. Wherin instead of looking for odd multiples and removing them, you keep an extra reference in the incomplete sieve and look those specific multiples of the next prime. The first step is the same, looking for 3* (3, 5, 7, 9, ...,) but the second step you only look for 5*(5, 7, 11, 13, 17, 19, 23 , 25, ...)

Overall a small speedup, but the upside is you will never be looking for a multiple that isn't either in the partial sieve or beyond it's range.

(define (primes-up2 n)
  (cons 2 (unfold null?
                  car
                  (lambda (x) (remove-partial-seive-multiples (car x) (cdr x)))
                  (iota (quotient (- n 1) 2) 3 2))))

(define (remove-partial-seive-multiples n partials)
  (let recur ((xs partials)
              (ys partials)
              (i (* n n)))
    (cond ((null? xs) '())
          ((< (car xs) i) (cons (car xs) (recur (cdr xs) ys i)))
          ((= (car xs) i) (recur (cdr xs) (cdr ys) (* (car ys) n)))
          (else (error "you ought not hit this")))))

This works as any multiple of n still in the partial sieve must also be coprime with all primes less than N. (of which the partial sieve is a complete list)

Now while not purely functional instead of keeping the result and reversing it in the end on your second iteration, you could implement it as a tail-modulo cons. While not stickly functional on the inside, it is functional from the outside.

Another optimization would be to pass the maximum canidate to the remover of multiples and just return the xs when i is smaller than the last search.