primes.eratosthenes-version.scm

;;
;; Gauche のジェネレータと区間篩を使った無限素数列
;;

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

(define isprime 1)
(define iscomposite 0)
(define segment-size #e1e5)
(define segment-start (ceiling->exact (expt segment-size 0.6))) ;; sqrt より少し大きめに取っておく

;; 初項 a0 交差 d の等差数列で lower-bound 以上となる初めての項を返す
(define (take-first-term a0 d lower-bound)
  (if (<= lower-bound a0)
      a0
      (+ a0 (* d (div (- lower-bound a0 (- d) 1) d)))))

;; 区間 [start, end] を small-primes によって篩に掛ける
;; また、small-primes には番兵として (sqrt end) より大きな素数を 1 つ要求する
(define (segment-sieve! bytevec start end small-primes)
  (let ([root (floor->exact (sqrt end))]
        [segment-size (- end start -1)])
    (let sieve! ([ps small-primes])
      (receive (p qs) (car+cdr ps)
        (when (<= p root)
	  (do-ec (:range i (- (take-first-term (* p p) p start) start) segment-size p)
	    (u8vector-set! bytevec i iscomposite))
	  (sieve! qs))))))

;; 篩に掛け終えた bytevec から素数を収集する
(define (collect-primes bytevec start end)
  (let recur ([i 0] [n start])
    (cond
     [(< end n) (list)]
     [(= iscomposite (u8vector-ref bytevec i)) (recur (+ i 1) (+ n 1))]
     [else (cons n (recur (+ i 1) (+ n 1)))])))

;; 区間 [start, end] を篩に掛け素数を列挙し、素数列を返す
;; segment-sieve! と同じく、small-primes には番兵として (sqrt end) より大きな素数を 1 つ要求する
(define (get-primes/segment start end small-primes)
  (let1 bytevec (make-u8vector (- end start -1) isprime)
    (segment-sieve! bytevec start end small-primes)
    (collect-primes bytevec start end)))

;; N までの素数列を返す。性能は殆ど要求されないので試し割りでも良い
(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 (- segment-start 1)))
	    (lambda () (prime-generator)))))

;; 無限素数列と区間篩を使った素数ジェネレータ
(define prime-generator
  (generate
   (lambda (yield)
     (let recur ([start segment-start] [ps (list)])
       (dolist (p ps) (yield p))
       (recur (+ start segment-size)
	      (get-primes/segment start (+ start segment-size -1) primes))))))

;; 試しに 100 万と 1 番目の素数を持ってくる
(print (list-ref primes #e1e6))