; pritchard's wheel sieve

(define (isqrt n)
  (if (not (and (positive? n) (integer? n)))
      (error 'isqrt "must be positive integer")
      (let loop ((x n))
        (let ((y (quotient (+ x (quotient n x)) 2)))
          (if (< y x) (loop y) x)))))

(define (eratosthenes n)
  (let* ((len (quotient (- n 1) 2))
         (bits (make-vector len #t)))
    (let loop ((i 0) (p 3) (ps (list 2)))
      (cond ((< n (* p p))
              (let loop ((i i) (p p) (ps ps))
                (cond ((= i len) (reverse ps))
                      ((vector-ref bits i)
                        (loop (+ i 1) (+ p 2) (cons p ps)))
                      (else (loop (+ i 1) (+ p 2) ps)))))
            ((vector-ref bits i)
              (let loop ((j (+ (* 2 i i) (* 6 i) 3)))
                (if (< j len)
                    (begin (vector-set! bits j #f)
                           (loop (+ j p)))))
              (loop (+ i 1) (+ p 2) (cons p ps)))
            (else (loop (+ i 1) (+ p 2) ps))))))

(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))))))

(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)))))

(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))))))

(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)))))))

(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))))))

(display (pritchard 100)) (newline)
(display (length (pritchard 1000))) (newline)

(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)
168