« 2009年3月 | トップページ | 2009年5月 »

2009年4月

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 に記載されていない手続きに関しては、過去のブログを探してみてください。

2009年4月23日 (木)

Project Euler - Problem 44

問題はこちらを参照してください。


特定の範囲の五角数を列挙して、その中で条件に合うものを探していく方法もありますが、それだと、

  • 範囲の上限を決めるのが難しい。
  • 差 D が最小であることを確認するのが難しい。

という問題があると思います

そこで今回は、次のような作戦を考えました。

  1. P1 = 1 を list に push しておく。
  2. 小さい順に五角数 (Pn) を求める。
  3. Pn と list の中の数を list の先頭から比較していく。
  4. 条件に合うペアが見つかれば終了。見つからなければ Pn を list の先頭に push して 2. に戻る。

この方法だと、上限を設定せずにペアが見つかるまで探索を続けられます。また、list の中の数は大きい順に並んでいるので、条件に合うペアが見つかった時点で、差 D は最小になっているはずです。


まずは Scheme 版です。

;;; ;;; (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 [integer? b] b #f)))) ;;; Problem 44 (define problem-044 (lambda () (define pentagon-lst '(1)) (define check (lambda (p) (exists (lambda (x) (if [and (pentagon-number? (- p x)) (pentagon-number? (+ p x))] (cons p x) #f)) pentagon-lst))) (let loop ([i 2]) (let* ([p (pentagon-number i)] [ans (check p)]) (if ans (let ([a (car ans)] [b (cdr ans)]) (printf "~d : (~d, ~d)~%" (- a b) a b)) (begin (set! pentagon-lst(cons p pentagon-lst)) (loop (+ i 1)))))))) (problem-044)

ちょっと時間はかかりますが、約 10 秒で答が出ます。


次に Ruby 版です。

class Integer def pentagon return ((3 * self * self - self) / 2) end def pentagon? x = (1 + Math.sqrt(1 + 24 * self)) / 6 return (x == x.to_i) end end penta_list = [1] 2.upto(1/0.0) do |i| p = i.pentagon penta_list.each do |j| if ((p - j).pentagon? and (p + j).pentagon?) then print"#{p - j}\n" exit end end penta_list.unshift(p) end

今回はちょっと驚くことがありました。

これまで同じアルゴリズムでコードを書いた場合、 Chez Scheme のほうが Ruby1.8 より数倍速く結果が出たのですが、なぜか今回はほぼ同じ時間で結果が出ました。Ruby1.9 を使った場合は、約 5 秒で答がでました。

原因ははっきりしませんが、Chez Scheme の list 処理よりも Ruby の配列の処理のほうが効率が良いのでしょうか?

2009年4月18日 (土)

Project Euler - Problem 4 再び

以前、「Problem 4」に関する記事を書いた際に、「二重ループで探すのは時間がかかる」と書いたのですが、それは私の思い違いであることが分かりました。

正確に書くと、「工夫次第でかなり速くなる」のが分かりました。

前回投稿した後、いろいろなブログを見ていて、次のようなコード見つけました(記憶を元に私が再構成したので、元のコードとは違う部分があると思いますが、アルゴリズムは同じはずです)。

class Integer def palindrome? s = self.to_s return s == s.reverse end end max = 0 999.downto(100) do |i| 999.downto(100) do |j| num = i * j break if (num <= max) max = num if (num.palindrome? and (num > max)) end end puts max

このコードでポイントとなるのは、大きい数から探していくことと、12 行目の "break" です。

その時点で分かっている最大値より小さい値の計算を飛ばすことによって、全体の計算量を減らしています。

調べてみたら、"num = i * j" の計算は 9200 回しか行われていませんでした。

アルゴリズムも、コードも自分が考えていたものよりもシンプルでわかりやすいと思います。

「ちょっとした工夫」を考えつけなかったのが非常に悔しいです……。

2009年4月17日 (金)

Project Euler - Problem 1 再び

以前解いた問題も、時間が経ってから見直すと、また違った方法を考えつくことがあります。

というわけで、再び "Problem 1" です。

以前投稿したコードは次のようなものでした。

(define problem-001 (lambda (n) (printf "~d~%" (apply + (filter (lambda (n) (or (zero? (remainder n 3)) (zero? (remainder n 5)))) (iota n)))))) (problem-001 1000)

リストのデータを加工していくという点で、Scheme らしいかな?と思っています。

以前、考えていたもう一つのコードは反復的プロセスを使うパターンでした。

(define problem-001 (lambda (n-max) (let loop ([i 1] [sum 0]) (cond ([= i n-max] (printf "~d~%" sum)) ([or (zero? (remainder i 3)) (zero? (remainder i 5))] (loop (+ i 1) (+ sum i))) (else (loop (+ i 1) sum)))))) (problem-001 1000)

これは手続き言語っぽいやり方ですね。


今回紹介するアルゴリズムは、Ruby でコードを書いていたときに考えついたものです。

ということで、まずは元になった Ruby のコードから……

MAX = 999 a = Array.new(MAX + 1, 0) step = [3, 5] 2.times do |i| 0.step(MAX, step[i]) do |j| a[j] = j end end puts(a.reduce(:+))

自作の「エラトステネスの篩」の応用で、割り算を使わないので結構速いです。

これを Scheme で書くと……

(define problem-001 (lambda (n-max) (define vct (make-vector n-max 0)) (define value-set! (lambda (step) (do ([i 0 (+ i step)]) ([>= i n-max]) (vector-set! vct i i)))) (for-each value-set! '(3 5)) (printf "~d~%" (apply + (vector->list vct))))) (problem-001 1000)

Scheme で代入の "vector-set!" を使うのはあまり好きではないのですが、こういうやり方もあるんだなぁ……ということで。

2009年4月15日 (水)

Project Euler - Problem 43 : Pandigital 数を作ってしまえ

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

数 1406357289 は 0 から 9 の Pandigital 数である ( 0 から 9 が 1 度ずつ現れるので).この数は部分語が面白い性質を持っている.

d1 を 1 桁目, d2 を 2 桁目の数とし, 以下順に dn を定義する. この記法を用いると次のことが分かる.

  • d2d3d4 = 406 は 2 で割り切れる
  • d3d4d5 = 063 は 3 で割り切れる
  • d4d5d6 = 635 は 5 で割り切れる
  • d5d6d7 = 357 は 7 で割り切れる
  • d6d7d8 = 572 は 11 で割り切れる
  • d7d8d9 = 728 は 13 で割り切れる
  • d8d9d10 = 289 は 17 で割り切れる

このような性質をもつ 0 から 9 の Pandigital 数の総和を求めよ.

まず考えるのは、10 桁の Pandigital 数を順列を使って求めて、条件に合うものをピックアップしていく方法でしょう。

でも 10 桁の Pandigital 数は、d1 に 0 が来るときも含めると、10P10 = 3628800 通りになります。まともにループを回すのはちょっと大変なので、別な方法を考えてみます。

Problem 37 の考え方を使って、条件に合う数を作ってしまいましょう。(だって、1000 未満の 17 の倍数なんて、58 個しかないんですよ)

ただし、今回は途中に "0" が来てもいいので、データを数値のまま使っていくと、桁落ちをする可能性があります。そこで、今回はデータをリストの形で保持して、必要に応じて数値に変換する方法をとりました。

(なお、「自作の手続きたち」で紹介した手続きを使用していますので、こちらも参照してみてください)

;; * lst1 から lst2 の要素を取り除く (define diff-list (lambda (lst1 lst2) (let loop ([lst1 lst1] [lst2 lst2]) (if [null? lst2] lst1 (loop (remove (car lst2) lst1) (cdr lst2)))))) (define problem-043 (lambda () (define seq-0-9 (iota 10)) ;; d8d9d10 の候補をりすとにまとめて返す。 ;; => ((0 1 7) (0 3 4) ... (9 6 9) (9 8 6)) (define make-first-list (lambda () (filter all-different? (map (lambda (n) (if [< n 100] (cons 0 (integer->list n)) (integer->list n))) (sequence 17 999 17))))) ;; 以下の条件に合う新しいリストを返す。 ;; * すべての数字が異なる ;; * リストの前から 3 個の数字をつなげたものが d で ;; 割り切れる。 (define make-new-list (lambda (lst d) (filter (lambda (ls) (zero? (remainder (list->integer (list-head ls 3)) d))) (map (lambda (n) (cons n lst)) (diff-list seq-0-9 lst))))) ;; d1 を決めて数字を完成させる。 (define finish (lambda (lst) (let ([rest (- 45 (apply + lst))]) (if [zero? rest] 0 (list->integer (cons rest lst)))))) (let loop ([lst (make-first-list)] [d-lst '(13 11 7 5 3 2)]) (if [null? d-lst] (printf "~d~%" (apply + (map finish lst))) (loop (apply append (map (lambda (ls) (make-new-list ls (car d-lst))) lst)) (cdr d-lst)))))) (problem-043)

自宅のパソコンで実行したら、REPL 上では約 0.002 秒、Shell 上でスクリプトとして実行した場合には約 0.16 秒で答えが出ました。


ちなみに、その後、順列を使って 10 桁の Pandigital 数を作っていく方法も試して見ましたが、7 秒で答えが出ました。案外速くてビックリでした。


※追記 : Ruby でコードを組んだときの手法をフィードバックしてコードを若干改良しました。(2009/08/04)

2009年4月13日 (月)

Project Euler - Problem 37 : 素数を作ってしまえ

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

3797は面白い性質を持っている. まずそれ自身が素数であり, 左から右に桁を 除いたときに全て素数になっている (3797, 797, 97, 7). 同様に右から左に桁を除いたときも全て素数である (3797, 379, 37, 3).

右から切り詰めても左から切り詰めても素数になるような素数は11個しかない. 総和を求めよ.

注: 2, 3, 5, 7を切り詰め可能な素数とは考えない.

まず初めに考えたのが、とにかく条件に合う素数を 11 個探すという方法。

ただし探索する素数の上限が分からないので、奇数を順番に数え上げていって、その数が素数かどうかを判定してから、切り詰め可能かを調べています。

そのため、結果が出るまでかなり時間がかかります (実際、7.6秒かかりました)。

(define problem-037 (lambda () ;; 左右から切り詰めた時、全ての数が素数になるか? (define check (lambda (p) (let loop ([d 10]) (let ([a (quotient p d)] [b (remainder p d)]) (cond ([zero? a] #t) ([and (prime? a) (prime? b)] (loop (* d 10))) (else #f)))))) (let ([result '()] [count 0]) (let loop ([i 11]) (cond ([= count 11] (printf "~d : ~a~%" (apply + result) (reverse result))) ([and (prime? i) (check i)] (set! result (cons i result)) (set! count (+ count 1)) (loop (+ i 2))) (else (loop (+ i 2)))))))) (problem-037)

このコードでは、条件に合う素数が 11 個しかないことを前提にして問題を解いていますが、そもそも 11 個しかないことはどうやって確認されたのでしょうか?

そう考えると、他にもアプローチの仕方があるはずです。


ということで、考え方を逆転させて、条件に合う素数を作ってみることにします。どういうことかというと、一桁の素数の右か左に数字を一つずつ付けていって素数になるものを探してみます。

では、右と左、どちらに付けていくほうがいいのでしょうか?


まず、左に付けていく場合。

実は一桁目が素数の場合、左に付ける数が 1 〜 9 のいずれでも素数になる可能性があります。


次に、右に付けていく場合。

二桁以上の数では、一桁目が偶数や 5 の場合素数ではなくなるので、右に付けられる数は 1, 3, 7, 9 に限定されます。

これだけ限定されれば、どんどん右に数を付け足していっても、いつかは素数が作れなくなるはずです

そこで、素数が作れなくなるまで右に数を付け足していって、できた素数が左詰め可能かをチェックしていきます。

(define problem-037 (lambda () ;; 引数の下に数字を一つ加えて素数を作る。 ;; 右から切り詰める場合の逆の動作。 (define make-new-prime (lambda (p) (filter prime? (map (lambda (x) (+ (* p 10) x)) '(1 3 7 9))))) ;; 引数を左から切り詰める場合の逆の動作でチェック。 (define check (lambda (n) (let loop ([d 10]) (let ([m (remainder n d)]) (cond ([= m n] #t) ([prime? m] (loop (* d 10))) (else #f)))))) (let loop ([lst '(2 3 5 7)] [ans '()]) (if [null? lst] (let ([result (filter check ans)]) (printf "~d : ~a~%" (apply + result) (reverse result))) (let ([new-lst (apply append (map make-new-prime lst))]) (loop new-lst (append new-lst ans))))))) (problem-037)

このコードの場合、探索範囲がかなり限定されるので、0.04 秒しかかかりませんでした。

2009年4月 5日 (日)

『数学ガール』が面白い

つい先日から読み始めた『数学ガール』が非常に面白い。
新しい発見、好奇心をそそる話題など、もっと早く読んどくんだった……。

例えば、自分が Scheme でプログラミングする際に作った階乗の特殊形が、「下降階乗冪」という名前で出てきたのを見て、「自分以外にも、組合せの計算で律儀に階乗を使って計算するのが面倒だと考えた人がいたのかな?」と思うと、ちょっとうれしくなったり……。

『オイラーの贈り物』で二項定理の説明を読んだときには、「そんなものなんだ……」という具合に流してしまったんだけど、『数学ガール』の中で、「n 個の x の中から何個の x をとるか」という説明をしてあるのを読んで、「これは確かに組合せになるなぁ……」とすっきりと納得ができたり……。

今のところ、難しい計算は飛ばしてとりあえず物語を読んで行く方針だけど、その次には本の中に出てきた計算を自分で確かめながらじっくりと読んでみようと思っている。

2009年4月 3日 (金)

数学ガール --- ωのワルツ

Omega1_5


遅ればせながら、昨日から『数学ガール』を読んでいます。

何となく買いそびれていたのですが、昨日、意を決して買ってみました。

いや〜おもしろい。しかも読みやすい。


ところでこの本の中に「ωのワルツ」という章がありますよね。

初項 C0 = 0, 公比 ω = ( -1 + √3 i ) / 2

の等比数列が、1, ω, ω² を繰り返すというやつです。

これが複素平面上で正三角形を形作るという話ですが、本当にそうなのか確かめたくて、 自作の「Turtle Graphics on Ruby/SDL」を使った、下記のようなスクリプトを作ってみました。

# # ωのワルツ # # from 『数学ガール』 # require 'turtle_graphics' require 'complex' class Turtle def move_to(x, y) dx = x.to_f - @x dy = y.to_f - @y case when (dx < 0) @angle = atan(dy / dx) / TO_RADIAN - 180 when (dx == 0) case when (dy < 0) @angle = 270 when (dy > 0) @angle = 90 end when (dx > 0) @angle = atan(dy / dx) / TO_RADIAN end len = sqrt(dx * dx + dy * dy) forward(len) end self end s = Screen.create s.set_zoom(160).grid_on s.set_speed(4) s.draw_circle(0, 0, 1, Blue) t = Turtle.new(1, 0) # Turtle.verbose_on t.set_color(Yellow) omega = Complex.new(-0.5, 0.5 * sqrt(3)) c = 1 loop do c = c * omega x = c.real y = c.image t.move_to(x, y) end s.main_loop

これは、横軸を実数軸、縦軸を虚数軸に見立てた複素平面で、ωの作る等比数列の示す点を 亀がたどって行くものです。青い円は複素平面の単位円です。

単に正三角形を描いているように見えますが、内部では虚数の計算をしながら、実際にωつくる 等比数列の項を順々に求めています。

こうして実際に亀が動く軌跡を見てみると、ちょっとした感動があります。


ちなみに

omega = Complex.new(0, 1)

とすると、『数学ガール』のp44の正方形上の移動を再現できます。

« 2009年3月 | トップページ | 2009年5月 »

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            
フォト

最近のトラックバック

無料ブログはココログ