primes.odd-sieve-version.scm

;;
;; odd sieve + segment sieve + generator -> primes-lseq
;;

(use srfi-1)
(use srfi-42)
(use gauche.uvector)
(use gauche.generator)

(define isprime 1)
(define iscomposite 0)

(define (take-first-term a0 d lower-bound)
  (if (<= lower-bound a0)
      a0
      (+ a0 (* d (div (- lower-bound a0 (- d) 1) d)))))
(define (->odd x)
  (if (odd? x) x (+ x 1)))
(define (->even x)
  (if (even? x) x (+ x 1)))
(define (integer->index x)
  (ash x -1))

(define segment-size (->even #e1e5))
(define first-segment-start (->odd (ceiling->exact (expt segment-size 0.6))))
(define first-segment-end (+ first-segment-start segment-size -1))
(define vecsize (/ segment-size 2))

(define (segment-sieve! bytevec start end small-primes)
  (let ([root (floor->exact (sqrt end))]
        [start-index (integer->index start)])
    (let sieve! ([ps (cdr small-primes)])
      (receive (p qs) (car+cdr ps)
        (when (<= p root)
	  (do-ec (:let i (- (take-first-term (integer->index (* p p)) p start-index) start-index))
		 (:range j i vecsize p)
	    (u8vector-set! bytevec j iscomposite))
	  (sieve! qs))))))

(define (foreach-primes proc bytevec start end)
  (let recur ([i 0] [n start])
    (when (< n end)
      (when (= isprime (u8vector-ref bytevec i))
	(proc n))
      (recur (+ i 1) (+ n 2)))))

(define (segment-sieve&foreach-primes proc start end small-primes)
  (let1 bytevec (make-u8vector vecsize isprime)
    (segment-sieve! bytevec start end small-primes)
    (foreach-primes proc bytevec start end)))

(define (get-primes N)
  (list-ec (:range n 2 (+ N 1))
    (if (every?-ec (:range j 2 (+ 1 (floor->exact (sqrt n))))
	  (< 0 (mod n j))))
    n))

(define primes
  (generator->lseq
   (gappend (list->generator (get-primes (- first-segment-start 1)))
	    (lambda () (prime-generator)))))

(define prime-generator
  (generate
   (lambda (yield)
     (let recur ([start first-segment-start] [end first-segment-end])
       (segment-sieve&foreach-primes yield start end primes)
       (recur (+ start segment-size) (+ end segment-size))))))