Project Euler: Problem 10
The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17.
Find the sum of all the primes below two million.
http://projecteuler.net/index.php?section=problems&id=10
2,000,000 以下の自然数列の中に現われる素数の合計を求めよとのこと。
とりあえず、Problem 7 までに作った素数判定関数を使って深く考えずに Gauche で書いてみたが…死ぬほど遅い。
#!/usr/bin/env gosh ;; prime test of fermat algorithm (define (prime-fermat? n) (case n [(1) #f] [(2) n] [else (if (even? n) #f (= (remainder (expt 2 (- n 1)) n) 1))])) ; (gcd n 2) ==> 1 for all odd number except 1 ;; fast but inexact prime test: ;; divide number with known small primes (define (prime-inexact? n) ; fast but inexact (case n [(1) #f] [(2) n] [else (if (even? n) #f (let loop [(lis '(3 5 7 11 13 17 19 23 29 31))] ; magic of first 10 odd prime (if (pair? lis) (if (= n (car lis)) n (if (= (remainder n (car lis)) 0) #f (loop (cdr lis)))) n)))])) ;; slow but exact prime test ;; loop odd number from 3 until sqrt(number) and test if the number is dividable (define (prime-exact? n) ; slow but exact (case n [(1) #f] [(2) n] [else (if (even? n) #f (let loop [(i 3) ; minumum odd prime (m (floor->exact (sqrt n)))] (if (<= i m) (if (= (remainder n i) 0) #f (loop (+ i 2) m)) n)))])) (define (prime? n) ; try fast test first (and (prime-inexact? n) ; almost of non-primes are rejected here (prime-fermat? n) (prime-exact? n))) (define (main args) (define threshold 1000) (let loop [(n 3) ; start from first odd prime (r 2)] ; result (2 is a prime below n) (if (<= n 2000000) [begin (if (< threshold n) [begin (set! threshold (+ threshold 1000)) (display ".") ; to display progress (flush)]) (loop (+ n 2) ; loop only for odds (if (prime? n) (+ n r) r))] [begin (newline) (newline) (print #`"the answer is ,r")])) 0)
3 から 2,000,000 までの素数を数え上げて素数だったら加算していって表示。安直。実行してみたら以下のような感じ。
% time gosh problem10-1.scm ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ................................................................................. ....................................................... the answer is ************ gosh problem10-1.scm 3294.73s user 16.61s system 85% cpu 1:04:48.67 total
死ぬほど遅い。とりあえず答えは出たものの明かに何かが間違っていると思われてならない。
そんなわけで、Project Euler の解説を元にエラトステネスのふるい (Sieve of Eratosthenes) を使って作り直してみた。
#!/usr/bin/env gosh (use srfi-1) (define (prime-list n) ; sieve of eratosthenes (list version) (case n [(1) '()] [else (let loop [(i 0) (m (floor->exact (sqrt n))) (p (iota (- n 1) 2))] (if (< i m) (let* [(e (list-ref p i)) (p (filter (lambda (x) (or (= x e) ; 素数リストの先頭そのものであるか? (not (= (remainder x e) 0)))) ; 倍数でないか? p))] (loop (+ i 1) m p)) ; 素数の倍数を除いたものを持って次のループに p))])) (use gauche.collection) (define (prime-vector n) ; sieve of eratosthenes (vector version) ; ベクターを使うように作ってみたけど遅いしリストと併用してるし微妙 (case n [(1) '()] [else (let iloop [(i 0) (m (floor->exact (sqrt n))) (p (list->vector (iota (- n 1) 2)))] (if (< i m) (let jloop [(j 1) (m (vector-length p))] (if (and (vector-ref p i) (< j m)) [begin (if (and (not (= (+ j 2) (+ i 2))) (= (remainder (+ j 2) (+ i 2)) 0)) (vector-set! p j #f)) (jloop (+ j 1) m)] (iloop (+ i 1) m p))) (filter identity p)))])) (define (main args) (print #`"the answer is ,(fold + 0 (prime-list 2000000))") 0)
とりあえず実行してみた。
% time gosh problem10-2.scm the answer is ************ gosh problem10-2.scm 149.73s user 0.74s system 83% cpu 3:00.71 total
まだちょっと遅いものの 22 倍くらいは速い。最初から最後までまるまるメモリの上に乗っけてしまうのも無駄な気がしないでもないけど、イマドキそんな些細なことばっかり気にしていてもしょうがないですね。とりあえず、今回の教訓としては適切なアルゴリズムの選択は死ぬほど重要ってことですね。当たり前すぎるけど。
