Hatena::ブログ(Diary)

思いついたら書く日記 Ex このページをアンテナに追加 RSSフィード

2009-12-07

Project Euler Problem 31

| 22:40

問題31: イギリス硬貨 1ペンス〜2ポンドを使い 2ポンドをつくる方法は何通りあるか。[=>問題文]

問題31の解答

;; Emacs lisp
(require 'cl)

(defun problem031a (n)
  (labels
      ((repeat (coins acc)
               (if (null coins)
                   0
                 (do ((add (car coins))
                      (next-coins (cdr coins))
                      (count 0))
                     ((>= acc n) (if (= acc n) (1+ count) count))
                   (incf count (repeat next-coins acc))
                   (incf acc add)))))
    (repeat '(200 100 50 20 10 5 2 1) 0)))

(problem031a 200)

 組み合わせパターンの生成は簡単でしたが、count の置き方、受け渡し方に迷いました。do の終了処理周りがいまいちスッキリしませんが・・・とりあえず答えは出たのでこれで・・・。

 ところでいろいろコードをいじってて思ったんですが、なんか labels って遅くないですか? というわけで実験。

 ロジックは全部同じです

(require 'cl)

;; flet 版
(defun problem031b (n)
  (flet
      ((repeat (coins acc)
               (if (null coins)
                   0
                 (do ((add (car coins))
                      (next-coins (cdr coins))
                      (count 0))
                     ((>= acc n) (if (= acc n) (1+ count) count))
                   (incf count (repeat next-coins acc))
                   (incf acc add)))))
    (repeat '(200 100 50 20 10 5 2 1) 0)))

;; let + lambda 版
(defun problem031c (n)
  (let ((repeat
         (lambda (coins acc)
           (if (null coins)
               0
             (do ((add (car coins))
                  (next-coins (cdr coins))
                  (count 0))
                 ((>= acc n) (if (= acc n) (1+ count) count))
               (incf count (funcall repeat next-coins acc))
               (incf acc add))))))
    (funcall repeat '(200 100 50 20 10 5 2 1) 0)))

 例によって Core2 Quad 6600 にて計測。

Function      Avg (sec)
===========   ==========
problem031a   3.822200  labels
problem031b   0.959400  flet
problem031c   1.037200  let + labmda

 やっぱり labels は遅かった。labels と flet はどちらも cl パッケージのマクロです。レキシカルスコープの実装がオーバーヘッドになっているんでしょうか?

コンパイル時の警告

 let とflet はコンパイル時(compile-defun)に警告が出ます。

;; flet
Warning: repeat called with 2 arguments, but accepts only 1
Warning: repeat called with 2 arguments, but accepts only 1

;; let + lambda
Warning: reference to free variable `repeat'

 後者は、let + lambda で再帰を書くといつもでる警告ですが、前者はどういうことなんでしょう? 引数が 1 つしか accept されない? よくわかりません。

追記: repeat は Emacs Lisp で定義済みの関数でした。直前のコマンドを繰り返す機能です。上記 3 つの関数では flet だけが、グローバルな関数を参照するため、既存の repeat と引数が合わないという警告を出しているようです。

2009-12-06

Project Euler Problem 30

| 03:16

問題30: 2以上の自然数について、各桁を5乗した数の和が元の数と一致する数をすべて求め、その和を答えよ。[=>問題文]

 上限が与えられない問題です。左辺の値いくつまで調べればよいのか。注目するのは「桁数」です。

 k 桁の数を考えます。各桁の最大値は 9 なので右辺 a のとりうる最大値 a-max は次のようになります。

[右辺 a 最大値]
a-max(1): 9^5
a-max(2): 9^5 + 9^5
a-max(3): 9^5 + 9^5 + 9^5
:
a-max(k-1): 9^5 + 9^5 + 9^5 ... 9^5
a-max(k)  : 9^5 + 9^5 + 9^5 ... 9^5 + 9^5

 つまり、a の最大値は桁が増える毎に 9^5 ずつ増えていきます。次は左辺 b の最小値 b-min について考えます。

[左辺 b 最小値]
b-min(1): 1 (0)
b-min(2): 10
b-min(3): 100
:
b-min(k-1): 1000...00  (k-1桁)
b-min(k)  : 1000...000 (k桁)

 k=2 の時は2桁の最小値なので10、k=3の bの最小値は 100 です。b の最小値は桁が増えるごとに 10 倍になります。1桁の最小値は 0 ですが便宜的にここでは 1 とします。

 もう一度 b-min と a-max を並べて書いてみます。

k b-min  a-max
1    1    59049 (9^5)
2   10   118098 (9^5 + 9^5)
3  100   177147 (9^5 + 9^5 + 9^5)
4 1000   236196 (9^5 + 9^5 + 9^5 + 9^5)
:

 このように、はじめのうちは b-min(左辺最小)< a-max (右辺最大)ですが、桁が増える毎の値の増加量は b-min が指数的増加量(というか指数関数そのものですが)であるのに対して、a-max は単調増加です。だから、k を増やしていけばいずれどこかで b-min > a-max となる瞬間があるはずです。そして大小が逆転した後はその差は開いていく一方です。この問題は b = a となる数を探すわけですから、b-min > a-max となった後は調べる必要がありません。

 b-min > a-max となった瞬間の a-max が、調べる値の上限となります。

問題30の解答その1

;;Emacs Lisp
(require 'cl)

(defun problem030a (n)
  (let
      ;; 0^n 〜 9^n を格納した配列。
      ;; n乗計算を省略するためにメモ化しておく。
      ((expt-n-vect
        (do ((v (make-vector 10 0))
             (i 1 (1+ i)))
            ((= i 10) v)
          (aset v i (expt i n))))
       ;; 計算対象の上限。
       ;; k 桁での左辺最小値を bk。 右辺最大値を ak とし、
       ;; bk が ak を超えたとき、その時点の ak の値を
       ;; 計算すべき上限値とする。 
       (limit (do ((k 1 (1+ k))        ; 実は不要だけど、一応
                   (a-max (expt 9 n))  ; a(k)の最大値
                   (b-min 1)           ; b(k)の最小値
                   (a-inc (expt 9 n))) ; a増加量
                  ((<= a-max b-min) a-max)
                (incf a-max a-inc)
                (incf b-min (* b-min 10)))))
    (let*
        ;; 桁分割リスト作成関数。 123 => (1 2 3)
        ((split-digit
          (lambda (x &optional lis)
            (if (zerop x)
                lis
              (funcall split-digit (/ x 10)
                       (cons (mod x 10) lis)))))
         ;; 左辺 b から 右辺 a を計算する関数
         (calc-a
          (lambda (b)
            (reduce
             #'(lambda (acc y)
                 (+ acc (aref expt-n-vect y)))
             (funcall split-digit b)
             :initial-value 0))))
      ;; メインルーチン
      ;; 左辺 b (limit 〜 2) に対する 右辺 a を計算し
      ;; b = a となる b をリストに保存し、合計する。
      (do ((lis) (b limit (1- b)))
          ((= b 1) (apply #'+ lis))
        (when (= b (funcall calc-a b))
          (push b lis))))))

(problem030a 5)

 やってることは、最初の説明のとおりです。最後のメインルーチンでは、左辺 b の値を上限 limit から 2 までデクリメントし、すべての b に対する a を計算して b = a となるものを得ています。

 このやり方だとn=5の場合、約47万回の右辺計算をしています。しかもその中には b = 1234,1342,2341... のように右辺 a の値が同じになるものが多数含まれています。この無駄をなんとかしてみましょう。

問題30の解答その2

(require 'cl)

(defun problem030b (n)
  ;; --- 前半は、problem030a と同じ ----------
  (let ((expt-n-vect                    ; 0^n 〜 9^n を格納した配列
         (do ((v (make-vector 10 0))
              (i 1 (1+ i)))
             ((= i 10) v)
           (aset v i (expt i n))))
       ;; 計算対象の上限。
       (limit (do ((k 1 (1+ k))        ; 実は不要だけど、一応
                   (a-max (expt 9 n))  ; a(k)の最大値
                   (b-min 1)           ; b(k)の最小値
                   (a-inc (expt 9 n))) ; a増加量
                  ((<= a-max b-min) a-max)
                (incf a-max a-inc)
                (incf b-min (* b-min 10)))))
    (let
        ;; 桁分割リスト作成関数。 123 => (1 2 3)
        ((split-digit
          (lambda (x &optional lis)
            (if (zerop x)
                lis
              (funcall split-digit (/ x 10)
                       (cons (mod x 10) lis)))))
         ;; --- ここまで、problem030a と同じ ----------
         ;; リスト要素のn乗和計算関数。 (1 2 3) => 1^n+2^n+3^n
         (sum-of-expt-n
          (lambda (lis) (reduce #'(lambda (acc y)
                                    (+ acc (aref expt-n-vect y)))
                                lis
                                :initial-value 0))))
      (let*
          ;; limitの桁数
          ((limit-digits (1+ (truncate (log limit 10))))
           ;; 2〜limit までの数を桁分割し sortしたリストを key として
           ;; 保存した hash tableを作成。hash table の key なので重複は
           ;; 除去される。
           (ht (do ((ht (make-hash-table :test #'equal))
                    (x limit (1- x))
                    (border (expt 10 (1- limit-digits)))
                    (zero-list))
                   ((= x 1) ht)
                 ;; すべてのリストが同じ長さになるよう zero-list を
                 ;; append して調整する。
                 (let ((lis (append zero-list
                                    (funcall split-digit x))))
                   (puthash (sort lis #'<) t ht))
                 ;; x の桁が減ったら、zero-list に 0 を追加
                 (when (= x border)
                   (push 0 zero-list)
                   (setq border (/ border 10)))))
           ;; 題意に該当する数を保存する為のリスト
           (correspond-numbers))
        ;; 作った hash table (ht) のすべてのキーに対して題意のテスト
        ;; を行い合致するものを correspond-numbers に保存する
        (maphash
         #'(lambda (a-list v)
             ;; a-list は右辺の基数となる数値のリスト。
             ;; a-list の要素の n乗和を 左辺 bとし、
             ;; bを桁にばらして sortしたリストを b-list とする.
             (let* ((b (funcall sum-of-expt-n a-list))
                    (b-list (sort (funcall split-digit b) #'<))
                    (b-list-len (length b-list)))
               ;; a-list の桁数はすべて limit-digits なので
               ;; b-list の桁数も 0 埋めしてあわせる
               (if (< b-list-len limit-digits)
                   (setq b-list
                         (append
                          (make-list (- limit-digits b-list-len) 0)
                          b-list)))
               ;; a-list と b-list が一致したら、それが求める数。
               (if (equal a-list b-list)
                   (push b correspond-numbers))))
         ht)
        ;; correspond-numbers には 1 も含まれるので、
        ;; sort して car を除いてから合算。
        (apply #'+ (cdr (sort correspond-numbers #'<))))))))

(problem030b 5)

 1234,1342,2341...のような重複要素を除去してからメインの計算ループに入ろうという発想です。上限 limit までの数をすべて (1 2 3 4) (1 3 4 2) (2 3 4 1) のようなリストに変換し、sort した上で hash table に key として保存しています。ソートすることにより、上記 3 つのような数はすべて同じパターンとなり hash table によって重複要素は除去されます。

 n=5の場合、すべてのリストを ht に登録し重複除去された時点でパターン数は 4542個まで絞り込まれました。右辺 a の計算回数は、problem030a の 100分の1にまで減っています。その代わり前準備がかなり複雑になっています。

プロファイル

Function     Avg (sec)
===========  =========
problem030a  4.156500
problem030b  2.421500

 hash table 登録や リストのソートなど、オーバーヘッドになりそうな処理をしているにもかかわらず、後者の方式のほうが速いようです。ループ回数 100倍の差が効いているのでしょうか?

2009-12-04

Project Euler Problem 29

| 03:11

問題29: 2<=a<=100, 2<=b<=100 において、a^b で作成される数は重複を除いて何個あるか[=>問題文]

問題29の解答 (list + delete-dups)

;; Emacs Lisp
(require 'cl)
(require 'bigint) ;integer-to-bigint, bigint*

(defun problem029a (limit)
  (length
   (delete-dups
    (do ((a 2 (1+ a)) (lis))
        ((> a limit) lis)
      (let ((n (integer-to-bigint a))
            (acc (integer-to-bigint a)))
        (do ((b 2 (1+ b)))
            ((> b limit))
          (setq acc (bigint* acc n))
          (push acc lis)))))))

(problem029a 100)

 特別なことはしていません。作成される数をリストにすべて保存し、delete-dups で重複を除去しています。100の100乗は明らかに多倍長数なので、bigintライブラリを使っています。

 bigintライブラリに新たに integer-to-bigint を追加しました。整数からbigintへ文字列を経由することなくダイレクトに変換できるようになりました。また掛け算関数 bigint* にバグがあったので修正しました。

 ところで、delete-dups は遅いらしいです。対象データが多いときは hash-table を使ったほうが良いらしい。というわけで。

問題29の解答その2 (hash-table + hash-table-count)

(defun problem029b (limit)
  (hash-table-count
   (do ((a 2 (1+ a)) (ht (make-hash-table :test #'equal)))
       ((> a limit) ht)
     (let ((n (integer-to-bigint a))
           (acc (integer-to-bigint a)))
       (do ((b 2 (1+ b)))
           ((> b limit))
         (setq acc (bigint* acc n))
         (puthash acc t ht))))))

(problem029b 100)

 hash-table はデフォルトではキーの比較に eq を使いますが、今回はbigint、つまりリストをキーにするので equal で比較するよう、make-hash-table に指定しています。ちなみに、delete-dups は equal 比較で固定されています。

どれくらい違うか 

Function     Avg (sec)
===========  ========= 
problem029a  1.615900  delete-dups 版
problem029b  0.196600  hash-table 版

 これはまた、想像以上に差がつきました。

問題29の解答その3

 実際に指数計算しなくても、素因数分解結果を保持していれば生成される数のパターンは数えられることに気づきました。これならば多倍長演算不要です。

(require 'cl)
(require 'math-lib) ;factorize

(defun problem029c (limit)
  (let* ((primes (primes-sieve (truncate (sqrt limit))))
	 (parse
	  (lambda (a b)
	    (mapcar #'(lambda (pair)
			(cons (car pair) (* b (cdr pair))))
		    (factorize a primes)))))
    (hash-table-count
     (do ((a 2 (1+ a)) (ht (make-hash-table :test #'equal)))
	 ((> a limit) ht)
       (do ((b 2 (1+ b)))
	   ((> b limit))
	 (puthash (funcall parse a b) t ht))))))

(problem029c 100)

 parse にて (a b) の組から a^b の計算結果を素因数分解したリストを作成します。factorize は素因数分解の結果を 素因数の小さいほうから並べたリストで返しますので sort 不要です。(順序が保証されなければここで sort が必要になります。

2009-12-03

Project Euler Problem 28

| 00:10

問題28:渦巻状に自然数を並べた1001x1001数方陣の対角線上にある数の合計を求める。方陣の詳細はリンク先参照[=>問題文]

問題28の解答

;; Emacs Lisp
(require 'cl)

(defun problem028 (N)
  (do ((n 3.0 (+ n 2)) (result 1))
      ((> n N) result)
    (let ((d (1- n)) (x (* n n)))
      (incf result x)
      (incf result (decf x d))
      (incf result (decf x d))
      (incf result (decf x d)))))

(problem028 1001)

 これは簡単ですね。外周の一辺のマス目の数を n とすると右上隅は n の二乗。その他の隅は n-1 ずつ減じた数になる。これを合わせれば外周上の四隅の合計が求まり、外周を増やしながらそれらを加算していけば、最終的には対角線上の総計が求まります。

2009-12-02

Project Euler Problem 27

| 00:32

問題27: 二次式 n^2 + an + n = p (|a| < 1000, |b| < 1000, p = 素数) において、n に 0から連続する整数を代入したときもっとも大きな n まで式が成り立つ a b の組を求め、その積を答えよ。[=>問題文]

 効率的な素数判定方法を考えていたら2日経ってしまいました。最終的に、処理時間は 1.89秒になりました。これは速いなのか遅いほうなのか、これから他の人の解答を見るのが楽しみです。

問題27の解答

;; Emacs Lisp
(require 'cl)
(require 'math-lib)

(defun problem027 (limit)
  (let* ((b-list (primes-sieve limit))
         (prime-p (make-prime-number-p b-list))
         (pairs (let ((lis) (b))
                  (dolist (b (reverse b-list) lis)
                    (do ((a (1- limit) (1- a)) (end (* -1 limit)))
                        ((= a end))
                      (push (list a b) lis))))))
    (do ((n 1 (1+ n)))
        ((>= 1 (length pairs)) (apply #'* (car pairs)))
      (let ((lis) (n2 (* n n)) (continue t))
        (while (and (consp pairs) continue)
          (let ((pair (pop pairs)))
            (destructuring-bind (a b) pair
              (if (>= n b)
                  (setq continue nil)
                (let ((x (+ n2 (* a n) b)))
                  (if (and (> x 1) (funcall prime-p x))
                      (push pair lis)))))))
        (setq pairs lis)))))

(problem027 1000)

 n^2+an+b=p に n=0 を代入すると b=p となり、bは素数で無ければならないことがわかります。また、n=b を代入すると左辺は b^2+ab+b となりこれは明らかに b で割り切れます。したがって 連続する nにおいて、式が成り立つ n の最大値は高々 b - 1 であることがわかります。

 problem027 では、はじめに a, b の組み合わせをすべて作成 (pairs) しています。但しbは1000未満の素数だけです。メインのループでは n の値を増やしながら、成り立たないものを除去ししてゆき リスト pairs の要素が 1 つになったところで、処理を完了としています。

 なお素数判定関数 prime-p は「素数判定関数作成関数」によって取得しています。

素数判定関数作成関数

(defun make-prime-number-p (primes)
  (lexical-let*
      ((primes (if (> 2 (length primes)) [2 3] (vconcat primes)))
       (capacity (length primes))
       (index (1- capacity))
       (expand-capacity
        (lambda ()
          (setq primes (vconcat primes (make-vector capacity 0)))
          (setq capacity (* 2 capacity))))
       (iterate)
       (iterate-2
        (lambda ()
          (when (= index (1- capacity))
            (funcall expand-capacity))
          (let ((n (aref primes index)) (r 0))
            (while (zerop r)
              (incf n 2)
              (setq r 1)
              (do ((i 0 (1+ i)) (p 0) (end (truncate (sqrt n))))
                  ((or (>= p end) (zerop r)))
                (setq p (aref primes i))
                (setq r (mod n p))))
            (aset primes (incf index) n))))
       (iterate-1
        (lambda ()
          (if (< index (1- capacity))
              (aref primes (incf index))
            (progn
              (setq iterate iterate-2)
              (funcall iterate-2)))))
       (binary-search
        (lambda (x)
          (let ((i 0) (j index) (exists))
            (while (and (<= i j) (not exists))
              (let* ((k (/ (+ i j) 2)) (m (aref primes k)))
                (cond ((< m x) (setq i (1+ k)))
                      ((> m x) (setq j (1- k)))
                      (t (setq exists t)))))
            exists)))
       (prime-number-p
        (lambda (x)
          (let ((last-prime (aref primes index)))
            (cond
             ((< x last-prime) (funcall binary-search x))
             ((> x last-prime)
              (do ((p (funcall iterate) (funcall iterate)))
                  ((>= p x) (= p x))))
             (t t))))))

    (setq iterate iterate-1)
    (lambda (x) (if (evenp x) 
                    nil
                  (funcall prime-number-p x)))))

 最大値のわからない素数判定に使える関数が欲しくて考えた末できたのがこれです。はじめに 素数を格納した vector primes を作成し、これに収まる整数 x についてはバイナリ探索で該当する素数があるかを判定します。primes よりもおおきな素数が渡されたときは、primes を拡張して対応できるようにしています。一度拡張するとそれよりさらに大きい数がくるまではバイナリ探索ですむので、大量のランダムな数を素数判定するには有利なんじゃないかと。(でも条件によるか・・・小さな素因数を持つ合成数の判定はむしろ遅くなるな、すごくダメな気がしてきた。普通に試し割りの方式に修正するかも)。

 無駄に複雑に見えるかもしれませんが、実はもう少し発展させて便利なものをつくろうと画策しています。うまくいくかな。

ツムジツムジ 2009/12/03 08:38 試し割方式なら、2 と √n 以下の奇数があれば n 以下の整数の素数判定ができるので、Project Euler の問題程度なら、素数判定は試し割方式で十分ではないのでしょうか?

ancient-vancient-v 2009/12/04 03:12 やっぱりそうですよね・・・投稿しててそんな気がしました。ちょっと週末書き直そうと思います。