Gaucheでエラトステネスの篩で素数ストリーム(世間一般のやつより速い)

Haskellなどで人気の1行エラトステネスの篩は、アルゴリズム本来の計算量より遅いらしい。

sieve (p : xs) = p : sieve [x | x <- xs, x ‘mod‘ p > 0]

• LtU:The Genuine Sieve of Eratosthenes
• Wikipedia:en:Sieve of Eratosthenes

それを受けて、LtUで以下の速いHaskellコードがあった。

primes = 2 : xor [3..] (composites primes)

composites (p:ps) = cs
    where cs = (p*p) : xor (map (p*) (xor ps cs)) (composites ps)

xor (x:xs) (y:ys)  | x==y          = xor xs ys
                   | x<y           = x : xor xs (y:ys)
                   | x>y           = y : xor (x:xs) ys

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

実行。25番めの素数は97。

(print (stream->list (stream-take primes 25)))
=> (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)

Project Eulerに使えるかも。