篩法のメモ
参考: 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しかないから通常の方が速いということらしい. 成る程.