Scheme:数遊び:素数列

素数列を求める

関連: Scheme:数遊び:素因数分解

エラトステネスの篩い

gemma(2008/03/21 13:32:37 PDT)追記
普通にリストで。

(use srfi-1)
;; n must be greater than 3.
(define (primes n)
  (let loop ((l (unfold (cut > <> n) values (cut + <> 2) 3))
             (prime-list '(2)))
    (let1 m (car l)
      (if (> (expt m 2) n)
        (append (reverse prime-list) l)
        (loop (remove (lambda (x) (zero? (modulo x m))) l) (cons m prime-list))))))

(primes 100) outputs (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)
ストリーム版(リスト版に比べて遅くなるのは避けられない)

(use util.stream)

(define (sieve xs)
  (stream-delay (stream-cons (stream-car xs)
                             (sieve (stream-remove (lambda (n)
                                                     (zero? (modulo n (stream-car xs))))
                                                   (stream-cdr xs))))))

(define prime (stream-cons 2 (sieve (stream-iota -1 3 2))))

(stream->list (stream-take prime 25)) outputs (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)

gemma(2006/06/19 20:15:30 PDT)

(use srfi-1)

(define (make-tousa start end step)
  (unfold (cut > <> end) values (cut + <> step) start))

(define (sieve-multiplies m l)
  (remove! (lambda (x)
             (eq? (modulo x m) 0))
           l))

(define (make-prime-list n)
  (let1 sqrt-n (sqrt n)
    (let loop ((l (make-tousa 2 n 1))
               (prime-list '()))
      (if (null? l)
          (reverse! prime-list)
          (let1 m (car l)
            (if (> m sqrt-n)
                (append! (reverse! prime-list) l)
                (loop (sieve-multiplies m l) (cons m prime-list))))))))

(define (prime? n)
  (cond 
   ((eq? n 1) #f)
   ((eq? n (car (last-pair (make-prime-list n)))) #t)
   (else #f)))

篩いにかける前に、make-tousaの時点であらかじめ偶数を除いておけば、少し速くできる。

(define (make-prime-list n)
  (case n
    ((1) '())
    ((2) (list 2))
    (else 
     (let1 sqrt-n (sqrt n)
       (let loop ((l (make-tousa 3 n 2))
                  (prime-list (list 2)))
         (let1 m (car l)
           (if (> m sqrt-n)
               (append! (reverse! prime-list) l)
               (loop (sieve-multiplies m l) (cons m prime-list)))))))))

skimu (2006/06/19 22:27:19 PDT): ちょっとトリッキーだけど、

(define (make-prime-list n)
  (let ((sqrt-n (sqrt n))
        (lis    (make-tousa 2 n 1)))
    (let loop ((l lis))
      (cond ((null? l) lis)
            ((> (car l) sqrt-n) lis)
            (else
             (loop (sieve-multiplies (car l) (cdr l))))))))

• SRFI:1 によると、

  Many procedures in this library have "pure" and "linear update" variants. A "pure" procedure has no side-effects, and in particular does not alter its arguments in any way. A "linear update" procedure is allowed -- but not required -- to side-effect its arguments in order to construct its result.

  なので、remove! が副作用なしで実装されているとまずくないですか？
  • skimu 突っ込みサンクス。

Shiro (2006/06/19 23:47:53 PDT):

短さを優先するならcomprehension。これはLL Ringのコメントにも書いたけど。

(use srfi-42)
(list-ec (: n 2 101) (if (not (any?-ec (: j 2 (sqrt n)) (zero? (modulo n j))))) n)

長さを決めないで、「最初のn個の素数」を求めるならstreamが使える。 ちょっとさぼって奇数のみのstreamからフィルタしている。

(use util.stream)
(use srfi-42)

(define (primes)
  (stream-filter (lambda (k)
                   (not (any?-ec (: j 3 (sqrt k) 2) (zero? (modulo k j)))))
                 (stream-cons 2 (stream-iota -1 3 2))))

gosh> (stream->list (stream-take (primes) 20))
(2 3 5 7 9 11 13 17 19 23 25 29 31 37 41 43 47 49 53 59)

いやまてよ、これは各数についてチェックをかけてるから篩とはちょっと違うか…

leque(2006/06/20 03:13:37 PDT) エラトステネスの篩の「sqrt(n) まで調べればよい」というのを使っていないけれど

(define primes
  (let1 f (lambda (strm)
            (stream-cons
             (stream-car strm)
             (stream-remove (lambda (k)
                              (zero? (remainder k (stream-car strm))))
                            (stream-cdr strm))))
    (stream-cons 2 (f (stream-iota -1 3 2)))))

gosh> (stream->list (stream-take primes 20))
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71)

skimu(2006/06/20 03:28:09 PDT): sqrt(n) まで判定すれば良いのは n が素数か判定するときなので n 個の素数を求めるのとはちと違います。

(define (prime? n)
  (let lp ((ps primes))
    (cond ((> (stream-car ps) (sqrt n)) #t)
          ((zero? (remainder n (stream-car ps))) #f)
          (else 
           (lp (stream-cdr ps))))))

しかし、shiro さんバージョンで平方数が消えないのはなぜ?

Shiro(2006/06/20 04:10:13 PDT): あっしまった。srfi-42のイテレーションは上限がexclusive だから、(+ (floor (sqrt n)) 1) までにしないとまずいか。

齊藤(2009/12/17 21:25:07 PST):

(use util.stream)

(define prime-numbers
  (stream-cons 2
    (stream-filter prime? (stream-iota -1 3 2))))

(define (prime? n)
  (let1 limit (floor (sqrt n))
    (not
     (stream-any
      (lambda(x) (zero? (modulo n x)))
      (stream-break (pa$ < limit) prime-numbers)))))

gemma(2009/12/18 05:03:36 PST): LtU:The Genuine Sieve of EratosthenesのHaskellコードをGaucheに移植したもの。これまでの素数ストリームよりも性能がよい。

(use util.stream)
    
(define (composites s)
  (stream-delay
   (let1 p (stream-car s)
     (stream-cons (expt p 2)
                  (xor (stream-map (pa$ * p) (xor (stream-cdr s) (composites s)))
                       (composites (stream-cdr s)))))))

(define (xor s0 s1)
  (stream-delay
   (let ((x (stream-car s0))
         (y (stream-car s1)))
     (cond
      ((= x y) (xor (stream-cdr s0) (stream-cdr s1)))
      ((< x y) (stream-cons x (xor (stream-cdr s0) s1)))
      ((> x y) (stream-cons y (xor s0 (stream-cdr s1))))))))

(define primes (stream-cons 2 (xor (stream-iota -1 3) (composites primes))))

Tags: Puzzle, util.stream