素数判定の高速化
素数判定を使うことが多いProject Eulerなので、Problem 6,Problem 10の回答例を見ながら試し割りの高速化を図ります。
試し割り編
今までのアルゴリズムは、
2以外の素数は奇数。
というルールに基づき、3から3,5,7,9,11,13...と偶数で割ってみて√2まで割り切れなければ素数。としてました。
回答例を読むと、
「2,3以外の素数は6k±1の中にある。」
6個づつ表にしてみると、
1 2 3 4 5 6
2 3 5
7 11
13 17
19 23
29ホントだ。綺麗に素数が現れる。
次の素数ループは30なので
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
2 3 5 7 11 13 17 19 23 29
31 37 41 43 47 53 59
61 67 71 73 79 83 89
2,3,5以外の素数は30k+{-23,-19,-17,-13,-11,-7,-1,1}の中にある。
30はちと複雑。しかも劇的に減るわけでもない。6±1は現実的なラインです。
こうやって素数表を作ってみるとオイラーが41の素数式を見つけられた訳がわかりますね。
というわけで修正。
(define (prime? n) (letrec ((e (sqrt n)) (trial-divisor (lambda (t) (or (> t e) (and (not (or (zero? (modulo n t)) (zero? (modulo n (+ t 2))))) (trial-divisor (+ t 6))))))) (and (> n 1) (or (<= n 3) (and (odd? n) (not (zero? (modulo n 3))) (trial-divisor 5)))))) (time (apply + (filter prime? (iota 2000000)))) ; cpu time: 7042 real time: 7043 gc time: 194 ; 142913828922 (time (apply + (filter old-prime? (iota 2000000)))) ; cpu time: 14802 real time: 14805 gc time: 22 ; 142913828922
200万までの素数を求めるのに、14秒かかっていたのに、7秒に短縮。半分になりました。
最大ループは √200万 / 飛ばし回数 なので、700回から、240回に削減。条件分岐が増えたので1/3とまではいきませんが相当効果アリっす。
試し割りでも結構早くなるもんだなぁと関心しました。
もうちょい高速化
2,3以外の素数は6k±1なので、あらかじめ6k±1だけ作っておく。という手法。
200万の母数を70万に絞れば高速化できる！！
どうかな・・・。
(define (make-primes n) (letrec ((iter (lambda (k) (if (>= k n) '() (cons k (if (>= (+ k 2) n) '() (cons (+ k 2) (iter (+ k 6))))))))) (filter prime? (cons* 2 3 (iter 5))))) (time (apply + (make-primes 2000000))) ; cpu time: 6236 real time: 6239 gc time: 35 ; 142913828922
7秒から,6.2秒。 0.8秒早くなりました。
び・微妙。
考えてみればあまり早くならないのは当然で、2,3の倍数の繰り返し回数は1。1を0にした事でちょっと早くなったけど、「巨大素数の判定には時間がかかる」。
素数判定の道はまだまだ遠いな。
エラトステネスのふるい
1から200万までの素数表を作り上げる。
ベクタを使って破壊的にチェックを入れていくらしい。
とりあえずナイーブにやってみる。
(define (sieve-of-eratothenes n) (letrec ((sieve (make-vector (+ n 1) #f)) (limit (sqrt n)) (sum (lambda (i) (if (> i n) 0 (+ (if (not (vector-ref sieve i)) i 0) (sum (+ i 1)))))) (check (lambda (i step) (if (<= i n) (begin (vector-set! sieve i #t) (check (+ i step) step))))) (iter (lambda (i) (cond ((> i limit) sieve) ((vector-ref sieve i) (iter (+ i 1))) (else (check (+ i i) i) (iter (+ i 1))))))) (iter 2) (sum 2))) (time (sieve-of-eratothenes 2000000)) ; cpu time: 877 real time: 877 gc time: 276 ; 142913828922
1秒！！
はやっ。
求める素数の最大値がわかっていれば素数表を作り上げて参照したほうが早い。
2,3の倍数を特殊化したらとんでもなく早いぞ。
エラトステネスのふるい(2の倍数の特殊化)
2の倍数を特殊化します。
テーブル数が半分になってメモリにもやさしい。
(define (sieve-of-eratothenes n) (letrec ((sieve (make-vector (quotient (+ n 1) 2) #f)) (limit (sqrt n)) (sum (lambda (i) (if (> i n) 0 (+ (if (not (vector-ref sieve (quotient i 2))) i 0) (sum (+ i 2)))))) (check (lambda (i step) (if (< i n) (begin (vector-set! sieve (quotient i 2) #t) (check (+ i step) step))))) (iter (lambda (i) (cond ((> i limit) sieve) ((vector-ref sieve (quotient i 2)) (iter (+ i 2))) (else (check (+ i i i) (+ i i)) (iter (+ i 2))))))) (iter 3) (+ 2 (sum 3)))) (time (sieve-of-eratothenes 2000000)) ; cpu time: 464 real time: 464 gc time: 109 ; 142913828922
200万までの素数を全て判定できるテーブルを作ったのに、0.5秒。
あ・り・え・な・い。
充分過ぎるくらい早いっす。
ライブラリ化
テーブルをキャッシュ化しておいて高速化を狙います。
(define (make-sieve-of-eratothenes n) (letrec ((sieve (make-vector (quotient (+ n 1) 2) #t)) (limit (sqrt n)) (check (lambda (i step) (if (< i n) (begin (vector-set! sieve (quotient i 2) #f) (check (+ i step) step))))) (make-table (lambda (i) (cond ((> i limit) sieve) ((not (vector-ref sieve (quotient i 2))) (make-table (+ i 2))) (else (check (+ i i i) (+ i i)) (make-table (+ i 2)))))) (prime? (lambda (i) (or (= i 2) (and (> i 1) (odd? i) (vector-ref sieve (quotient i 2)))))) (fold (lambda (i f acc) (if (>= i n) acc (fold (+ i 2) f (if (vector-ref sieve (quotient i 2)) (f i acc) acc)))))) (make-table 3) (lambda (m) (case m ('prime? prime?) ('fold (lambda (f acc) (fold 3 f (f 2 acc)))))))) (define *sieve* (make-sieve-of-eratothenes 2000000)) (define (prime? i) ((*sieve* 'prime?) i)) (define (sieve-fold f acc) ((*sieve* 'fold) f acc)) (time (sieve-fold + 0)) ; cpu time: 96 real time: 96 gc time: 0 ; 142913828922
前と変わらずですが、素数テーブル生成に0.5秒。
キャッシュを使っているので当たり前に早いです。
あらかじめ素数の最大値を割り出す必要がありますが、かなり使えそうです。
改造しながら高めていこう。
