Scheme

2009年11月23日 (月)

Project Euler - Problem 3 : 素因数分解

問題はこちらをご覧ください。


この問題は要するに素因数分解の問題です。自作の手続き "factorize" を使うのなら、

(define problem-003-1 (lambda (n) (caar (last-pair (factorize n))))) > (time (problem-003-1 600851475143)) (time (problem-003-1 600851475143)) no collections 1 ms elapsed cpu time 1 ms elapsed real time 24272 bytes allocated

これだけで終わってしまいます。

でも、わざわざ素因数分解の手続きを作らなくても、2 で繰り返し割った後、奇数の小さい方から順に割り続けるだけでも答えは出ます。

(define problem-003-2 (lambda (n) (define divide (lambda (n i) (cond ([= n 1] i) ([zero? (remainder n i)] (divide (/ n i) i)) (else n)))) (let loop ([n (divide n 2)] [i 3]) (if [> (* i i) n] n (loop (divide n i) (+ i 2)))))) > (time (problem-003-2 600851475143)) (time (problem-003-2 600851475143)) no collections 1 ms elapsed cpu time 1 ms elapsed real time 144 bytes allocated

2009年11月 9日 (月)

Project Euler - Problem 2 : 再帰的アルゴリズムと反復的アルゴリズム

問題はこちらをご覧ください。


この問題は、Scheme のコードを載せたことがあるのですが、改めて……

フィボナッチ数を求める場合、一番単純なのは定義に基づいて再帰的に求める方法です。(以下のコードはすべて Chez Scheme です)

;;; n 番目の Fibonacci 数を再帰的に求める。 (define fibonacci (lambda (n) (cond ([= n 1] 1) ([= n 2] 2) (else (+ (fibonacci (- n 2)) (fibonacci(- n 1)))))))

この手続きを使って、二つの方法を考えてみました。

;;; n 未満の偶数の Fibonacci 数の合計を出す。 (define problem-002-1 (lambda (n) (let loop ([i 1] [sum 0]) (let ([f (fibonacci i)]) (cond ([>= f n] sum) ([even? f] (loop (+ i 1) (+ sum f))) (else (loop (+ i 1) sum))))))) ;; > (time (problem-002-1 4000000)) ;; (time (problem-002-1 4000000)) ;; no collections ;; 3344 ms elapsed cpu time ;; 3411 ms elapsed real time ;; 2432 bytes allocated ;;; n 未満の Fibonacci 数列を作り、偶数だけを合計する。 (define problem-002-2 (lambda (n) (let loop ([i 1] [result '()]) (let ([f (fibonacci i)]) (if [>= f n] (apply + (filter even? result)) (loop (+ i 1) (cons f result))))))) ;; > (time (problem-002-2 4000000)) ;; (time (problem-002-2 4000000)) ;; no collections ;; 3324 ms elapsed cpu time ;; 3389 ms elapsed real time ;; 2848 bytes allocated

フィボナッチ数を求める場合、再帰的なアルゴリズムでは同じ計算を繰り返し行うので、非常に遅くなります。

実際、どちらも答えが出るまで約 3 秒かかっています。


もう一つの方法として、反復的なアルゴリズムでフィボナッチ数列を作ってから、そのリストを処理するというものがあります。

;;; n 未満の Fibonacci 数列を反復的に生成する。 (define fibonacci-list (lambda (n) (let loop ([a 1] [b 1] [result '()]) (if [>= a n] result (loop b (+ a b) (cons a result)))))) ;;; n 未満の Fibonacci 数列を求め、偶数だけを合計する。 (define problem-002-3 (lambda (n) (apply + (filter even? (fibonacci-list n))))) ;; > (time (problem-002-3 4000000)) ;; (time (problem-002-3 4000000)) ;; no collections ;; 0 ms elapsed cpu time ;; 0 ms elapsed real time ;; 496 bytes allocated

この場合、同じ計算を繰り替えさないので、非常に速く答えが出ます。

なお、"fibonacci-list" は内部処理を末尾再帰にしたかったのと、答えを出すのにリスト内の順序は関係なかったので、返ってくるリストはあえて「降順」のままになっています。


以前も書きましたが、繰り返しの処理には、"do" や "内部定義"、"letrec" などいろいろな方法があります。

でも、今回は自分の好みで、"named let" ばかりになってしまいました。たぶん、今後もそうなるんだろうなぁ……。

2009年11月 1日 (日)

Project Euler - Problem 1 : 今度は Scheme で……

問題はこちらをご覧ください。


この問題はいろいろな解き方があります。

例えば繰り返しを使って和を求めるもの。

Scheme では繰り返しは再帰を使って表現するのですが、それにもいくつか方法があります。

(なお、以下の手続きはすべて、「n 未満の和」を求めるものです。(problem-001-1 1000) といった具合に、引数に "1000" を与えるとこの問題の答えが出ます。また、自作の手続きについてはこちらをご覧ください)

内部定義を使うもの。

;; 内部定義を使って繰り返し(末尾再帰) (define problem-001-1 (lambda (n) (define iter (lambda (m sum) (cond ([>= m n] sum) ([or (zero? (remainder m 3)) (zero? (remainder m 5))] (iter (+ m 1) (+ sum m))) (else (iter (+ m 1) sum))))) (iter 1 0)))

"retlec" や "named let" を使っても同じことができます。

;; letrec を使って繰り返し (define problem-001-2 (lambda (n) (letrec ([iter (lambda (m sum) (cond ([>= m n] sum) ([or (zero? (remainder m 3)) (zero? (remainder m 5))] (iter (+ m 1) (+ sum m))) (else (iter (+ m 1) sum))))]) (iter 1 0))))

;; named let を使って繰り返し (define problem-001-3 (lambda (n) (let iter ([m 1] [sum 0]) (cond ([>= m n] sum) ([or (zero? (remainder m 3)) (zero? (remainder m 5))] (iter (+ m 1) (+ sum m))) (else (iter (+ m 1) sum))))))

条件に合うリストを作っておいて、それを加工していく方法もあります。

;;; list を加工して計算する 1 (define problem-001-5 (lambda (n) (apply + (filter (lambda (x) (or (zero? (remainder x 3)) (zero? (remainder x 5)))) (seq 1 (- n 1))))))

これは 1 から n までの数を並べたリストから条件に合う数だけを抜き出して、足していく手続きです。("seq" は等差数列を作る自作の手続きです)


;;; list を加工して計算する 2 (define problem-001-6 (lambda (n) (let ([lst1 (seq 0 (- n 1) 3)] ; 3 の倍数のリスト [lst2 (seq 0 (- n 1) 5)]) ; 5 の倍数のリスト (apply + (unique (append lst1 lst2))))))

これは予め「3 の倍数のリスト」と「5 の倍数のリスト」を連結して、重複したものを取り除いたうえで和を求めています。("unique" はリストの要素から重複したものを取り除く自作の手続きです)


ちなみに、こんな解き方もあります。

(define problem-001-7 (lambda (n) (let* ([m (- n 1)] [sum1 (apply + (seq 0 m 3))] ; 3 の倍数の和 [sum2 (apply + (seq 0 m 5))] ; 5 の倍数の和 [sum3 (apply + (seq 0 m 15))]) ; 15 の倍数の和 (- (+ sum1 sum2) sum3))))

一つの問題でもアプローチの仕方やコードの書き方にはいろいろな方法があります。解き方のバリエーションを考えるのもなかなか面白いものです。

2009年10月31日 (土)

自作の手続き for Project Euler

今まで Ruby で解いてきた Project Euler の問題を今度は Scheme で解いていくことにしました。

以前、いくつかの問題を Scheme で解いたものをピックアップして投稿したことがありましたが、Ruby で解いたときの経験を生かして、改めてコードを書いてみようと思っています。

基本的なアルゴリズムは、Ruby で解いたときと同じものになっていくと思うので、アルゴリズムの説明は Ruby 版の記事を読んでもらうことになると思います。


まずは、Project Euler を解くために使っていく自作の手続きをまとめて書いておきます。

以前にもほぼ同じものを投稿してあるのですが、いじっているうちに内容が少し変わってしまった手続きもあるので、今回改めて紹介しておきます。

なお、一部のコードはこちらのホームページのコード使用させていただいています。

また、Chez Scheme で独自に拡張された手続きやマクロを一部使用していますが、不明な点はこちらの "Chez Scheme Version 7 User's Guide" をご覧ください。(なお、"case-lambda" に関してはここにも情報があります)


;;;; ;;;; math.ss --- 数学関連の関数群 ;;;; ;;;; for Chez Scheme by Tsumuji ;;;; ;;; ;;; (sequence start end) ;;; (sequence start end step) ;;; ;;; * start 〜 end の範囲の等差数列のリストを返す。 ;;; * step が省略された場合、step = 1 とみなされる。 ;;; * 引数が少数の場合、計算誤差がでる場合があるので注意。 ;;; (define sequence (case-lambda ([start end] (sequence start end 1)) ([start end step] (if [or (zero? step) (< (* (- end start) step) 0)] '() (let ([op (if [> step 0] > <)]) (let loop ([val start] [result '()]) (if [op val end] (reverse result) (loop (+ val step) (cons val result))))))))) (define seq sequence) ; alias ;;; ;;; (prime? n) ;;; ;;; * n が「素数」なら #t ;;; (define prime? (lambda (n) (cond ([= n 1] #f) ([< n 4] #t) ([or (even? n) (zero? (remainder n 3))] #f) (else (let ([i-max (isqrt n)]) (let loop ([i 5] [add 2]) (cond ([> i i-max] #t) ([zero? (remainder n i)] #f) ([> i 10000] (fast-prime? n)) (else (loop (+ i add) (- 6 add)))))))))) ;;; ;;; (fast-prime? n) ;;; ;;; * n が「素数」なら #t ;;; * Miller-Rabin test によって判定している。 ;;; (define fast-prime? (lambda (n) (define miller-rabin (lambda (b s t) (if [= (expt-mod b t n) 1] #t (let ([n (- n)]) (let loop ([m (expt-mod b t n)] [i 0]) (cond ([= i s] #f) ([= m -1] #t) (else (loop (expt-mod m 2 n) (+ i 1))))))))) (cond ([< n 2] #f) ([= n 2] #t) ([even? n] #f) (else (let loop1 ([t (- n 1)] [s 0]) (if [even? t] (loop1 (/ t 2) (+ s 1)) (let loop2 ([i 5]) ; 5 回テストする。 (cond ([zero? i] #t) ([miller-rabin (+ 2 (random (- n 2))) s t] (loop2 (- i 1))) (else #f))))))))) ;;; ;;; (prime-list n) ;;; ;;; * 自然数 n 以下の素数のリストを返す。 ;;; * Ruby の"sample/sieve.rb" のアルゴリズムを改良したものを使用 ;;; している。 ;;; (define prime-list (lambda (n) (if [= n 2] '(2) ;; エラトステネスの篩 (let* ([p-vct (list->vector (cons 0 (sequence 3 n 2)))] ;; p-vct = #(3 5 7 ...) [size (vector-length p-vct)] [max (isqrt n)]) (let loop1 ([i 1]) (let ([j (vector-ref p-vct i)]) (cond ([zero? j] (loop1 (+ i 1))) ([> j max] (cons 2 (remv 0 (vector->list p-vct)))) (else (let loop2 ([k (* 2 i (+ i 1))]) (if [>= k size] (loop1 (+ i 1)) (begin (vector-set! p-vct k 0) (loop2 (+ j k))))))))))))) ;;; ;;; (factorize n) ;;; ;;; * 自然数 n を素因数分解した要素のリストを返す。 ;;; * 引数が 1 の場合、本来はエラーとすべきだが……。 ;;; ex : (factorize 12) => ((2 . 2) (3 . 1)) ;;; ex : (factorize 1) => ((1 . 0)) ;;; (define factorize (lambda (n) (define result '()) (define f-count (lambda (n i) (let loop ([n n] [c 0]) (cond ([zero? (remainder n i)] (loop (/ n i) (+ c 1))) ([> c 0] (set! result (cons (cons i c) result)) n) (else n))))) (cond ([= n 1] '((1 . 0))) (else (set! n (f-count n 2)) (set! n (f-count n 3)) (let loop ([n n] [i 5] [add 2]) (cond ([= n 1] (reverse result)) ([> (* i i) n] (reverse (cons (cons n 1) result))) (else (loop (f-count n i) (+ i add) (- 6 add))))))))) ;;; ;;; (divisor n) ;;; ;;; * 自然数 n の約数のリストを返す。 ;;; * n が大きくなると iota がエラーを出すことがある。 ;;; (define divisor (lambda (n) (let* ([a (filter (lambda (x) (zero? (remainder n x))) (cdr (iota (+ (isqrt n) 1))))] [b (reverse (map (lambda (x) (/ n x)) a))]) (if [= n (expt (car b) 2)] (append a (cdr b)) (append a b))))) ;;; ;;; (divisor2 n) ;;; ;;; * 自然数 n の約数をリストを返す。 ;;; * 素因数分解をもとに約数を探す。 ;;; (define divisor2 (lambda (n) ;; ex : (2 . 3) => (1 2 4 8) (define spread (lambda (lst) (let ([b (car lst)] [i-lst (iota (+ 1 (cdr lst)))]) (map (lambda (x) (expt b x)) i-lst)))) ;; ex : (1 2 4) (1 3) => (1 3 2 6 4 12) (define distribute (lambda (lst1 lst2) (apply append (map (lambda (x) (map (lambda (y) (* x y)) lst2)) lst1)))) (let ([lst (map spread (factorize n))]) (sort < (fold-left distribute '(1) lst))))) ;;; ;;; (count-divisor n) ;;; ;;; * 自然数 n の約数の個数を素因数分解を基にして数える。 ;;; (define count-divisor (lambda (n) (apply * (map (lambda (lst) (+ 1 (cdr lst))) (factorize n))))) ;;; ;;; (factorial n) ;;; (factorial n m) ;;; ;;; * 引数が 1 個の場合 : n! ;;; * 引数が 2 個の場合 : 下降階乗冪 ;;; (define factorial (case-lambda ([n] (factorial n n)) ([n m] (apply * (sequence (- n m -1) n))))) ;;; ;;; (perm n r) ;;; ;;; * 順列 nPr ;;; (define perm (lambda (n r) (factorial n r))) ;;; ;;; (permutations lst) ;;; (permutations lst size) ;;; ;;; * 順列をリストにして返す。 ;;; ex : (permutations '(1 2 3)) ;;; => ((1 2 3) (1 3 2) (2 1 3) (2 3 1) (3 1 2) (3 2 1)) ;;; ex : (permutations '(1 2 3) 2) ;;; => ((1 2) (1 3) (2 1) (2 3) (3 1) (3 2)) ;;; ;;; from 『M.Hiroi's Home Page』, modified by Tsumuji ;;; (define permutations (case-lambda ([lst] (permutations lst (length lst))) ([lst size] (define perm (lambda (lst a b size) (if [zero? size] (cons (reverse a) b) (fold-right (lambda (x y) (perm (remove x lst) (cons x a) y (- size 1))) b lst)))) (perm lst '() '() size)))) ;;; ;;; (permutations-for-each proc lst) ;;; (permutations-for-each proc lst size) ;;; ;;; * 順列を生成し、proc を適用する。 ;;; ;;; from 『M.Hiroi's Home Page』, modified by Tsumuji ;;; (define permutations-for-each (case-lambda ([proc lst] (permutations-for-each proc lst (length lst))) ([proc lst size] (define perm (lambda (lst a size) (if [zero? size] (proc (reverse a)) (for-each (lambda (n) (perm (remove n lst) (cons n a) (- size 1))) lst)))) (perm lst '() size)))) ;;; ;;; (comb n r) ;;; ;;; * 組合せ nCr ;;; (define comb (lambda (n r) (/ (factorial n r) (factorial r)))) ;;; ;;; (combinations ls n) ;;; ;;; * 組み合わせをリストにして返す。 ;;; ex : (combinations '(1 2 3 4) 2) ;;; => ((1 2) (1 3) (1 4) (2 3) (2 4) (3 4)) ;;; ;;; from 『M.Hiroi's Home Page』 ;;; (define combinations (lambda (ls n) (define comb (lambda (n ls a b) (cond ([zero? n] (cons (reverse a) b)) ([= (length ls) n] (cons (append (reverse a) ls) b)) (else (comb (- n 1) (cdr ls) (cons (car ls) a) (comb n (cdr ls) a b)))))) (if [> n (length ls)] #f (comb n ls '() '())))) ;;; ;;; (combinations-for-each proc ls n) ;;; ;;; * 組合せを生成し、proc を適応する。 ;;; ;;; from 『M.Hiroi's Home Page』 ;;; (define combinations-for-each (lambda (proc ls n) (define comb (lambda (ls n a) (cond ([zero? n] (proc (reverse a))) ([= (length ls) n] (proc (append (reverse a) ls))) (else (comb (cdr ls) (- n 1) (cons (car ls) a)) (comb (cdr ls) n a))))) (if [> n (length ls)] #f (comb ls n '())))) ;;; ;;; (real->integer n) ;;; ;;; * 実数の少数部分を切り捨てて、整数に変換する。 ;;; (define real->integer (lambda (n) (exact (truncate n)))) ;;; ;;; (triangle-number n) ;;; ;;; * n 番目の三角数 ;;; * Tn = n * (n + 1) / 2 ;;; (define triangle-number (lambda (n) (/ (+ (* n n) n) 2))) ;;; ;;; (triangle-number? n) ;;; ;;; * n が三角数なら「何番目か」を返す。 ;;; (define triangle-number? (lambda (n) (let* ([a (sqrt (+ 1 (* 8 n)))] [b (/ (+ -1 a) 2)]) (if [positive-integer? b] b #f)))) ;;; ;;; (square-number n) ;;; ;;; * n 番目の四角数(平方数) ;;; * Sn = n * n ;;; (define square-number (lambda (n) (* n n))) ;;; ;;; (square-number? n) ;;; ;;; * n が四角数(平方数)なら「何番目か」を返す。 ;;; (define square-number? (lambda (n) (let ([a (sqrt n)]) (if [positive-integer? a] a #f)))) ;;; ;;; (pentagon-number n) ;;; ;;; * n 番目の五角数 ;;; * Pn = n * (3 * n - 1) / 2 ;;; (define pentagon-number (lambda (n) (/ (- (* 3 n n) n) 2))) ;;; ;;; (pentagon-number? n) ;;; ;;; * n が五角数なら「何番目か」を返す。 ;;; (define pentagon-number? (lambda (n) (let* ([a (sqrt (+ 1 (* 24 n)))] [b (/ (+ 1 a) 6)]) (if [positive-integer? b] b #f)))) ;;; ;;; (hexagon-number n) ;;; ;;; * n 番目の六角数 ;;; * Hn = n * (2 * n - 1) ;;; (define hexagon-number (lambda (n) (- (* 2 n n) n))) ;;; ;;; (hexagon-number? n) ;;; ;;; * n が六角数なら「何番目か」を返す。 ;;; (define hexagon-number? (lambda (n) (let* ([a (sqrt (+ 1 (* 8 n)))] [b (/ (+ 1 a) 4)]) (if [positive-integer? b] b #f)))) ;;; ;;; (heptagon-number n) ;;; ;;; * n 番目の七角数 ;;; (define heptagon-number (lambda (n) (/ (- (* 5 n n) (* 3 n)) 2))) ;;; ;;; (heptagon-number? n) ;;; ;;; * n が七角数なら「何番目か」を返す。 ;;; (define heptagon-number? (lambda (n) (let* ([a (sqrt (+ 9 (* 40 n)))] [b (/ (+ 3 a) 10)]) (if [positive-integer? b] b #f)))) ;;; ;;; (octagon-number n) ;;; ;;; * n 番目の八画数 ;;; (define octagon-number (lambda (n) (- (* 3 n n) (* 2 n)))) ;;; ;;; (octagon-number? n) ;;; ;;; * n が八角数なら「何番目か」を返す。 (define octagon-number? (lambda (n) (let* ([a (sqrt (+ 4 (* 12 n)))] [b (/ (+ 2 a) 6)]) (if [positive-integer? b] b #f)))) ;;; ;;; (sum-of-divisors n) ;;; ;;; * 「真の約数の和」を求める。(約数に n 自身は含まれない) ;;; * n < 1 の場合エラーになる。 ;;; * 詳しくは『数学ガール』参照。 ;;; (define sum-of-divisors (lambda (n) (let ([lst (map (lambda (ls) (let ([f (car ls)] [i (cdr ls)]) (/ (- (expt f (+ i 1)) 1) (- f 1)))) (factorize n))]) (- (apply * lst) n)))) ;;; ;;; (perfect-number? n) ;;; ;;; * n が「完全数」なら #t ;;; (define perfect-number? (lambda (n) (= n (sum-of-divisors n)))) ;;; ;;; (deficient-number? n) ;;; ;;; * n が「不足数」なら #t ;;; (define deficient-number? (lambda (n) (> n (sum-of-divisors n)))) ;;; ;;; (abundant-number? n) ;;; ;;; * n が「過剰数」なら #t ;;; (define abundant-number? (lambda (n) (< n (sum-of-divisors n)))) ;;; ;;; (amicable-pair n) ;;; ;;; * n が「友愛数」の場合、ペアとなる数を返す。 ;;; * n が「友愛数でない場合は、#f を返す。 ;;; (define amicable-pair (lambda (n) (let ([sum (sum-of-divisors n)]) (cond ([= n sum] #f) ([= n (sum-of-divisors sum)] sum) (else #f))))) ;;; ;;; (pythagorean-triples n) ;;; (pythagorean-triples n flag) ;;; ;;; * a² + b² = c² かつ a + b + c = n (a < b < c) を満たす ;;; a, b, c の組み合わせを返す。 ;;; * flag が #f 以外の場合は、既約であるものだけを返す。 ;;; ex : (pythagoras-triples 120) ;;; => ((20 48 52) (24 45 51) (30 40 50)) ;;; ex : (pythagoras-triples 120 #t) => () ;;; (define pythagorean-triples (case-lambda ([n] (pythagorean-triples n #f)) ([n flag] (define check-b (lambda (a) (let ([b (- n (/ (* n n) (- n a) 2))]) (if [and (integer? b) (< a b)] (list a b (- n a b)) #f)))) (if [odd? n] '() (let* ([lst-a (sequence 1 (real->integer (/ n 3)))] [lst-b (remq #f (map check-b lst-a))]) (if flag (filter (lambda (lst) (= 1 (apply gcd lst))) lst-b) lst-b))))))

;;;; ;;;; euler.ss --- Project Euler のための関数群 ;;;; ;;;; for Chez Scheme by Tsumuji ;;;; ;;; ;;; (integer->list n) ;;; (integer->list n radix) ;;; ;;; * 整数の各桁を要素とするリストを返す。 ;;; * radix を指定すると、radix を基数にしたリストを返す。 ;;; ex : (integer->list 123) => (1 2 3) ;;; ex : (integer->list 123 2) => (1 1 1 1 0 1 1) ;;; (define integer->list (case-lambda ([n] (integer->list n 10)) ([n radix] (do ([n n (quotient n radix)] [result '() (cons (remainder n radix) result)]) ([< n radix] (cons n result)))))) ;;; ;;; (list->integer lst) ;;; (list->integer lst radix) ;;; ;;; * リストの要素を各桁とする整数を返す。 ;;; * radix を指定すると、radix を基数とした整数を返す。 ;;; ex : (list->integer '(1 2 3)) => 123 ;;; ex : (list->integer '(1 1 1 1 0 1 1) 2) => 123 ;;; (define list->integer (case-lambda ([lst] (list->integer lst 10)) ([lst radix] (fold-left (lambda (a b) (+ b (* radix a))) 0 lst)))) ;;; ;;; (sum-of-digits n) ;;; ;;; * 整数の各桁の合計を求める。 ;;; ex : (sum-of-digits 123) => 6 ;;; (define sum-of-digits (lambda (n) (apply + (integer->list n)))) ;;; ;;; (unique lst) ;;; ;;; * lst の重複する要素を取り除く。 ;;; ex : (unique '(a b a b c)) => (a b c) ;;; (define unique (lambda (lst) (if [null? lst] '() (let ([a (car lst)]) (cons a (unique (remove a (cdr lst)))))))) ;;; ;;; (pandigital? obj) ;;; ;;; * obj が 「Pandigital」ならば #t ;;; * obj は正の整数またはリスト。 ;;; (define pandigital? (lambda (obj) (cond ([positive-integer? obj] (pandigital? (integer->list obj))) ([list? obj] (let* ([lst-a (list-head '(1 2 3 4 5 6 7 8 9) (length obj))] [lst-b (sort < obj)]) (equal? lst-a lst-b))) (else #f)))) ;;; ;;; (all-different? lst) ;;; ;;; * lst の要素がすべて異なれば #t ;;; (define all-different? (lambda (lst) (cond ([null? lst] #t) ([member (car lst) (cdr lst)] #f) (else (all-different? (cdr lst)))))) ;;; ;;; (palindromic? obj) ;;; ;;; * obj が「回文」なら #t ;;; * obj は正の整数またはリスト。 ;;; (define palindromic? (lambda (obj) (cond ([positive-integer? obj] (palindromic? (integer->list obj))) ([list? obj] (equal? obj (reverse obj))) (else #f)))) ;;; ;;; (make-palindromic-number n) ;;; (make-palindromic-number n #t) ;;; ;;; * 第 1 引数 n を上位の桁とする回文数を作る。 ;;; * 第 2 引数がない場合は奇数桁の回文数を作る。 ;;; * 第 2 引数が #t の場合、偶数桁の回文数を作る。 ;;; ex : (make-palindromic-number 123) => 12321 ;;; ex : (make-palindromic-number 123 #t) => 123321 ;;; (define make-palindromic-number (case-lambda ([n] (make-palindromic-number n #f)) ([n flag] (let* ([a (integer->list n)] [b (reverse a)]) (list->integer (append a (if flag b (cdr b)))))))) ;;; ;;; (partial-list lst start length) ;;; ;;; * lst の start の位置から、長さ length の部分リストを返す。 ;;; (define partial-list (lambda (lst start length) (list-head (list-tail lst start) length)))

2009年10月18日 (日)

九九表を作る。

プラプラ〜と Scheme 関連のブログを眺めていたら、"九九表を作る"というのを見つけました。

かなり前の記事なんだけど、"名前付き let" を使ったループを回して頑張っているようです。

私も "名前付き let" は大好きで、良く使うけど、この問題の場合はそんなことしなくても、かなり単純化できることに気づきました。

こんな感じ……。

(define kuku (lambda () (let ([lst '(1 2 3 4 5 6 7 8 9)]) (map (lambda (x) (map (lambda (y) (* x y)) lst)) lst))))
> (kuku) ((1 2 3 4 5 6 7 8 9) (2 4 6 8 10 12 14 16 18) (3 6 9 12 15 18 21 24 27) (4 8 12 16 20 24 28 32 36) (5 10 15 20 25 30 35 40 45) (6 12 18 24 30 36 42 48 54) (7 14 21 28 35 42 49 56 63) (8 16 24 32 40 48 56 64 72) (9 18 27 36 45 54 63 72 81))

2009年6月22日 (月)

Project Euler - Problem 14

問題はこちらをご覧ください。


有名な「Collatz の問題」です。

以前、この「Collatz の問題」を知ったときに、大変興味深かったので、こんなスクリプトを書いて、いろんな数がどうやって 1 に収束していくかを実験したことがありました。

;; Chez Scheme 版 (define collatz (lambda (n) (let loop ((s 1) (i n)) (printf "step ~3d : ~8d~%" s i) (cond ((= i 1) (printf "終了~%")) ((even? i) (loop (+ s 1) (/ i 2))) (else (loop (+ s 1) (+ (* i 3) 1)))))))

今回は生成される数列の長さだけが必要なので、まずこんなスクリプトを書いてみました。

MAX = 100_0000 max = [1, 1] 2.upto(MAX - 1) do |n| v = n c = 1 begin if v.even? then v = v / 2 else v = v * 3 + 1 end c = c + 1 end until v == 1 max = [n, c] if c > max[1] end puts max[0]

このスクリプトは真面目にすべての数を計算するので、非常に遅いです。

Ruby 1.9 で 34 秒ほどかかりました。


ところで、この「Collatz の問題」ではすべての数が 1 に収束すると仮定されています。この仮定が正しいとするならば、すべての数はどこか途中から同じプロセスをたどって 1 に収束していくはずです。

そこで、小さい数の順に生成される数列の長さを記録していき、数列の中に既に記録された数が現れた場合はそれ以上計算しないようにしてみました。

いわゆる Memoize というやつです。『まつもとゆきひろ コードの世界』で「空間で時間を買う」と表現された高速化の一種です。

MAX = 100_0000 collatz = Array.new(MAX, 1) max = [1, 1] 2.upto(MAX - 1) do |n| v = n c = 0 while v >= n if v.even? then v = v / 2 else v = v * 3 + 1 end c = c + 1 end count = c + collatz[v] collatz[n] = count max = [n, count] if count > max[1] end puts max[0]

これだと、Ruby 1.9 で約 2 秒ほどで答が出ました。効果は絶大です。

しかし N88-BASIC の時代からコンピュータに触れている私にとって、「空間で時間を買う」という表現はとても贅沢な響きを持っています。

なんたって、私が初めて買った「PC-8801 MH」のメモリは 128KB でしたからねぇ……。ところが今や、1GB や 2GB は当たり前……。

時代の流れを感じますねぇ……。

2009年5月29日 (金)

素因数分解

あるブログで「素因数分解の速度」を競っているような話題が出ていました。

二人が仲がよくてただじゃれあっているのか、真剣に競い合っているのかは、私の知るところではないのですが、「素因数分解の速度」については思うところがあるので、私も記事を書いてみようかな……。

もし、二人が「素数列をつくる速度」を競っていて、その延長線上で「因数分解の速度」の話になったのなら的外れになるかもしれませんが……。


「素因数分解」の速度を上げようと思ったら、素数を探して割っていくより、対象の数を 2 で割り続けて奇数にした後に、単純に小さい順に奇数で割っていくほうが速いと思います。

なにせ、素数を求めること自体にけっこう時間がかかるのですから。

(さらに、5 以上の素数 p は必ず p = 6 * n - 1 または p = 6 * n + 1 と表せる事を利用するとさらに割る数が少なくなります)


と、いうことで奇数で割り続ける「素因数分解」の Ruby 版です。

# = factorize.rb class Integer def factorize return [] if self < 2 ans = Array.new n = self [2, 3].each do |i| c = 0 while n % i == 0 c = c + 1 n = n / i end ans.push([i, c]) if c > 0 end i = 5 add = 2 while true c = 0 while n % i == 0 c = c + 1 n = n / i end ans.push([i, c]) if c > 0 i = i + add break if n < i * i add = 6 - add end ans.push([n, 1]) if n > 1 return ans end end n = ARGV[0].to_i print "#{n.factorize}\n"
$ ruby -v ruby 1.9.1p129 (2009-05-12 revision 23412) [i686-linux] $ time ruby factorize.rb 111111111 [[3, 2], [37, 1], [333667, 1]] real 0m0.078s user 0m0.020s sys 0m0.012s

Scheme 版はこちら。(普段は "Chez Scheme" を使っているのですが、同じ条件で比較するために "Gauche" で動くようにしてみました)

(define factorize (lambda (n) (define result '()) (define factor-count (lambda (n i) (let loop ([n n] [c 0]) (cond ([zero? (remainder n i)] (loop (/ n i) (+ c 1))) (else (if [> c 0] (set! result (cons (cons i c) result))) n))))) (set! n (factor-count n 2)) (set! n (factor-count n 3)) (let loop ([n n] [i 5] [add 2]) (cond ([= n 1] (reverse result)) ([> (* i i) n] (reverse (cons (cons n 1) result))) (else (loop (factor-count n i) (+ i add) (- 6 add))))))) (define (main args) (print (factorize (string->number (cadr args)))))
$ gosh -V Gauche scheme interpreter, version 0.8.13 [utf-8,pthreads] $ time gosh factorize.scm 111111111 ((3 . 2) (37 . 1) (333667 . 1)) real 0m0.060s user 0m0.012s sys 0m0.016s

では、Gemma さんのコードを見てみましょう。

(ちなみにここに書いてあるコードは終了判定に問題があったようで、こちらに改訂版が書いてあります。以下の記事では改訂版をもとに話をします)

$ gosh -V Gauche scheme interpreter, version 0.8.13 [utf-8,pthreads] $ time gosh factStrm.scm 111111111 ((3 . 2) (37 . 1) (333667 . 1)) real 0m0.174s user 0m0.092s sys 0m0.008s

このコードのスピードは Ruby 版よりは遅いですが、まだまだ実用の範囲にあると思います。


ところで、こちらのコードはというと……

まず "primes.py" というファイルが無いとエラーが出ます。しかも "primes.py" には 997 までしか素数が書き込まれないので、理論的には 994009 以上の数は素因数分解できません。

実際に実効してみると…… 7884 以上では "RuntimeError: maximum recursion depth exceeded in cmp" が出てまともな答えが出ません。

そこでこのコードが使える範囲で速度を見てみると……

$ python -V Python 2.6.2 $ time python fact.py 7882 [2, 7, 563] real 0m0.383s user 0m0.332s sys 0m0.004s $ time python fact.py 7882 [2, 7, 563] real 0m0.326s user 0m0.252s sys 0m0.012s

"primes.py" が空の一回目も "primes.py" が書き込まれた後の二回目もあまりスピードが変わりません。

同じ数字を私の Ruby 版で試してみると……

$ time ruby factorize.rb 7882 [[2, 1], [7, 1], [563, 1]] real 0m0.062s user 0m0.024s sys 0m0.004s

Python 2.6 と Ruby 1.9 では処理系自体のスピードにも差があるのかもしれませんが、Python 版は Ruby 版の約 6 倍時間がかかっています。

使用できる範囲も狭いし、速度も遅い……少なくとも Project Euler クラスの問題を解くにはとても使えません。


※ Gemma さんやみずぴーさんのところのトラックバックをたどって来られる方がいらっしゃるようなので、新しく分かった事実などを元に、以前書いた記事を書き直しました。(2009/09/02)

※ 記事を書き直した際に、Scheme 版を書き忘れていたようなので、改めて最新版を載せました。(2009/10/19)

2009年5月22日 (金)

Project Euler - Problem 35

* 今回の内容は、以前投稿した内容を一部修正して、Ruby のコードを加えたものです。


Problem 35 は次のような問題でした。

197は巡回素数と呼ばれる. 桁を回転させたときに得られる数 197, 971, 719 が全て素数だからである.

100未満には巡回素数が13個ある: 2, 3, 5, 7, 11, 13, 17, 31, 37, 71, 73, 79, および 97 である. 100万未満の巡回素数は何個か?

100 万未満の素数の全てで巡回素数かどうかを調べると、

(length (prime-list 1000000)) ; => 78498

78498 個もの素数を調べる必要があります。


そこで、もう少し考えてみましょう。

2 桁以上の素数で、偶数または 5 を含むものは、巡回しているうちに偶数または 5 が一の位になった時点で、素数でなくなります。ならば偶数や 5 を含む素数は最初から除外しておけば、調べる数が減るはずです。

ということで、私は以下のようなコードを書いてみました。いかがでしょうか?

(define problem-035 (lambda (n-max) ;; 偶数または 5 を含んでいるか? (define include_even_or_5? (lambda (lst) (if [null? (cdr lst)] #f (exists (lambda (x) (memv x lst)) '(0 2 4 5 6 8))))) ;; 循環素数となるリストか? (define circular? (lambda (lst) (let* ([len (length lst)] [lst-a (append lst lst)]) (let loop ([lst-a (cdr lst-a)]) (cond ([= (length lst-a) len] #t) ([prime? (list->integer (list-head lst-a len))] (loop (cdr lst-a))) (else #f)))))) (let* ([lst-a (map integer->list (prime-list n-max))] [lst-b (remp include_even_or_5? lst-a)] [lst-c (filter circular? lst-b)]) (printf "~d~%~a~%" (length lst-c) (map list->integer lst-c))))) (problem-035 1000000)


ちなみに、Ruby だとこんなかんじでしょうか?

# -*- coding: utf-8 -*- include Math class Integer # # 素数の判定 # def prime? return false if (self <= 1) return true if (self == 2) return false if (self.even?) 3.step(sqrt(self).to_i, 2) do |i| return false if (self.modulo(i) == 0) end return true end # # 整数を配列に変換する。 # ex : 123.to_a => [1, 2, 3] # ex : 123.to_a(2) => [1, 1, 1, 1, 0, 1, 1] # def to_a(radix = 10) n = self a = Array.new while (n >= radix) do a.unshift(n.modulo(radix)) n = n / radix end a.unshift(n) end end # # n 以下の素数を要素とする配列を返す。 # ex : prime_list(20) => [2, 3, 5, 7, 11, 13, 17, 19] # def prime_list(n) return([]) if (n <= 1) size = n / 2 - 1 max = sqrt(n).to_i a = Array.new(size) size.times do |i| a[i] = 2 * i + 3 end a.unshift(nil) max.times do |i| j = a[i] next unless j (2 * i * (i + 1)).step(size, j) do |k| a[k] = nil end end return a.unshift(2).compact end class Array # # 配列を整数に変換する。 # ex : [1, 2, 3].to_i => 123 # ex : [1, 1, 1, 1, 0, 1, 1].to_i(2) => 123 # def to_i(radix = 10) sum = 0 self.each do |i| sum = sum * radix + i if (i < radix) end return sum end # # 偶数または 5 を含むか? # def check_even_and_5 return false if self.length == 1 if [0, 2, 4, 5, 6, 8].any? {|x| self.include?(x)} then return true else return false end end # # 循環素数になるリストか? # def circular? a_size = self.size a = self * 2 a_size.times do |i| next if i.zero? return false if not a[i, a_size].to_i.prime? end return true end end p_list = prime_list(100_0000).map {|x| x.to_a} p_list = p_list.delete_if {|x| x.check_even_and_5} p_list = p_list.select {|x| x.circular?} print("#{p_list.size}\n") p p_list.map {|x| x.to_i}

2009年5月 6日 (水)

Project Euler - Problem 45

Problem 45 です。

この問題は、ちょっといやらしいです。

n = 2 * m - 1 の時、Tn = Hm になるので、六角数は完全に三角数に含まれてしまいます。つまり、六角数であれば、必ず三角数になるのです。

ということで、この問題を解くのに三角数の判定は必要ありません。


実際のコードは 144 番目以降の六角数を探していって、五角数かどうかを判定しています。

(define pentagon-number? (lambda (n) (let* ([a (sqrt (+ 1 (* 24 n)))] [b (/ (+ 1 a) 6)]) (if [positive-integer? b] b #f)))) (define hexagon-number (lambda (n) (- (* 2 n n) n))) (define problem-045 (lambda (n) (do ([i n (+ i 1)]) ([pentagon-number? (hexagon-number i)] (printf "~d~%" (hexagon-number i)))))) (problem-045 144)

ruby で書くとこんな感じでしょうか?

class Integer # = 五角数であれば、何番目かを返す。 def pentagon? x = (1 + Math.sqrt(1 + 24 * self)) / 6 if (x == x.to_i) then return x.to_i else return false end end # = 六角数を求める。 def hexagon return (2 * self * self - self) end end 144.upto(1.0/0.0) do |i| h = i.hexagon if (h.pentagon?) then puts h break end end

2009年4月26日 (日)

Project Euler - Problem 39 : ループの抽象化

問題はこちら

今回は、以前紹介した pythagorean-triples を使って、1000 未満の偶数を総当たりしていますが、できるだけ Scheme らしいコードを書くつもりで、for-each を使ってループを抽象化してみました。(でも、set! を使っているのはちょっと無粋かな?)

さらに、pythagorean-triples 自体もループの抽象化を試みてみました。

(define pythagorean-triples (case-lambda ([n] (pythagorean-triples n #f)) ([n flag] (define check-b (lambda (a) (let ([b (- n (/ (* n n) (- n a) 2))]) (if [and (integer? b) (< a b)] (list a b (- n a b)) #f)))) (if [odd? n] '() (let* ([lst-a (sequence 1 (real->integer (/ n 3)))] [lst-b (remq #f (map check-b lst-a))]) (if flag (filter (lambda (lst) (= 1 (apply gcd lst))) lst-b) lst-b)))))) (define problem-039 (lambda (n-max) (let ([v-max '(0 . 0)]) (for-each (lambda (x) (let ([len (length (pythagorean-triples x))]) (if [> len (cdr v-max)] (set! v-max (cons x len))))) (seq 4 n-max 2)) (printf "~d~%" (car v-max))))) (problem-039 999)

P.S. これまでにこのブログで紹介した自作の手続きは、説明抜きで使用しています。R5RS や R6RS に記載されていない手続きに関しては、過去のブログを探してみてください。

2016年7月
          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
31            
フォト

最近のトラックバック

無料ブログはココログ