篩法のメモ

参考: d:id:odz:20080229:1204287954
Project Euler用のメモ
100万までの各数について真の約数の和を求めることを考える.
C系だとこんな感じ.

long[] a;
a.length = 1000000;
for (long i = 2; i < 1000000; i++){
  for (j = 2; j * i < 1000000; j++){
     a[j * i] = a[j * i] + i;
  }
}

gaucheで篩をどうやって書いたものかと思って最初に書いたのが以下.

(define (divisors-sieve n)
  (define (divisors-sieve-in m lis)
    (let loop ((i 1) (cur lis) (res '()))
      (if (null? cur) (reverse res)
        (if (and (not (= i m))
                 (= 0 (remainder i m)))
            (loop (+ i 1) (cdr cur) (cons (+ (car cur) m) res))
          (loop (+ i 1) (cdr cur) (cons (car cur) res))))))
  (let loop ((i 2) (lis (make-list n 1)))
    (if (> i n) lis
      (loop (+ i 1) (divisors-sieve-in i lis)))))
;(time (divisors-sieve 10000))
; real  43.203
; user  40.094
; sys    3.000
;(1 1 1 3 1 6 1 7 4 8 1 16...

当然のごとく遅い. こういうときのためにvectorがあると思い出して書き直した.

(define (divisors-sieve-by-vector n vec)
  (let loop ((i 2))
    (if (> i n) vec
      (let loop-in ((j 2))
        (if (> (* i j) n)
            (loop (+ i 1))
          (begin
            (vector-set! vec (* i j) (+ (vector-ref vec (* i j)) i))
            (loop-in (+ j 1))))))))

(time (divisors-sieve-by-vector 10000 (make-vector 10001 1)))
; real   0.016
; user   0.015
; sys    0.000
; #(1 1 1 1 3 1 6 1 7 4 8 1 16 ...

(time (vector-length (divisors-sieve-by-vector 1000000 (make-vector 1000001 1))))
; real   5.016
; user   5.016
; sys    0.000
1000001

速いっていいことですよね. Project Euler自体はちょくちょく時間を見つけては解いていて, 現在89問.

(define (divisors-sieve-by-vector n)
  (let1 vec (make-vector (+ n 1) 1)
        (let loop ((i 2))
          (if (> i n) vec
            (let loop-in ((j 2))
              (if (> (* i j) n)
                  (loop (+ i 1))
                (begin
                  (vector-set! vec (* i j) (+ (vector-ref vec (* i j)) i))
                  (loop-in (+ j 1)))))))))

vectorを引数に取る前に内部で作ってそれを返すべきだろう.

もう一個メモ

;通常のvector
(define (prime-sieve-by-vector n vec)
  (let loop ((i 2))
    (cond ((> i n) vec)
          ((= 1 (vector-ref vec i))
           (let loop-in ((j 2))
             (if (> (* i j) n)
                 (loop (+ i 1))
               (begin
                 (vector-set! vec (* i j) 0)
                 (loop-in (+ j 1))))))
          (else
           (loop (+ i 1))))))
;(time (vector-length (prime-sieve-by-vector (expt 10 7) (make-vector (+ ...
; real   8.375
; user   8.281
; sys    0.031

;(time (vector-length (prime-sieve-by-vector (expt 10 8) (make-vector (+ ...
;out of memory (400000012).  aborting...

;s8vectorを使用
(define (prime-sieve-by-s8vector n vec)
  (let loop ((i 2))
    (cond ((> i n) vec)
          ((= 1 (s8vector-ref vec i))
           (let loop-in ((j 2))
             (if (> (* i j) n)
                 (loop (+ i 1))
               (begin
                 (s8vector-set! vec (* i j) 0)
                 (loop-in (+ j 1))))))
          (else
           (loop (+ i 1))))))
;(time (s8vector-length (prime-sieve-by-vector (expt 10 7) (make-s8vecto ...
; real  13.891
; user  13.625
; sys    0.031

0と1しかないから通常の方が速いということらしい. 成る程.