; segmented sieve of eratosthenes

(define (isqrt n)
  (let loop ((x n) (y (quotient (+ n 1) 2)))
    (if (<= 0 (- y x) 1) x
      (loop y (quotient (+ y (quotient n y)) 2)))))

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

(define (primes-range l r b)
  (let* ((ps (cdr (primes (+ (isqrt r) 1))))
         (qs (map (lambda (p) (modulo (* -1/2 (+ l 1 p)) p)) ps))
         (zs '()) (z (lambda (p) (set! zs (cons p zs)))))
    (do ((t l (+ t b b))
         (qs qs (map (lambda (p q) (modulo (- q b) p)) ps qs)))
        ((= t r) (reverse zs))
      (let ((bs (make-vector b #t)))
        (do ((qs qs (cdr qs)) (ps ps (cdr ps))) ((null? qs))
          (do ((j (car qs) (+ j (car ps)))) ((<= b j))
            (vector-set! bs j #f)))
        (do ((j 0 (+ j 1))) ((= j b))
          (if (vector-ref bs j) (z (+ t j j 1))))))))

(display (length (primes #e1e6))) (newline)
(display (length (primes #e2e6))) (newline)
(display (- 148933 78498)) (newline)
(display (length (primes-range #e1e6 #e2e6 #e1e5))) (newline)

78498
148933
70435
70435